MODULE****************************************************************** WARNING: DO NOT EDIT THIS FILE. IT WAS GENERATED MECHANICALLY. See the Makefile for more details. ******************************************************************; R4
*********************************************************** 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; PROCEDUREFromCoords (x0, x1, x2, x3: REAL): T = VAR rr: T; BEGIN rr[0] := x0; rr[1] := x1; rr[2] := x2; rr[3] := x3; RETURN rr END FromCoords; PROCEDUREUnit (i: Axis): T = VAR rr: T; BEGIN rr := Origin; rr[i] := 1.0; RETURN rr END Unit; PROCEDUREEqual (READONLY x, y: T): BOOL = BEGIN RETURN (x[0] = y[0]) AND (x[1] = y[1]) AND (x[2] = y[2]) AND (x[3] = y[3]) END Equal; PROCEDUREIsZero (READONLY x: T): BOOL = BEGIN RETURN (x[0] = 0.0) AND (x[1] = 0.0) AND (x[2] = 0.0) AND (x[3] = 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]; rr[2] := x[2] + y[2]; rr[3] := x[3] + y[3]; 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]; rr[2] := x[2] - y[2]; rr[3] := x[3] - y[3]; RETURN rr END Sub; PROCEDUREMinus (READONLY x: T): T = VAR rr: T; BEGIN rr[0] := - x[0]; rr[1] := - x[1]; rr[2] := - x[2]; rr[3] := - x[3]; 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]; rr[2] := alpha * x[2]; rr[3] := alpha * x[3]; 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; rr[2] := x[2] + delta; rr[3] := x[3] + 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; rr[2] := x[2]*alpha + y[2]*beta; rr[3] := x[3]*alpha + y[3]*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]; rr[2] := x[2] * y[2]; rr[3] := x[3] * y[3]; 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]); rr[2] := F(x[2]); rr[3] := F(x[3]); RETURN rr END FMap; PROCEDURESum (READONLY x: T): REAL = VAR dd: LONGREAL; BEGIN dd := FLOAT(x[0], LONGREAL) + FLOAT(x[1], LONGREAL) + FLOAT(x[2], LONGREAL) + FLOAT(x[3], LONGREAL) ; RETURN FLOAT(dd) END Sum; PROCEDUREMax (READONLY x: T): REAL = BEGIN RETURN MAX(MAX(x[0], x[1]), MAX(x[2], x[3])) END Max; PROCEDUREMin (READONLY x: T): REAL = BEGIN RETURN MIN(MIN(x[0], x[1]), MIN(x[2], x[3])) END Min; PROCEDURESumSq (READONLY x: T): REAL = BEGIN RETURN x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3] END SumSq; PROCEDUREL1Norm (READONLY x: T): REAL = BEGIN RETURN ABS(x[0]) + ABS(x[1]) + ABS(x[2]) + ABS(x[3]) END L1Norm; PROCEDURELInfNorm (READONLY x: T): REAL = BEGIN RETURN MAX(MAX(ABS(x[0]), ABS(x[1])), MAX(ABS(x[2]), ABS(x[3]))) END LInfNorm; PROCEDURELInfDist (READONLY x, y: T): REAL = BEGIN RETURN MAX( MAX(ABS(x[0] - y[0]), ABS(x[1] - y[1])), MAX(ABS(x[2] - y[2]), ABS(x[3] - y[3])) ) END LInfDist; PROCEDUREL1Dist (READONLY x, y: T): REAL = BEGIN RETURN ABS(x[0] - y[0]) + ABS(x[1] - y[1]) + ABS(x[2] - y[2]) + ABS(x[3] - y[3]) END L1Dist; PROCEDUREDist (READONLY x, y: T): REAL = BEGIN RETURN FLOAT(Math.hypot( Math.hypot(FLOAT(x[0] - y[0], LONGREAL), FLOAT(x[1] - y[1], LONGREAL)), Math.hypot(FLOAT(x[2] - y[2], LONGREAL), FLOAT(x[3] - y[3], LONGREAL)) )) END Dist; PROCEDUREL2Dist (READONLY x, y: T): REAL = BEGIN RETURN FLOAT(Math.hypot( Math.hypot(FLOAT(x[0] - y[0], LONGREAL), FLOAT(x[1] - y[1], LONGREAL)), Math.hypot(FLOAT(x[2] - y[2], LONGREAL), FLOAT(x[3] - y[3], 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; t := x[2] - y[2]; dd := dd + t*t; t := x[3] - y[3]; 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 3 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) + FLOAT(x[2], LONGREAL) * FLOAT(y[2], LONGREAL) + FLOAT(x[3], LONGREAL) * FLOAT(y[3], 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 := LInfNorm(x); my := LInfNorm(y); (* Now collect dot product and lengths (rescaled): *) xx := 0.0d0; yy := 0.0d0; xy := 0.0d0; FOR i := 0 TO 3 DO tx := x[i]/mx; ty := y[i]/my; xx := xx + FLOAT(tx*tx, LONGREAL); yy := yy + FLOAT(ty*ty, LONGREAL); xy := xy + FLOAT(tx, LONGREAL)*FLOAT(ty, LONGREAL); END; 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( Math.hypot(FLOAT(x[0], LONGREAL), FLOAT(x[1], LONGREAL)), Math.hypot(FLOAT(x[2], LONGREAL), FLOAT(x[3], LONGREAL)) )) END Length; PROCEDUREL2Norm (READONLY x: T): REAL = BEGIN RETURN FLOAT(Math.hypot( Math.hypot(FLOAT(x[0], LONGREAL), FLOAT(x[1], LONGREAL)), Math.hypot(FLOAT(x[2], LONGREAL), FLOAT(x[3], 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( Math.hypot(FLOAT(x[0], LONGREAL), FLOAT(x[1], LONGREAL)), Math.hypot(FLOAT(x[2], LONGREAL), FLOAT(x[3], LONGREAL)) )); rr[0] := d*x[0]; rr[1] := d*x[1]; rr[2] := d*x[2]; rr[3] := d*x[3]; RETURN rr END Direction; PROCEDUREDet (READONLY p0, p1, p2, p3: T): REAL = (* !!!!!! Should use double precision !!!!!! *) BEGIN RETURN + (p0[0]*p1[1] - p0[1]*p1[0]) * (p2[2]*p3[3] - p2[3]*p3[2]) - (p0[0]*p1[2] - p0[2]*p1[0]) * (p2[1]*p3[3] - p2[3]*p3[1]) + (p0[0]*p1[3] - p0[3]*p1[0]) * (p2[1]*p3[2] - p2[2]*p3[1]) + (p0[1]*p1[2] - p0[2]*p1[1]) * (p2[0]*p3[3] - p2[3]*p3[0]) - (p0[1]*p1[3] - p0[3]*p1[1]) * (p2[0]*p3[2] - p2[2]*p3[0]) + (p0[2]*p1[3] - p0[3]*p1[2]) * (p2[0]*p3[1] - p2[1]*p3[0]) END Det; PROCEDURECross (READONLY p1, p2, p3: T): T = (* !!!!!! Should use double precision !!!!!! *) VAR rr: T; BEGIN rr[0] := + p1[1] * (p2[2]*p3[3] - p2[3]*p3[2]) - p1[2] * (p2[1]*p3[3] - p2[3]*p3[1]) + p1[3] * (p2[1]*p3[2] - p2[2]*p3[1]); rr[1] := - p1[0] * (p2[2]*p3[3] - p2[3]*p3[2]) + p1[2] * (p2[0]*p3[3] - p2[3]*p3[0]) - p1[3] * (p2[0]*p3[2] - p2[2]*p3[0]); rr[2] := + p1[0] * (p2[1]*p3[3] - p2[3]*p3[1]) - p1[1] * (p2[0]*p3[3] - p2[3]*p3[0]) + p1[3] * (p2[0]*p3[1] - p2[1]*p3[0]); rr[3] := - p1[0] * (p2[1]*p3[2] - p2[2]*p3[1]) + p1[1] * (p2[0]*p3[2] - p2[2]*p3[0]) - p1[2] * (p2[0]*p3[1] - p2[1]*p3[0]); RETURN rr END Cross; PROCEDUREInit () = BEGIN FOR i := 0 TO 3 DO Origin[i] := 0.0; Ones[i] := 1.0 END END Init; BEGIN Init(); END R4.