## arithmetic/src/misc/approximation/Interpolation.mg

```GENERIC MODULE Interpolation(R, RT, V);
```
Arithmetic for Modula-3, see doc for details
```IMPORT Arithmetic AS Arith;

CONST Module = "Interpolation.";

PROCEDURE CheckSizes
(READONLY xa: ARRAY OF R.T; READONLY ya: ARRAY OF V.T; ) =
BEGIN
<* ASSERT NUMBER(xa) = NUMBER(ya),
"The number of interpolation nodes and the number of node values must match." *>
END CheckSizes;

PROCEDURE Linear (READONLY xa: ARRAY OF R.T;  (* interpolation nodes *)
READONLY ya: ARRAY OF V.T;  (* interpolation values *)
x : R.T;                                     ):
V.T =
(* Given an interpolation table with xa input and ya output, do linear
interpolation for x. *)
VAR
n                         := NUMBER(xa);
n1                        := 0;
nn                        := n - 1;
diffbest, diff : R.T;
ndx, ndx1, ndx2: CARDINAL;
x1, x2, x12    : R.T;
y1, y2         : V.T;

BEGIN
CheckSizes(xa, ya);

(*---find the best start point---*)
ndx := n1;                   (* this is arbitrary, but fix the FOR loop
if you change *)
diffbest := ABS(x - xa[ndx]);
FOR i := n1 + 1 TO nn DO
diff := ABS(x - xa[i]);
IF diff < RT.Tiny THEN
(* quick victory *)
RETURN ya[i];
ELSIF diff < diffbest THEN
ndx := i;
diffbest := diff;
END;
END;

(*---find the best partner---*)
IF ndx = n1 THEN
ndx1 := n1;
ndx2 := n1 + 1;
ELSIF ndx = nn THEN
ndx1 := nn - 1;
ndx2 := nn;
ELSIF ABS(x - xa[ndx - 1]) < ABS(x - xa[ndx + 1]) THEN
ndx1 := ndx - 1;
ndx2 := ndx;
ELSE
ndx1 := ndx;
ndx2 := ndx + 1;
END;

(*---compute the y value---*)
x1 := xa[ndx1];
y1 := ya[ndx1];
x2 := xa[ndx2];
y2 := ya[ndx2];
x12 := x1 - x2;
RETURN V.Add(V.Scale(y1, (x - x2) / x12), V.Scale(y2, (x1 - x) / x12));
END Linear;

PROCEDURE Newton (READONLY xa: ARRAY OF R.T;  (* interpolation nodes *)
READONLY ya: ARRAY OF V.T;  (* interpolation values *)
x : R.T;
VAR      dy: V.T;                                     ):
V.T RAISES {Arith.Error} =
(* Given an interpolation table with xa input and ya output, do Newton
access: Give the starting index and the length to be used. *)
<* UNUSED *>
CONST
ftn = Module & "Newton";
VAR
col_n: CARDINAL;
c               := NEW(REF ARRAY OF V.T, NUMBER(xa));
d               := NEW(REF ARRAY OF V.T, NUMBER(xa));
ndx             := LAST(xa); (* get a starter x *)
y    : V.T;

BEGIN
CheckSizes(xa, ya);

VAR
difftmp: R.T;
(*---find starting y---*)
diff := ABS(x - xa[ndx]);  (* and its difference from true x *)
BEGIN
FOR i := 0 TO LAST(xa) DO
difftmp := ABS(x - xa[i]);
IF difftmp < RT.Tiny THEN
y := ya[i];
dy := V.Sub(y, y);
(*dy:=V.NewCompliantZero(y);*)
RETURN y;
ELSIF difftmp < diff THEN (* found a better one *)
ndx := i;
diff := difftmp;
END;
c[i] := ya[i];           (* c and d are 1..NUMBER(xa) *)
END;
END;
d^ := c^;                    (* load d from c, thus from ya *)

y := ya[ndx];                (* use the best ndx to get starting y *)

(*---compute and use c and d---*)
col_n := NUMBER(xa);         (* originally there are n in the col *)
FOR m := 1 TO LAST(xa) DO
DEC(col_n);                (* each col recalc loses 1 cell *)
FOR i := 0 TO col_n - 1 DO
VAR
xi    := xa[i];
xim   := xa[i + m];
den   := xi - xim;
delta := V.Sub(c[i + 1], d[i]);
BEGIN
IF ABS(den) < RT.Tiny THEN
RAISE Arith.Error(NEW(Arith.ErrorDivisionByZero).init());
END;
d[i] := V.Scale(delta, (xim - x) / den);
c[i] := V.Scale(delta, (xi - x) / den);
END;
END;
(*---which correction to use?---*)
IF ndx * 2 >= col_n THEN
(* we are at or below the center, need to move up *)
DEC(ndx);
dy := d[ndx];
ELSE
(* we are above the center, need to move down *)
dy := c[ndx];
(* don't need to adjust ndx, because it is effectively moved down
when we slide the next col up *)
END;
(*---update y---*)
END;
RETURN y;
END Newton;

PROCEDURE CubicHermite
(READONLY xa: ARRAY OF R.T;    (* interpolation nodes *)
READONLY ya: ARRAY OF V.T;    (* interpolation values *)
x : R.T;             (* the function argument *)
): V.T =

READONLY yb: ARRAY [0 .. 2] OF V.T  ):
V.T =
(* for some datatypes no Error can occur *)
VAR
x01 := xb[0] - xb[1];
x12 := xb[1] - xb[2];
x02 := xb[0] - xb[2];
xx0 := x - xb[0];
xx1 := x - xb[1];
xx2 := x - xb[2];
sum := V.Scale(yb[0], xx1 * xx2 / (x01 * x02));
BEGIN
sum := V.Sub(sum, V.Scale(yb[1], xx0 * xx2 / (x01 * x12)));
sum := V.Add(sum, V.Scale(yb[2], xx0 * xx1 / (x12 * x02)));
RETURN sum;

(* probably not very efficient *)
PROCEDURE InterpolateHalf (READONLY xb: ARRAY [0 .. 2] OF R.T;
READONLY yb: ARRAY [0 .. 2] OF V.T  ): V.T =
(* for some datatypes no Error can occur *)
CONST Three = FLOAT(3, R.T);
VAR
x01    := xb[0] - xb[1];
x12    := xb[1] - xb[2];
x02    := xb[0] - xb[2];
xin12  := (x - xb[2]) / x12;
hermy1 := xin12 * xin12 * (Three - R.Two * xin12);
(* p(x[1])=1, p(x'[1])=0, p(x[2])=0, p(x'[2])=0 *)
hermdy1 := xin12 * xin12 * (x - xb[1]);
(* p(x[1])=0, p(x'[1])=1, p(x[2])=0, p(x'[2])=0 *)
sum := V.Scale(yb[1], hermdy1 * (x01 - x12) / (x01 * x12) + hermy1);
BEGIN
sum := V.Add(sum, V.Scale(yb[0], hermdy1 * x12 / (x01 * x02)));
sum := V.Sub(sum, V.Scale(yb[2], hermdy1 * x01 / (x12 * x02)));
RETURN sum;
END InterpolateHalf;

PROCEDURE InterpolatePiece (READONLY xb: ARRAY [0 .. 3] OF R.T;
READONLY yb: ARRAY [0 .. 3] OF V.T  ): V.T =
BEGIN
RETURN V.Add(InterpolateHalf(SUBARRAY(xb, 0, 3), SUBARRAY(yb, 0, 3)),
InterpolateHalf(ARRAY OF R.T{xb[3], xb[2], xb[1]},
ARRAY OF V.T{yb[3], yb[2], yb[1]}));
END InterpolatePiece;

BEGIN
CheckSizes(xa, ya);

IF x <= xa[1] THEN
RETURN InterpolateQuadratic(SUBARRAY(xa, 0, 3), SUBARRAY(ya, 0, 3));
ELSE
FOR j := 2 TO LAST(xa) - 1 DO
IF (* xa[j-1]<x AND *) x <= xa[j] THEN
RETURN InterpolatePiece(
SUBARRAY(xa, j - 2, 4), SUBARRAY(ya, j - 2, 4));
END;
END;
RETURN InterpolateQuadratic(SUBARRAY(xa, LAST(xa) - 2, 3),
SUBARRAY(ya, LAST(xa) - 2, 3));
END;
END CubicHermite;

BEGIN
END Interpolation.
```
```

```