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