read recden;
read PGCD;
read AGCD;
read NGCD;

m1 := x^4-10*x^2+1;
m1 := x^8-40*x^6+352*x^4-960*x^2+576;
alias( alpha = RootOf(m1,x) );
m2 := y^3-11*y-13;
alias( beta = RootOf(m2,y) );
g := z^2+12*beta*z+alpha*z/3+31*alpha^3-19;
g := z^2+123*beta*z+alpha*z/13+531*alpha^3-199;
c1 := z^2+alpha*z/12+23*beta-25;#*alpha^3+25;
c2 := z^2+beta/21+23*alpha*z+17;#*alpha^3-17;
g := convert(g,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
c1 := convert(c1,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
c2 := convert(c2,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
N := 5;
f1 := powrpoly(c1,N):
f2 := powrpoly(c2,N):
for k from 1 to N do
lprint(expanding);
  f1 := mulrpoly(quorpoly(f1,c1),g):
  f2 := mulrpoly(quorpoly(f2,c2),g):
lprint(gcding);
  st := time(); G := NGCD(f2,f1): print(k,time()-st);
od:

# Test AGCD
m1 := x^8-40*x^6+352*x^4-960*x^2+576;
alias( alpha = RootOf(m1,x) );
m2 := y^3-11*y-13;
alias( beta = RootOf(m2,y) );
g := z^2+1234*beta*z+alpha*z/13+531*alpha^3-1991;
c1 := z^2+alpha*z/12+23*beta-25;#*alpha^3+25;
c2 := z^2+beta/21+23*alpha*z+17;#*alpha^3-17;
g := convert(g,POLYNOMIAL,[z,a,b],[a=alpha,b=beta]);
c1 := convert(c1,POLYNOMIAL,[z,a,b],[a=alpha,b=beta]);
c2 := convert(c2,POLYNOMIAL,[z,a,b],[a=alpha,b=beta]);
N := 5;
for k from 1 to N do
lprint(expanding);
  f1 := mulrpoly(powrpoly(g,k),powrpoly(c1,N-k)):
  f2 := mulrpoly(powrpoly(g,k),powrpoly(c2,N-k)):
lprint(gcding);
  st := time(); G := NGCD(f1,f2): print(k,time()-st);
od:

# Test integer reconstruction
m1 := x^8-40*x^6+352*x^4-960*x^2+576;
alias( alpha = RootOf(m1,x) );
m2 := y^3-11*y-13;
alias( beta = RootOf(m2,y) );
g := z^2+1234*beta*z+13*alpha*z+531*alpha^3-1991;
c1 := z^2+alpha*z/12+23*beta-25;#*alpha^3+25;
c2 := z^2+beta/21+23*alpha*z+17;#*alpha^3-17;
g := convert(g,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
c1 := convert(c1,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
c2 := convert(c2,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
N := 5;
for k from 1 to N do
lprint(expanding);
  f1 := mulrpoly(powrpoly(g,k),powrpoly(c1,N-k)):
  f2 := mulrpoly(powrpoly(g,k),powrpoly(c2,N-k)):
lprint(gcding);
  st := time(); G := NGCD(f1,f2): print(k,time()-st);
od:


# Test reducible fields
m1 := x^8-40*x^6+352*x^4-960*x^2+576:
alias( alpha = RootOf(m1,x) ):
m2 := expand( (y+1234/5671)*(y^5-2341/1311) ):
alias( beta = RootOf(m2,y) ):
c1 := (alpha*5671*beta+alpha*1234)*z^2+alpha+1:
c2 := beta*z^3+alpha*z+1:
f1 := convert(c1,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
f2 := convert(c2,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
NGCD(f1,f2);
lasterror;
# Repeat with the order of the field extensions reversed.
f1 := convert(c1,POLYNOMIAL,[z,a,b],[a=alpha,b=beta]);
f2 := convert(c2,POLYNOMIAL,[z,a,b],[a=alpha,b=beta]);
NGCD(f1,f2);
lasterror;
# Check that it works if zero divisor not hit
g := z^2+12*beta*z+alpha*z/13+51*alpha^3-19;
c1 := z^2+alpha*z/12+23*beta-25;#*alpha^3+25;
c2 := z^2+beta/21+23*alpha*z+17;#*alpha^3-17;
g := convert(g,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
c1 := convert(c1,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
c2 := convert(c2,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
f1 := mulrpoly(powrpoly(g,3),powrpoly(c1,2));
f2 := mulrpoly(powrpoly(g,3),powrpoly(c2,2));
G := NGCD(f1,f2):

# Test reducible fields with AGCD in the middle
# AGCD needs to recover the zero divsor over the right ring
# It does this by trapping an error and calling PGCD
m1 := y^2+1;
alias( alpha = RootOf(m1,y) ):
m2 := expand( (x-3/2)*(x^4-10*x^2+1) );
alias( beta = RootOf(m2,x) ):
c1 := (2*alpha*beta-3*alpha)*z^2+alpha+1:
c2 := beta*z^3+alpha*z+1:
f1 := convert(c1,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
f2 := convert(c2,POLYNOMIAL,[z,b,a],[b=beta,a=alpha]);
NGCD(f1,f2);
lasterror;


# Test the non-separable case where m(y) = y^3
# Don't try to split y^3-0 because the roots are not separable
g := rpoly( x^3+(y-2)*x^2+1, [x,y], y^3 );
f := rpoly( x^2+(y+31)*x+1, [x,y], y^3 );
h := rpoly( x^3-(11-y)*x^3+31, [x,y], y^3 );
F := mulrpoly(f,g);
H := mulrpoly(h,g);
gcdrpoly(F,H);
NGCD(F,H); # should not try to split y^3