unit R3;

  { Linear algebra in R^3. }
  { Last edited by J. Stolfi on 93-04-20. }

{$N+}
interface

  type T = array [0..2] 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: T; var res: T);
    { Cross-product: res := p x q x r. }

  procedure orthize(x, u: T; var res: T);
    { Returns in res the component of x that is orthogonal to u. }

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

implementation

  procedure add(x, y: T; var res: T);
    var i: integer;
    begin
      for i := 0 to 2 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 2 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 2 do if abs(x[i]) > m then m := abs(x[i]);
      for i := 0 to 2 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 2 do m := m + x[i]*x[i];
      m := sqrt(m);
      for i := 0 to 2 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]
    end;

  procedure cross(p, q: T; var res: T);
    var d01, d02, d12: 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];

      res[0] := + d12;
      res[1] := - d02;
      res[2] := + d01;
    end;

  procedure orthize(x, u: T; var res: T);
    var i: integer;
        xi, ui, sxu, suu, c: double;
    begin
      suu := 0.0;
      sxu := 0.0;
      for i := 0 to 2 do
        begin
          xi := x[i]; ui := u[i];
          suu := suu + ui*ui;
          sxu := sxu + xi*ui;
        end;
      c := sxu/suu;
      for i := 0 to 2 do
        begin
          xi := x[i]; ui := u[i];
          res[i] := xi - c*ui;
        end;
    end;

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

end.