Critical Mass Modula-3: cube/src/R4.m3

cube/src/R4.m3

Copyright (C) 1992, Digital Equipment Corporation All rights reserved. See the file COPYRIGHT for a full description. Created on Sep 15 1988 by Jorge Stolfi Last modified on Wed Dec 22 17:07:55 PST 1993 by mhb modified on Tue Jun 16 18:29:58 PDT 1992 by muller modified on Fri Nov 22 20:20:30 PST 1991 by stolfi modified on Wed Jan 3 03:13:11 1990 by harrison

MODULE R4;

****************************************************************** WARNING: DO NOT EDIT THIS FILE. IT WAS GENERATED MECHANICALLY. See the Makefile for more details. ******************************************************************

*********************************************************** 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;

PROCEDURE FromCoords(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;

PROCEDURE Unit(i: Axis): T =
  VAR rr: T;
  BEGIN
    rr := Origin;
    rr[i] := 1.0;
    RETURN rr
  END Unit;

PROCEDURE Equal(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;

PROCEDURE IsZero(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;

PROCEDURE Add(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;

PROCEDURE Sub(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;

PROCEDURE Minus(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;

PROCEDURE Scale(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;

PROCEDURE Shift(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;

PROCEDURE Mix(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;

PROCEDURE Weigh(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;

PROCEDURE FMap(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;

PROCEDURE Sum (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;

PROCEDURE Max(READONLY x: T): REAL =

  BEGIN

    RETURN MAX(MAX(x[0], x[1]), MAX(x[2], x[3]))

  END Max;

PROCEDURE Min(READONLY x: T): REAL =

  BEGIN

    RETURN MIN(MIN(x[0], x[1]), MIN(x[2], x[3]))

  END Min;

PROCEDURE SumSq(READONLY x: T): REAL =

  BEGIN

    RETURN
        x[0] * x[0]
      + x[1] * x[1]

      + x[2] * x[2]

      + x[3] * x[3]

  END SumSq;

PROCEDURE L1Norm(READONLY x: T): REAL =

  BEGIN

    RETURN
        ABS(x[0])
      + ABS(x[1])

      + ABS(x[2])

      + ABS(x[3])

  END L1Norm;

PROCEDURE LInfNorm(READONLY x: T): REAL =

  BEGIN

    RETURN MAX(MAX(ABS(x[0]), ABS(x[1])), MAX(ABS(x[2]), ABS(x[3])))

  END LInfNorm;

PROCEDURE LInfDist(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;

PROCEDURE L1Dist(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;

PROCEDURE Dist(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;

PROCEDURE L2Dist(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;

PROCEDURE L2DistSq(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;

PROCEDURE RelDist(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;

PROCEDURE Dot(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;

PROCEDURE Cos (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;

PROCEDURE Length(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;

PROCEDURE L2Norm(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;

PROCEDURE Direction(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;

PROCEDURE Det(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;

PROCEDURE Cross(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;

PROCEDURE Init() =
  BEGIN
    FOR i := 0 TO 3 DO
      Origin[i] := 0.0;
      Ones[i] := 1.0
    END
  END Init;

BEGIN
  Init();
END R4.