arithmetic/src/algebra/root/RootBasic.mg

```GENERIC MODULE RootBasic(R, P);
```
Arithmetic for Modula-3, see doc for details
```
IMPORT Arithmetic AS Arith;
IMPORT NumberTheory AS NT;

<* UNUSED *>
CONST
Module = "RootBasic.";
```
handles the case that the leading coefficient c is not 1 multiplies all roots with c to fix that, callers must remember this constant for interpreting the power sum
```PROCEDURE ScaleUpRoots (x: T; ): T =
BEGIN
IF R.Equal(x[LAST(x^)], R.One) THEN
RETURN x;
ELSE
VAR
c, pow: R.T;
y           := NEW(T, NUMBER(x^));
BEGIN
c := x[LAST(x^)];
pow := c;
y[LAST(y^)] := R.One;
y[LAST(y^) - 1] := x[LAST(x^) - 1];
FOR j := LAST(x^) - 2 TO 0 BY -1 DO
y[j] := R.Mul(x[j], pow);
pow := R.Mul(pow, c);
END;
RETURN y;
END;
END;
END ScaleUpRoots;

PROCEDURE ScaleUpRootsFull (VAR x: TBody; c: R.T; ) =
VAR pow: R.T;
BEGIN
pow := c;
FOR j := LAST(x) - 1 TO 0 BY -1 DO
x[j] := R.Mul(x[j], pow);
pow := R.Mul(pow, c);
END;
END ScaleUpRootsFull;
```
reverse ScaleUpRoots PROCEDURE ScaleDownRoots(x:T; c:R.T;):T= <*FATAL Arith.Error*> (* Err.indivisible should not occur, if x and c are chosen properly
```BEGIN
IF R.Equal(c,R.One) THEN
RETURN x;
ELSE
VAR
pow:R.T;
y := NEW(T,NUMBER(x^));
BEGIN
pow:=c;
y[LAST(y^)  ]:=c;
y[LAST(y^)-1]:=x[LAST(x^)-1];
FOR j:=LAST(x^)-2 TO 0 BY -1 DO
y[j]:=R.Div(x[j],pow);
pow:=R.Mul(pow,c);
END;
RETURN y;
END;
END;
END ScaleDownRoots;
*)
```
reverse ScaleUpRoots divide all roots by c and scale the coefficients by c^(n-exp) to achieve fraction free coefficients
```PROCEDURE ScaleDownRoots (VAR x: TBody; c: R.T; exp: INTEGER; ) =
<* FATAL Arith.Error *>        (* ErrorIndivisible should not occur, if x
and c are chosen properly *)
BEGIN
IF NOT R.Equal(c, R.One) THEN
VAR pow: R.T;
BEGIN
pow := c;
FOR j := exp + 1 TO LAST(x) DO
x[j] := R.Mul(x[j], pow);
pow := R.Mul(pow, c);
END;
pow := c;
FOR j := exp - 1 TO 0 BY -1 DO
x[j] := R.Div(x[j], pow);
pow := R.Mul(pow, c);
END;
END;
END;
END ScaleDownRoots;

PROCEDURE FromCardinal (n: CARDINAL; ): R.T =
VAR z: R.T;
BEGIN
z := R.Zero;
FOR j := 1 TO n DO z := R.Add(z, R.One); END;
RETURN z;
END FromCardinal;

PROCEDURE Add (x, y: T; ): T =
<* FATAL Arith.Error *>        (* 'indivisible' cannot occur *)
VAR
qx, qy, qz: T;               (* polynomials with scaled roots,
necessary for eliminating non-one
px, py, pz             : REF PowerSumSeq;
cx, cy                 : R.T;
nrx, nry               : R.T;             (* number of roots *)
binom, rn, rk, rnp, sum: R.T;
nx                     : CARDINAL        := LAST(x^);
ny                     : CARDINAL        := LAST(y^);
nz                     : CARDINAL        := nx * ny;
BEGIN
cx := x[nx];
cy := y[ny];
qx := ScaleUpRoots(x);
qy := ScaleUpRoots(y);
ScaleUpRootsFull(qx^, cy);
ScaleUpRootsFull(qy^, cx);
px := ToPowerSumSeq(qx, nz);
py := ToPowerSumSeq(qy, nz);
pz := NEW(REF PowerSumSeq, nz);
nrx := FromCardinal(nx);
nry := FromCardinal(ny);
rn := R.Zero;
FOR n := 0 TO LAST(pz^) DO
rk := R.Zero;
binom := R.One;
rnp := rn;
(* applying binomial theorem: (x+y)^n=x^n+n*x^(n-1)*y+...+y^n *)
sum := R.Add(R.Mul(nrx, py[n]), R.Mul(nry, px[n]));
FOR k := 1 TO n DO
binom := R.Div(R.Mul(binom, rnp), rk);
rnp := R.Sub(rnp, R.One);
sum := R.Add(sum, R.Mul(binom, R.Mul(px[k - 1], py[n - k])));
END;
pz[n] := sum;
END;
qz := FromPowerSumSeq(pz^);
ScaleDownRoots(qz^, cx, nz - ny);
ScaleDownRoots(qz^, cy, nz - nx);
RETURN qz;

PROCEDURE Sub (x, y: T; ): T =
BEGIN
END Sub;

PROCEDURE Neg (x: T; ): T =
VAR y := NEW(T, NUMBER(x^));
BEGIN
FOR i := LAST(x^) - 1 TO 0 BY -2 DO
y[i + 1] := x[i + 1];
y[i] := R.Neg(x[i]);
END;
IF NUMBER(x^) MOD 2 # 0 THEN y[0] := x[0]; END;
RETURN y;
END Neg;

PROCEDURE IsZero (x: T; ): BOOLEAN =
VAR nonzerofound := FALSE;
BEGIN
IF x = NIL OR NUMBER(x^) <= 1 THEN RETURN FALSE; END;
FOR j := 0 TO LAST(x^) DO
IF NOT R.IsZero(x[j]) THEN
IF nonzerofound THEN
(* two non-zero coefficients found, there must be roots different
from zero *)
RETURN FALSE;
ELSE
nonzerofound := TRUE;
END;
END;
END;
(* at least one non-zero coefficient must found, if we arrive here it
was exactly one no-zero coefficient *)
RETURN nonzerofound;
END IsZero;
```
the temporary values can become very large compared with the resulting values
```PROCEDURE Mul (x, y: T; ): T =
<* FATAL Arith.Error *>        (* 'indivisible' cannot occur *)
VAR
qx, qy: T;                   (* polynomials with scaled roots,
necessary for eliminating non-one
px, py: REF PowerSumSeq;
cx, cy: R.T;
nx    : CARDINAL        := LAST(x^);
ny    : CARDINAL        := LAST(y^);
nz    : CARDINAL        := nx * ny;
z     : T;
BEGIN
cx := x[LAST(x^)];
cy := y[LAST(y^)];
qx := ScaleUpRoots(x);
qy := ScaleUpRoots(y);
px := ToPowerSumSeq(qx, nz);
py := ToPowerSumSeq(qy, nz);
(* assume that py is an independent buffer which can be messed here *)
FOR j := 0 TO LAST(py^) DO py[j] := R.Mul(px[j], py[j]); END;
z := FromPowerSumSeq(py^);
(* more factors might be canceled in some cases (e.g.  1 as zero?)
ScaleDownRoots(z^,cx,nx+2*ny-4); ScaleDownRoots(z^,cy,ny+2*nx-4); *)
ScaleDownRoots(z^, cx, nz - ny);
ScaleDownRoots(z^, cy, nz - nx);
RETURN z;                    (* for testing *)
END Mul;

PROCEDURE Div (x, y: T; ): T RAISES {Arith.Error} =
BEGIN
RETURN Mul(x, Rec(y));
END Div;

PROCEDURE Rec (READONLY x: T; ): T RAISES {Arith.Error} =
VAR y := NEW(T, NUMBER(x^));
BEGIN
IF R.IsZero(x[0]) THEN
RAISE Arith.Error(NEW(Arith.ErrorDivisionByZero).init())
END;
FOR j := 0 TO LAST(x^) DO y[j] := x[LAST(x^) - j]; END;
RETURN y;
END Rec;

PROCEDURE DivMod (x, y: T; VAR r: T; ): T RAISES {Arith.Error} =
BEGIN
r := Zero;
RETURN Mul(x, Rec(y));
END DivMod;

PROCEDURE Mod (<* UNUSED *> x: T; y: T; ): T RAISES {Arith.Error} =
BEGIN
IF IsZero(y) THEN
RAISE Arith.Error(NEW(Arith.ErrorDivisionByZero).init())
END;
RETURN Zero;
END Mod;
```
simple and inefficient method: successively multiply linear factors
```PROCEDURE FromRoots (READONLY root: RootArray; ): T =
VAR
z   := P.One;
fac := P.New(1);
BEGIN
fac[1] := R.One;
FOR j := 0 TO LAST(root) DO
fac[0] := R.Neg(root[j]);
z := P.Mul(z, fac);
END;
RETURN z;
END FromRoots;
```

(* find common roots

```PROCEDURE GCD(x,y:T;):T=
VAR
z:T;
BEGIN
WHILE NOT IsZero(y) DO
```

Reduce should work like Mod, but is allowed to multiply the dividend x by a 'scalar' But do we really get a greatest common divisor by this? E.g. is GCD(3x+2,2x+1) equal to 1 or to 3 ? z:=P.Reduce(x,y);

```    x:=y;
y:=z;
END;
RETURN x;
END GCD;

PROCEDURE ElimMultRoots(x:T;):T=
BEGIN
(* we need a special GCD for this purpose *)
RETURN (GCD(x,P.Derive(x)));
END ElimMultRoots;
*)
```
given the sequence x of power sums, return the next one with respect to the polynomial y consider x as a cyclic buffer for successive sums of powers of the roots
```PROCEDURE GetNextPowerSum (READONLY x: PowerSumSeq; k: CARDINAL; y: T; ):
R.T =
VAR z: R.T;
BEGIN
z := R.Zero;
FOR j := LAST(x) TO 0 BY -1 DO
IF k = 0 THEN k := LAST(x); ELSE DEC(k); END;
END;
RETURN R.Neg(z);
END GetNextPowerSum;

PROCEDURE PowNSlow (x: T; y: CARDINAL; ): T =
(* select each n-th element from the power sum sequence *)
VAR
ps, psz      : REF PowerSumSeq;
j, k, l, m   : CARDINAL;
cp, cx, powcx: R.T;
qx, qz       : T;
<* FATAL Arith.Error *>        (*'indivisible' cannot occur *)
BEGIN
IF y = 0 THEN
RETURN One;
ELSE
cx := x[LAST(x^)];
qx := ScaleUpRoots(x);

ps := ToPowerSumSeq(qx, LAST(qx^));
psz := NEW(REF PowerSumSeq, NUMBER(ps^));
k := y - 1;
j := 0;
WHILE k < NUMBER(ps^) DO psz[j] := ps[k]; INC(j); INC(k, y); END;
l := 0;
m := NUMBER(ps^);
WHILE j < NUMBER(psz^) DO
WHILE m <= k DO
cp := GetNextPowerSum(ps^, l, qx);
ps[l] := cp;
INC(m);
INC(l);
IF l = NUMBER(ps^) THEN l := 0; END;
END;
psz[j] := cp;
INC(j);
INC(k, y);
END;
qz := FromPowerSumSeq(psz^);

(* powcx:=cx^n *)
powcx := cx;
FOR i := 2 TO y DO powcx := R.Mul(powcx, cx); END;
ScaleDownRoots(qz^, powcx, LAST(qz^) - 1);
RETURN qz;
END;
END PowNSlow;
```
speed up by factorizing the exponent into primes
```PROCEDURE PowN (x: T; y: CARDINAL; ): T =
VAR pr: NT.Array;
BEGIN
pr := NT.Factor(y);
FOR j := 0 TO LAST(pr^) DO x := PowNSlow(x, pr[j]); END;
RETURN x;
END PowN;
```

s[k] - the elementary symmetric polynomial of degree k p[k] - sum of the k-th power

s(t) := s[0] + s[1]*t + s[2]*t^2 + ... p(t) := p[1]*t + p[2]*t^2 + ...

Then it holds t*s'(t) + p(-t)*s(t) = 0 This can be proven by considering p as sum of geometric series and differentiating s in the root-wise factored form.

The differential equation can be separated for each power of t which leads to a recurrence equation: n*s[n] + Sum((-1)^j*p[j]*s[n-j],j from 1 to n) = 0

This is known as the Newton-Girard formula. http://mathworld.wolfram.com/Newton-GirardFormulas.html

Note that we index the coefficients the other way round and that the coefficients of the polynomial are not pure elementary symmetric polynomials of the roots but have alternating signs, too.

```PROCEDURE ToPowerSumSeq (x: T; m: CARDINAL; ): REF PowerSumSeq =
VAR
y        := NEW(REF PowerSumSeq, m);
sum: R.T;
div: R.T;
BEGIN
<* ASSERT R.Equal(x[LAST(x^)], R.One) *>
x[LAST(x^)] := R.One;
div := R.One;
FOR n := 1 TO MIN(m, LAST(x^)) DO
sum := R.Mul(x[LAST(x^) - n], div);
FOR j := 1 TO n - 1 DO
sum := R.Add(sum, R.Mul(y[j - 1], x[LAST(x^) - n + j]));
END;
y[n - 1] := R.Neg(sum);    (* R.Div(sum,x[LAST(x^)]), but it is only
divisible if the leading coefficient is
one, what we asserted at the
beginning *)
END;
IF m > LAST(x^) THEN
FOR n := LAST(x^) TO LAST(y^) DO
y[n] :=
GetNextPowerSum(SUBARRAY(y^, n - LAST(x^), LAST(x^)), 0, x);
END;
END;
RETURN y;
END ToPowerSumSeq;

PROCEDURE FromPowerSumSeq (READONLY x: PowerSumSeq; ): T
RAISES {Arith.Error} =
VAR
y        := NEW(T, NUMBER(x) + 1);
sum: R.T;
div: R.T;
BEGIN
y[LAST(y^)] := R.One;
div := R.One;
FOR n := 1 TO LAST(y^) DO
sum := R.Zero;
FOR j := 1 TO n DO
sum := R.Add(sum, R.Mul(x[j - 1], y[LAST(y^) - n + j]));
END;
y[LAST(y^) - n] := R.Neg(R.Div(sum, div));
END;
RETURN y;
END FromPowerSumSeq;
```

PROCEDURE Deflate( x:T; c:R.T; VAR rem:R.T)= VAR xl:=LAST(x^); b,psave:R.T; BEGIN b:=x[xl]; psave:=x[xl-1]; x[xl-1]:=b; FOR i:=xl-2 TO 1 BY -1 DO b:=psave+c*b; psave:=x[i]; x[i]:=b; END; rem:=x[0]+c*x[1]; END Deflate;

```BEGIN
Zero := NEW(T, 2);
Zero^ := TBody{R.Zero, R.One};
One := NEW(T, 2);
One^ := TBody{R.MinusOne, R.One};
END RootBasic.
```
```

```