read recden;
read modulargcd4;
read hensel4;



# Example where cgcd that will encounter a zero-divisor over QQ.
T := [expand( (z+1)*(z+2) )];
Z := [x,z]:
a := x^4 + x^3 + (z+3)*x^2 + (z+4)*x + 3*z+1;
b := x^2 + x + z;
rT := [rpoly(T[1],[z])]:
A := rpoly(a,Z,T);
B := rpoly(b,Z,T);
g := cgcd(A,B,rT);



# Example where cgcd will enounter a zero-divisor in ZZ/11 but not over QQ.
T := [expand( (z+12)*(z+2) )];
Z := [x,z];
a := x^4 + x^3 + (z+3)*x^2 + (z+4)*x + 3*z+1;
b := x^2 + x + z;
rT := [rpoly(T[1],[z])];
A := rpoly(a,Z,T);
B := rpoly(b,Z,T);
g := cgcd(A,B,rT);



# Example where cgcd that will encounter a zero-divisor with fractions.
Z := [x,z]:
a := x^2+z;
b := x^2-(z^2-2/7)*x+1;
# RECDEN represents T using primitive associates over ZZ
T := [numer(expand( (z^2-2/7)*(z+2/3) ))];
rT := [rpoly(T[1],[z])];
A := rpoly(a,Z,T);
B := rpoly(b,Z,T);
g := cgcd(A,B,rT);


# Example where cgcd that will encounter an unlucky zero-divisor mod 11.
t1 := z^4+z+7;
Factor(t1) mod 11; # note z+5 factor
T := [t1];
Z := [x,z]:
a := x^3+z*x;
b := x^3-(z+5+11)*x^2+x; # note z+5 factor mod 11
rT := [rpoly(T[1],[z])];
A := rpoly(a,Z,T);
B := rpoly(b,Z,T);
g := cgcd(A,B,rT);


# Example where cgcd that will require multiple splits.
Z := [x,z]:
m := expand( (z^2+1)*(z-1)*(z-2) );
a := x^3+z*x;
b := x^3-(z^2+1)*x^2+2*x+4;
T := [m]:
rT := [rpoly(T[1],[z])];
A := rpoly(a,Z,T);
B := rpoly(b,Z,T);
g := cgcd(A,B,rT);


# Example with multiple extensions and using 31 bit primes
t := [z^3 - 20*z^2 - 4*z - 89, 
      y^3 + 28*y^2 + (-42*z+69)*y - 33*z^2 + 80*z - 77];
T := construct_triangular_set(t, [y,z]):
X := [x,y,z];
a := x^3 + 89*x*y - 16*x*z + 59*y^2 - 69*z^2 + 30*y - 64*z;
b := x^2 - 46*x - 33*y + 87*z - 34;
g := randpoly(X,dense,degree=3,coeffs=rand(10^10));
A := rpoly(a*g, X):
B := rpoly(b*g, X):
A := reduce_by_tset(A, T):
B := reduce_by_tset(B, T):
G := cgcd(A,B,T,2^31);
divrpoly(A,G);
divrpoly(B,G);


# Example where cgcd that will encounter bad and unlucky primes.
T := [numer(expand( (z^2-2/17)*(z+2/3) ))];
Z := [x,z]:
a := (26*z+13)*x^2+z+1;
b := (26*z+13)*x^2+z+24;
g := x^2+z+1/3;
# primes 3,13,17 are bad and 19 is unlucky
rT := [rpoly(T[1],[z])];
A := rpoly(a,Z,T):
B := rpoly(b,Z,T):
G := rpoly(g,Z,T):
A := mulrpoly(A,G):
B := mulrpoly(B,G):
g := cgcd(A,B,rT);


# Hensel lifting example over an extension. The coefficients
# of the factors have denominators even though their product
# doesn't.

t1 := z^6+3*z^5+6*z^4+z^3-3*z^2+12*z+16;
f := x^3 - 3;
a := -1/12*z^5 - 1/4*z^4 - 1/2*z^3 - 5/12*z^2 + 1/4*z + x - 1;
b := -1/12*z^5 - 1/12*z^4 - 1/6*z^3 + 7/12*z^2 - 11/12*z + x - 4/3;
c := 1/6*z^5 + 1/3*z^4 + 2/3*z^3 - 1/6*z^2 + 2/3*z + x + 7/3; 
p := 13;
A := rpoly(a, [x,z], t1, p);
B := rpoly(b*c, [x,z], t1, p);
F := rpoly(f,[x,z], t1):
T := [rpoly(t1, [z])]:
mulrpoly(A,B); # Should be x^3 + 10 mod <t1, 13>.
hensellift_univar(A,B,F,p);