PROCEDURE Cof(m: T; ix, jx: CARDINAL): LONGREAL =
  VAR r, d: LONGREAL;
      t: R2x2.T;
      ii, jj: CARDINAL;
  BEGIN
    (* Extract minor *)
    ii := 0;
    FOR i := 0 TO 2 DO
      IF i # ix THEN
        jj := 0;
        FOR j := 0 TO 2 DO
          IF j # jx THEN
            t[ii,jj] := m[i,j];
            INC(jj)
          END
        END;
        INC(ii)
      END
    END;
    (* Compute its determinant *)
    d := R2x2.Det(t);
    IF (ix + jx) MOD 2 = 1 THEN d := -d END;
    RETURN r
  END Cof;
    
PROCEDURE AdjAux(READONLY m: T; adj, det: BOOLEAN; VAR a: T; VAR d: LONGREAL) =
  (* 
    If "adj" is true, returns the adjoint of "m" in "a".
    If "det" is true, returns the determinant of "m" in "d". *)

  BEGIN
    WITH
      m_00 = FLOAT(m[0,0], LONGREAL),
      m_01 = FLOAT(m[0,1], LONGREAL),
      m_02 = FLOAT(m[0,2], LONGREAL),
      
      m_10 = FLOAT(m[1,0], LONGREAL),
      m_11 = FLOAT(m[1,1], LONGREAL),
      m_12 = FLOAT(m[1,2], LONGREAL),
      
      m_20 = FLOAT(m[2,0], LONGREAL),
      m_21 = FLOAT(m[2,1], LONGREAL),
      m_22 = FLOAT(m[2,2], LONGREAL),
      
      m_01_01 = m_00 * m_11 - m_01 * m_10,
      m_01_02 = m_00 * m_12 - m_02 * m_10,
      m_01_12 = m_01 * m_12 - m_02 * m_11
    DO
      IF adj THEN
        WITH
          m_02_01 = m_00 * m_21 - m_01 * m_20,
          m_02_02 = m_00 * m_22 - m_02 * m_20,
          m_02_12 = m_01 * m_22 - m_02 * m_21,

          m_12_01 = m_10 * m_21 - m_11 * m_20,
          m_12_02 = m_10 * m_22 - m_12 * m_20,
          m_12_12 = m_11 * m_22 - m_12 * m_21
        DO
          a[0,0] := FLOAT(+ m_12_12);
          a[1,0] := FLOAT(- m_12_02);
          a[2,0] := FLOAT(+ m_12_01);

          a[0,1] := FLOAT(- m_02_12);
          a[1,1] := FLOAT(+ m_02_02);
          a[2,1] := FLOAT(- m_02_01);

          a[0,2] := FLOAT(+ m_01_12);
          a[1,2] := FLOAT(- m_01_02);
          a[2,2] := FLOAT(+ m_01_01);
        END;
      END;
      IF det THEN
        d := + m_01_01 * m_22 - m_01_02 * m_21 + m_01_12 * m_20
      END;
    END
  END AdjAux;

(* *)
PROCEDURE dst( frm,tto: Point ): REAL =
VAR r: REAL;
    dd, vi, fw, tw: REAL;
BEGIN
    dd := 0.0;
    fw := frm.c[0];
    tw := tto.c[0];
    FOR i:= 1 TO 3 DO
       vi := fw * tto.c[i] - tw * frm.c[i];
       dd := dd + vi * vi;
    END; 
    r := FLOAT( Math.sqrt( FLOAT( dd, LONGREAL ) ) ) / fw / tw;
    RETURN ABS(r);
END dst;

PROCEDURE Dist(READONLY a, b: Point): REAL =
  BEGIN
    WITH 
      aw = FLOAT(a.c[0], LONGREAL),
      bw = FLOAT(b.c[0], LONGREAL),
      
      ax = FLOAT(a.c[1], LONGREAL) / aw,
      bx = FLOAT(b.c[1], LONGREAL) / bw,
      dx = ax - bx,
      
      ay = FLOAT(a.c[2], LONGREAL) / aw,
      by = FLOAT(b.c[2], LONGREAL) / bw,
      dy = ay - by,
      
      az = FLOAT(a.c[3], LONGREAL) / aw,
      bz = FLOAT(b.c[3], LONGREAL) / bw,
      dz = az - bz,
      
     length = Math.sqrt(dx * dx + dy * dy + dz * dz)
    DO
      RETURN FLOAT(length)
    END;
  END Dist;    

(* Version of R4.T with "FOR" loops: *)

PROCEDURE Add(x, y: T): T =
  VAR r: T;
  BEGIN
    FOR i := 0 TO 3 DO  r[i] := x[i] + y[i] END;
    RETURN r
  END Add;

PROCEDURE Sub(x, y: T): T =
  VAR r: T;
  BEGIN
    FOR i := 0 TO 3 DO  r[i] := x[i] - y[i] END;
    RETURN r
  END Sub;

PROCEDURE Neg(x: T): T =
  VAR r: T;
  BEGIN
    FOR i := 0 TO 3 DO  r[i] := - x[i] END;
    RETURN r
  END Neg;

PROCEDURE Scale(x: T; s: REAL): T = 
  VAR r: T;
  BEGIN
    FOR i := 0 TO 3 DO r[i] := x[i] * s END;
    RETURN r
  END Scale;

PROCEDURE ReduceInf(x: T): T =
  VAR m: REAL;
      r: T;
  BEGIN
    m := 0.0;
    FOR i := 0 TO 3 DO 
       IF ABS(x[i]) > m THEN m := ABS(x[i]) END;
    END;
    FOR i := 0 TO 3 DO r[i] := x[i] / m END;
    RETURN r
  END ReduceInf;  

PROCEDURE Reduce(x: T): T =
  VAR r: T;
      m: LONGREAL;
  BEGIN
    m := 0.0D0;
    FOR i:= 0 TO 3 DO 
      WITH xi = FLOAT(x[i], LONGREAL) DO m := m + xi * xi END
    END;
    WITH sm = FLOAT(Math.sqrt(m)) DO
      FOR i:= 0 TO 3 DO r[i] := x[i] / m END
    END;
    RETURN r
  END Reduce; 

PROCEDURE Dot(x, y: T): LONGREAL =
  VAR m: LONGREAL;
  BEGIN
    FOR i:= 0 TO 3 DO 
      WITH 
        xi = FLOAT(x[i], LONGREAL),
        yi = FLOAT(x[i], LONGREAL)
      DO
        m := m + xi * yi
      END
    END;
    RETURN m
  END Dot;

PROCEDURE Cross(x, y: T): T =
  VAR r: T;
  BEGIN
    WITH
      x0 = FLOAT(x[0], LONGREAL),
      x1 = FLOAT(x[1], LONGREAL),
      x2 = FLOAT(x[2], LONGREAL),
      x3 = FLOAT(x[3], LONGREAL),
      y0 = FLOAT(y[0], LONGREAL),
      y1 = FLOAT(y[1], LONGREAL),
      y2 = FLOAT(y[2], LONGREAL),
      y3 = FLOAT(y[3], LONGREAL),
      d01 = x0*y1 - x1*y0,
      d02 = x0*y2 - x2*y0,
      d12 = x1*y2 - x2*y1,
      d03 = x0*y3 - x3*y0,
      d13 = x1*y3 - x3*y1,
      d23 = x2*y3 - x3*y2
    DO
      r[0] := FLOAT(- d12*z[3] + d13*z[2] - d23*z[1]);
      r[1] := FLOAT(+ d02*z[3] - d03*z[2] + d23*z[0]);
      r[2] := FLOAT(- d01*z[3] + d03*z[1] - d13*z[0]);
      r[3] := FLOAT(+ d01*z[2] - d02*z[1] + d12*z[0]);
    END;
    RETURN r
  END Cross;

PROCEDURE Orthize(x, u: T): T =
  VAR r: T;
      xi, ui, sxu, suu, c: LONGREAL;
  BEGIN
   suu := 0.0D0;
   sxu := 0.0D0;
   FOR i := 0 TO 3 DO 
     WITH 
       xi = FLOAT(x[i], LONGREAL),
       ui = FLOAT(u[i], LONGREAL)
     DO
       suu := suu + ui*ui;
       sxu := sxu + xi*ui; 
     END
   END;
   c := sxu/suu;
   FOR i := 0 TO 3 DO 
     WITH 
       xi = FLOAT(x[i], LONGREAL),
       ui = FLOAT(u[i], LONGREAL)
     DO
       r[i] := FLOAT(xi - c*ui);
   END;
   RETURN r
  END Orthize;
  
(* End FOR version *)

PROCEDURE Cof(READONLY m: T; ix, jx: CARDINAL): LONGREAL;   
  (*
    Returns the cofactor of element "[ix,jx]" in matrix "m". *)

PROCEDURE Cof(m: T; ix, jx: CARDINAL): LONGREAL =
  VAR t: R3x3.T;
      ii, jj: CARDINAL;
  BEGIN
    (* Extract minor *)
    ii := 0;
    FOR i := 0 TO 3 DO
      IF i # ix THEN
        jj := 0;
        FOR j := 0 TO 3 DO
          IF j # jx THEN
            t[ii,jj] := m[i,j];
            INC(jj)
          END
        END;
        INC(ii)
      END
    END;
    (* Compute its determinant *)
    WITH d = R3x3.Det(t) DO
      IF (ix + jx) MOD 2 = 1 THEN 
        RETURN FLOAT(-d)
      ELSE
        RETURN FLOAT(d)
      END
    END
  END Cof;
    
PROCEDURE Det(READONLY m: T): LONGREAL =
  BEGIN
    WITH
      m_00 = FLOAT(m[0,0], LONGREAL),
      m_01 = FLOAT(m[0,1], LONGREAL),
      m_02 = FLOAT(m[0,2], LONGREAL),
      
      m_10 = FLOAT(m[1,0], LONGREAL),
      m_11 = FLOAT(m[1,1], LONGREAL),
      m_12 = FLOAT(m[1,2], LONGREAL),
      
      m_20 = FLOAT(m[2,0], LONGREAL),
      m_21 = FLOAT(m[2,1], LONGREAL),
      m_22 = FLOAT(m[2,2], LONGREAL),
      
      m_01_01 = m_00 * m_11 - m_01 * m_10,
      m_01_02 = m_00 * m_12 - m_02 * m_10,
      m_01_12 = m_01 * m_12 - m_02 * m_11,

      d = + m_01_01 * m_22 - m_01_02 * m_21 + m_01_12 * m_20

    DO
      RETURN d
    END
  END Det;  

PROCEDURE Det(READONLY m: T): LONGREAL =
  BEGIN
    WITH
      m_00 = FLOAT(m[0,0], LONGREAL),
      m_01 = FLOAT(m[0,1], LONGREAL),
      m_02 = FLOAT(m[0,2], LONGREAL),
      m_03 = FLOAT(m[0,3], LONGREAL),
      
      m_10 = FLOAT(m[1,0], LONGREAL),
      m_11 = FLOAT(m[1,1], LONGREAL),
      m_12 = FLOAT(m[1,2], LONGREAL),
      m_13 = FLOAT(m[1,3], LONGREAL),
      
      m_20 = FLOAT(m[2,0], LONGREAL),
      m_21 = FLOAT(m[2,1], LONGREAL),
      m_22 = FLOAT(m[2,2], LONGREAL),
      m_23 = FLOAT(m[2,3], LONGREAL),
      
      m_30 = FLOAT(m[3,0], LONGREAL),
      m_31 = FLOAT(m[3,1], LONGREAL),
      m_32 = FLOAT(m[3,2], LONGREAL),
      m_33 = FLOAT(m[3,3], LONGREAL),
      
      m_01_01 = m_00 * m_11 - m_01 * m_10,
      m_01_02 = m_00 * m_12 - m_02 * m_10,
      m_01_12 = m_01 * m_12 - m_02 * m_11,
      m_01_03 = m_00 * m_13 - m_03 * m_10,
      m_01_13 = m_01 * m_13 - m_03 * m_11,
      m_01_23 = m_02 * m_13 - m_03 * m_12,

      m_23_01 = m_20 * m_31 - m_21 * m_30,
      m_23_02 = m_20 * m_32 - m_22 * m_30,
      m_23_12 = m_21 * m_32 - m_22 * m_31,
      m_23_03 = m_20 * m_33 - m_23 * m_30,
      m_23_13 = m_21 * m_33 - m_23 * m_31,
      m_23_23 = m_22 * m_33 - m_23 * m_32,
      
      d = + m_01_01 * m_23_23 - m_01_02 * m_23_13 + m_01_12 * m_23_03
          + m_01_03 * m_23_12 - m_01_13 * m_23_02 + m_01_23 * m_23_01
      
    DO
      RETURN d
    END
  END Det;