MODULE****************************************************************** WARNING: DO NOT EDIT THIS FILE. IT WAS GENERATED MECHANICALLY. See the Makefile for more details. ******************************************************************; R2
*********************************************************** Disclaimer: the numerical algorithms were quickly hacked from the Modula-2+ version. They are not suppose to be the best possible, not even close. There are surely gross blunders, especially in the choice of LONGREAL vs REAL for intermediary results. ***********************************************************
IMPORT Math, Random, Fmt; PROCEDUREUnit (i: Axis): T = VAR rr: T; BEGIN rr := Origin; rr[i] := 1.0; RETURN rr END Unit; PROCEDUREEqual (READONLY x, y: T): BOOLEAN = BEGIN RETURN (x[0] = y[0]) AND (x[1] = y[1]) END Equal; PROCEDUREIsZero (READONLY x: T): BOOLEAN = BEGIN RETURN (x[0] = 0.0) AND (x[1] = 0.0) END IsZero; PROCEDUREAdd (READONLY x, y: T): T = VAR rr: T; BEGIN rr[0] := x[0] + y[0]; rr[1] := x[1] + y[1]; RETURN rr END Add; PROCEDURESub (READONLY x, y: T): T = VAR rr: T; BEGIN rr[0] := x[0] - y[0]; rr[1] := x[1] - y[1]; RETURN rr END Sub; PROCEDUREMinus (READONLY x: T): T = VAR rr: T; BEGIN rr[0] := - x[0]; rr[1] := - x[1]; RETURN rr END Minus; PROCEDUREScale (alpha: REAL; READONLY x: T): T = VAR rr: T; BEGIN rr[0] := alpha * x[0]; rr[1] := alpha * x[1]; RETURN rr END Scale; PROCEDUREShift (READONLY x: T; delta: REAL): T = VAR rr: T; BEGIN rr[0] := x[0] + delta; rr[1] := x[1] + delta; RETURN rr END Shift; PROCEDUREMix (READONLY x: T; alpha: REAL; READONLY y: T; beta: REAL): T = VAR rr: T; BEGIN rr[0] := x[0]*alpha + y[0]*beta; rr[1] := x[1]*alpha + y[1]*beta; RETURN rr END Mix; PROCEDUREWeigh (READONLY x, y: T): T = VAR rr: T; BEGIN rr[0] := x[0] * y[0]; rr[1] := x[1] * y[1]; RETURN rr END Weigh; PROCEDUREFMap (READONLY x: T; F: Function): T = VAR rr: T; BEGIN rr[0] := F(x[0]); rr[1] := F(x[1]); RETURN rr END FMap; PROCEDURESum (READONLY x: T): REAL = VAR dd: LONGREAL; BEGIN dd := FLOAT(x[0], LONGREAL) + FLOAT(x[1], LONGREAL) ; RETURN FLOAT(dd) END Sum; PROCEDUREMax (READONLY x: T): REAL = BEGIN RETURN MAX(x[0], x[1]) END Max; PROCEDUREMaxAbsAxis (READONLY x: T): Axis = VAR a: Axis; BEGIN IF ABS(x[0]) > ABS(x[1]) THEN a := 0 ELSE a := 1 END; RETURN a; END MaxAbsAxis; PROCEDUREMin (READONLY x: T): REAL = BEGIN RETURN MIN(x[0], x[1]) END Min; PROCEDURESumSq (READONLY x: T): REAL = BEGIN RETURN x[0] * x[0] + x[1] * x[1] END SumSq; PROCEDUREL1Norm (READONLY x: T): REAL = BEGIN RETURN ABS(x[0]) + ABS(x[1]) END L1Norm; PROCEDURELInfNorm (READONLY x: T): REAL = BEGIN RETURN MAX(ABS(x[0]), ABS(x[1])) END LInfNorm; PROCEDURELInfDist (READONLY x, y: T): REAL = BEGIN RETURN MAX(ABS(x[0] - y[0]), ABS(x[1] - y[1])) END LInfDist; PROCEDUREL1Dist (READONLY x, y: T): REAL = BEGIN RETURN ABS(x[0] - y[0]) + ABS(x[1] - y[1]) END L1Dist; PROCEDUREDist (READONLY x, y: T): REAL = BEGIN RETURN FLOAT(Math.hypot(FLOAT(x[0] - y[0], LONGREAL), FLOAT(x[1] - y[1], LONGREAL))) END Dist; PROCEDUREL2Dist (READONLY x, y: T): REAL = BEGIN RETURN FLOAT(Math.hypot(FLOAT(x[0] - y[0], LONGREAL), FLOAT(x[1] - y[1], LONGREAL))) END L2Dist; PROCEDUREL2DistSq (READONLY x, y: T): REAL = VAR t, dd: REAL; BEGIN t := x[0] - y[0]; dd := t*t; t := x[1] - y[1]; dd := dd + t*t; RETURN dd END L2DistSq; PROCEDURERelDist (READONLY x, y: T; eps: REAL := 1.0e-37): REAL = VAR u, v: REAL; s, m: REAL; BEGIN s := 0.0; FOR i := 0 TO 1 DO u := x[i]; v := y[i]; m := MAX(MAX(ABS(u), ABS(v)), eps); s := MAX(ABS(u/m - v/m) - eps/m, s); END; RETURN s END RelDist; PROCEDUREDot (READONLY x, y: T): REAL = BEGIN RETURN FLOAT( FLOAT(x[0], LONGREAL) * FLOAT(y[0], LONGREAL) + FLOAT(x[1], LONGREAL) * FLOAT(y[1], LONGREAL) ) END Dot; PROCEDURECos (READONLY x, y: T): REAL = VAR xy, xx, yy: LONGREAL; tx, ty, mx, my: REAL; BEGIN (* Compute rescaling factors to avoid spurious overflow: *) mx := MAX(ABS(x[0]), ABS(x[1])); my := MAX(ABS(y[0]), ABS(y[1])); (* Now collect dot product and lengths (rescaled): *) tx := x[0]/mx; ty := y[0]/my; xx := FLOAT(tx*tx, LONGREAL); yy := FLOAT(ty*ty, LONGREAL); xy := FLOAT(tx, LONGREAL)*FLOAT(ty, LONGREAL); tx := x[1]/mx; ty := y[1]/my; xx := xx + FLOAT(tx*tx, LONGREAL); yy := yy + FLOAT(ty*ty, LONGREAL); xy := xy + FLOAT(tx, LONGREAL) * FLOAT(ty, LONGREAL); xx := 1.0d0 / Math.sqrt(FLOAT(xx*yy, LONGREAL)); xx := xx*xy; RETURN FLOAT(xx) END Cos; PROCEDURELength (READONLY x: T): REAL = BEGIN RETURN FLOAT(Math.hypot(FLOAT(x[0], LONGREAL), FLOAT(x[1], LONGREAL))) END Length; PROCEDUREL2Norm (READONLY x: T): REAL = BEGIN RETURN FLOAT(Math.hypot(FLOAT(x[0], LONGREAL), FLOAT(x[1], LONGREAL))) END L2Norm; PROCEDUREDirection (READONLY x: T): T = (* !!!!! Should try to avoid spurious overflow by prescaling x *) VAR d: REAL; rr: T; BEGIN d := 1.0/FLOAT(Math.hypot(FLOAT(x[0], LONGREAL), FLOAT(x[1], LONGREAL))); rr[0] := d*x[0]; rr[1] := d*x[1]; RETURN rr END Direction; PROCEDUREDet (READONLY p0, p1: T): REAL = (* !!!!!! Should use double precision !!!!!! *) BEGIN RETURN p0[0]*p1[1] - p0[1]*p1[0] END Det; PROCEDURECross (READONLY p1: T): T = (* !!!!!! Should use double precision !!!!!! *) VAR rr: T; BEGIN rr[0] := p1[1]; rr[1] := -p1[0]; RETURN rr END Cross; PROCEDUREThrow (lo, hi: REAL; src: Random.T := NIL): T = VAR rr: T; dd, xx: REAL; BEGIN IF src = NIL THEN src := NEW(Random.Default).init() END; dd := hi - lo; <* ASSERT dd > 0.0 *> <* ASSERT dd/MAX(1.0E-25, MAX(ABS(hi), ABS(lo))) > 1.0E-6 *> FOR i := 0 TO 1 DO REPEAT xx := src.real() UNTIL xx > 0.0; rr[i] := lo + dd*xx END; RETURN rr END Throw; PROCEDUREToText (READONLY x: T): TEXT = BEGIN RETURN "(" & Fmt.Real(x[0]) & " " & Fmt.Real(x[1]) & ")"; END ToText; BEGIN END R2.