unit R4;

  { Linear algebra in R^4. }
  { Last edited by J. Stolfi on 93-04-20. }
{$N+}
interface

  type T = array [0..3] of double;

  procedure add(x, y: T; var res: T);
    { Sets res := x + y. }

  procedure scale(x: T; a: double; var res: T);
    { Sets res := a * x. }

  procedure reduce_inf(var x: T);
    { Scales x to unit L_infinity norm. }

  procedure reduce(var x: T);
    { Scales x to unit Euclidean norm. }

  function dot(x, y: T): double;
    { Dot product: sum of x[i]*y[i]. }

  procedure cross(p, q, r: T; var res: T);
    { Cross-product: res := p x q x r. }

  procedure print(var f: TEXT; x: T);
    { Prints x on f. }

implementation

  procedure add(x, y: T; var res: T);
    var i: integer;
    begin
      for i := 0 to 3 do res[i] := x[i] + y[i];
    end;

  procedure scale(x: T; a: double; var res: T);
    var i: integer;
    begin
      for i := 0 to 3 do res[i] := a * x[i];
    end;

  procedure reduce_inf(var x: T);
    var i: integer;
        m: real;
    begin
      m := 0.0;
      for i := 0 to 3 do if abs(x[i]) > m then m := abs(x[i]);
      for i := 0 to 3 do x[i] := x[i] / m;
    end;

  procedure reduce(var x: T);
    var i: integer;
        m: double;
    begin
      m := 0.0;
      for i := 0 to 3 do m := m + x[i]*x[i];
      m := sqrt(m);
      for i := 0 to 3 do x[i] := x[i] / m;
    end;

  function dot(x, y: T): double;
    begin
      dot := x[0]*y[0] + x[1]*y[1] + x[2]*y[2] + x[3]*y[3]
    end;

  procedure cross(p, q, r: T; var res: T);
    var d01, d02, d12, d03, d13, d23: double;
    begin
      d01 := p[0]*q[1] - p[1]*q[0];
      d02 := p[0]*q[2] - p[2]*q[0];
      d12 := p[1]*q[2] - p[2]*q[1];
      d03 := p[0]*q[3] - p[3]*q[0];
      d13 := p[1]*q[3] - p[3]*q[1];
      d23 := p[2]*q[3] - p[3]*q[2];

      res[0] := - d12*r[3] + d13*r[2] - d23*r[1];
      res[1] := + d02*r[3] - d03*r[2] + d23*r[0];
      res[2] := - d01*r[3] + d03*r[1] - d13*r[0];
      res[3] := + d01*r[2] - d02*r[1] + d12*r[0];
    end;

  procedure print(var f: TEXT; x: T);
    var i: integer;
    begin
      for i := 0 to 3 do
        write(f, x[i]:16:8)
    end;

end.