MODULE SPPolyFunction;

IMPORT Wr, Rd, Fmt, Thread;
IMPORT SPFunction, SPHomoLabel;
IMPORT FileFmt, NGet, NPut, FGet, FPut;

TYPE
  Point = SPFunction.Point; 
  Gradient = SPFunction.Gradient;

CONST MaxDegree = 6;

REVEAL
  T = PublicT BRANDED OBJECT
    OVERRIDES
      eval := Eval;
      grad := Grad;
      slap := SLap;
      toMaple := ToMaple;
      read := Read;
      write := Write;
    END;

<* FATAL Wr.Failure, Thread.Alerted *>

CONST FileVersion = "98-04-09";

PROCEDURE Write(t: T; wr: Wr.T) =
  VAR l: CARDINAL;
  BEGIN
    FileFmt.WriteHeader(wr, "SPPolyFunction.T", FileVersion);
    WITH
      c0 = t.c0^,
      c1 = t.c1^,
      deg = t.deg
    DO
      NPut.Int(wr, "degree", deg); FPut.EOL(wr);
      l:= 0;
      FOR j:= 0 TO deg DO 
        FOR k := 0 TO j  DO 
          IF k > 0 THEN FPut.Space(wr) END;
          FPut.LongReal(wr, c0[l]);
          l:= l + 1;
        END;
        FPut.EOL(wr);
      END;
      l:= 0;
      FOR j:= 0 TO deg-1 DO 
        FOR k := 0 TO j  DO 
          IF k > 0 THEN FPut.Space(wr) END;
          FPut.LongReal(wr, c1[l]);
          l:= l + 1;
        END;
        FPut.EOL(wr);
      END;
    END;
    FileFmt.WriteFooter(wr, "SPPolyFunction.T")
  END Write;

PROCEDURE Read(t: T; rd: Rd.T) =
  VAR l: CARDINAL;
  BEGIN
    FileFmt.ReadHeader(rd, "SPPolyFunction.T", FileVersion);
    WITH 
      deg = NGet.Int(rd, "degree"),
      
      nc0 = (deg + 1) * (deg + 2) DIV 2,
      rc0 = NEW(REF Coeffs, nc0), c0 = rc0^,
      
      nc1 = deg * (deg + 1) DIV 2,
      rc1 = NEW(REF Coeffs, nc1), c1 = rc1^
    DO 
      t.deg := deg;
      t.c0 := rc0;
      t.c1 := rc1;
      FGet.EOL(rd);
      l:= 0;
      FOR j:= 0 TO deg DO 
        FOR k := 0 TO deg-j  DO 
          FGet.Skip(rd);
          c0[l] := FGet.LongReal(rd);
          l:= l + 1;
        END;
        FGet.EOL(rd);
      END;
      l:= 0;
      FOR j:= 0 TO deg-1 DO 
        FOR k := 0 TO deg-1-j  DO 
          FGet.Skip(rd);
          c1[l] := FGet.LongReal(rd);
          l:= l + 1;
        END;
        FGet.EOL(rd);
      END;
    END;
    FileFmt.ReadFooter(rd, "SPPolyFunction.T")
  END Read;

PROCEDURE Eval(t: T; READONLY p: Point): LONG = 
  VAR k,r: INTEGER;
      v: LONG;
      xp, yp, zp: ARRAY [0..MaxDegree] OF LONG; (* Coordinate powers *)
  BEGIN
    WITH
      d = t.deg,
      x = p[0], 
      y = p[1], 
      z = p[2],
      c0 = t.c0^,
      c1 = t.c1^
    DO
      (* Compute power tables: *)
      xp[0] := 1.0d0; yp[0] := 1.0d0; zp[0] := 1.0d0;
      IF d > 0 THEN xp[1] := x; yp[1] := y; zp[1] := z END;
      FOR k := 2 TO d DO 
        xp[k] := xp[k-1]*x;
        yp[k] := yp[k-1]*y;
        zp[k] := zp[k-1]*z
      END;
      (* Evaluate polynomial: *)
      v := 0.0d0;
      r := 0;
      FOR m := 0 TO d  DO 
        FOR  n := 0 TO d - m DO 
          k := d - m - n;
          v := v + c0[r] * xp[m] * yp[n] * zp[k];
          r := r + 1;
        END;
      END;
      r := 0;
      FOR m := 0 TO d - 1  DO 
        FOR  n := 0 TO d - 1 - m DO 
          k := d - 1 - m - n;
          v := v + c1[r] * xp[m] * yp[n] * zp[k];
          r := r + 1;
        END;
      END;
    END;
    RETURN v
  END Eval;

PROCEDURE Grad(<*UNUSED*> t: T; <*UNUSED*> READONLY p: Point): Gradient = 
  BEGIN
    <* ASSERT FALSE *>
  END Grad;

PROCEDURE SLap(<*UNUSED*> t: T; <*UNUSED*> READONLY p: Point): LONGREAL = 
  BEGIN
    <* ASSERT FALSE *>
  END SLap;

PROCEDURE ToMaple(t: T; f: Wr.T) =
  BEGIN
    ToMapleHomo(f, t.deg, t.c0^);
    Wr.PutText(f, " + ");
    ToMapleHomo(f, t.deg-1, t.c1^);
    Wr.Flush(f)
  END ToMaple;
 
PROCEDURE ToMapleHomo(f: Wr.T; deg: CARDINAL; READONLY c: Coeffs) =
  BEGIN
    Wr.PutText(f, "\n# Warning: coeff labels may be wrong for poly!\n");
    FOR l := 0 TO LAST(c) DO
      WITH
        label = SPHomoLabel.IndexToLabel(l, deg)
      DO
        IF l > 0 THEN 
          Wr.PutText(f, " + ");
          IF l MOD 5 = 0 THEN Wr.PutText(f, "\n") END;
        END; 
        Wr.PutText(f, Fmt.LongReal(c[0])); 
        Wr.PutText(f, "*x^"); Wr.PutText(f, Fmt.Int(label[2])); 
        Wr.PutText(f, "*y^"); Wr.PutText(f, Fmt.Int(label[1])); 
        Wr.PutText(f, "*z^"); Wr.PutText(f, Fmt.Int(label[0])); 
      END
    END
  END ToMapleHomo;

BEGIN
END SPPolyFunction.