restart;
with(PolynomialTools):
currentdir("/Users/mmonagan/Desktop");
read "recden";
#The recden codes are in the following Code Edit Region for convenience
#Find the leading monomial of Rpoly format
#X1 is the list of variables which we do want to find lm w.r.t them i.e lm(F,[x,y],[x,y,z],p)
#Note:X1 and X are not necessarily equal
lmrpoly:=proc(F,X1,X)
local f,lc,R,M,E,Rnew,m,lm1 ;
#switch to the simple polynomial to find its leading monomial
#rpoly command is back and forth
f:=rpoly(F);
lc:=lcoeff(f,X1,'lm');
lm1:=rpoly(lm,X);
R:=getring(F);
#reconstruct R
M:=op(1,R);
E:=op(3,R);
Rnew:=[M,X,E];
m:=POLYNOMIAL(Rnew,op(2,lm1));
return(m);
end:
#Find the smallest monomial in f
tmrpoly:=proc(F,X1,X)
local f,lc,lm1,m,R,E,Rnew,M;
f:=rpoly2maple(op(2,F),X);
lc:=tcoeff(f,X1,'lm');
lm1:=rpoly(lm,X);
R:=op(1,F);
M:=op(1,R);
E:=op(3,R);
Rnew:=[M,X,E];
m:=POLYNOMIAL(Rnew,op(2,lm1));
return(m);
end:
#Finding the denominator of f
#If f belongs to Q[x], then the denominator of f is the smallest positive integer such that den(f)*f belongs to Z[x]
#input: a rpoly F
#output: den(F)
den:=proc(F)
idenomrpoly(F);
end:
#Semi-Associate of f
#Semi-associate of f is defined as rf s.t r is the smallest rational for which den(f.r)=1
#Example: if f(x,y,z)=2/5*x^2+4/3*y^3+8/11*z^2 then den(f)=5*11*3 but r=den(f)/gcd(2,4,8)
semi:=proc(F)
mulrpoly(1/icontrpoly(F),F);
end:#Monomlist, RED, Basemaker,LAMinpoly
#LAminpoly
#Monomlist
#input: a polynomial f, a list of monomials called B, the list of variables called X
#Output:Coordinate vector of f w.r.t B
with(LinearAlgebra):
Monomlist:=proc(f,B1,X)
local m,v,i,mon,j,cof,c,d,B;
#At some examples B1 is Vector and the nops(B1) does not give the correct result so we need to convert it to list
B:=convert(B1,list);
d:=nops(B);
v:=Vector(d);
c:=coeffs(f,X,'m');
cof:=[c];
mon:=[m];
for j to nops(mon) do
if member(mon[j],B,'i') then
v[i]:=cof[j];
else error "monomial not in B"
fi;
od;
return v;
end:
#RED
#input: a polynomial, a, the list of minimal polynomials,m, and the list of variables,X.
#output: the polynomial ,a, mod m[1],m[2],...
RED:=proc(a,m,X)
local a1,i;
a1:=a;
for i to nops(m) do
a1:=rem(a1,m[-i],X[-i]);
od;
return(a1);
end:
#Basemaker
#input:The list of variables, the list of degree of each variable
#output: The list of all the possible monomials constructed by X mod m[i]s
Basemaker := proc(X,D::list(integer))
local x,i,M,m,Z;
if nops(X)=1 then
x := X[1];
seq( x^i,i=0..D[1]-1 );
else
M := Basemaker(X[2..-1],D[2..-1]);
x := X[1];
Z := [seq( x^i,i=0..D[1]-1 )];
[seq(seq( x*m, m in M ), x in Z)]
fi
end:
LAMinpoly:=proc(m,p,gamma,X,Z)
local n,deg,d,monombasis,g,A,B,det,q,i,j,b,M,st;
n:=nops(m); #nops(m)=nops(X)
deg:=[seq(degree(m[i],X[i]),i=1..n)];
d:=mul(deg[i],i=1..n);
monombasis:=Basemaker(X,deg,m);
#It is a basis for K=Z_p[z1,z2]/<z1^2-2,z2^2-3>: {1,z1,z2,z1z2}
g[0]:=rpoly(1,getring(gamma));
for i to d do
g[i]:=mulrpoly(gamma,g[i-1]);
od;
A:=Matrix(d);
#Creating the matrix of coeffs
for i from 1 to d do
g[i-1] := expand(rpoly(g[i-1]));
A[1..d,i]:=Monomlist(g[i-1],monombasis,X);
od;
g[d] := expand(rpoly(g[d]));
b:=Monomlist(g[d],monombasis,X);
# Construct the augmented matrix B = [A|I|b] and row reduce it
B:=Matrix(d,2*d+1,datatype=integer[8]);
B[1..d,1..d] := A;
for i to d do B[i,d+i] := 1; od;
for i to d do B[i,2*d+1] := b[i] od;
#Check for appropriate C here: If det(A)<>0, then A is invertible and we are good.
#Otherwise, MGCD goes back and chooses another C
#if Det(A) mod p=0 then return(FAIL); fi;
Modular:-RowReduce(p,B,d,2*d+1,2*d+1,'det',0,0,0,0,true);
if det=0 then return FAIL fi;
q := -B[1..d,2*d+1]; # Solve of Aq=-b for q
B := B[1..d,d+1..2*d]; # A^(-1)
#B:=Inverse(A) mod p;
#q:=-B.b mod p; #As det(A)<>0 we are sure that there is a solution
#Constructing the polynomial
M:=add(q[i]*Z^(i-1),i=1...d)+Z^d mod p;
return M,A,B;
end:
#phiminpoly, Phi isomorphism#The isomorphism phi
#Let m=[m1(z1),m2(z2),..,mn(zn)] be the list of minimal polynomials in maple version.
#To build the ring of Q[z1,..zn]/<m1..mn> in RECDEN format,we need to present m2 as a polynomial
#over Q[z1]/<m1>, m3 as a polynomial over Q[z1,z2]/<m1,m2>,..,and m_n as a polynomial
#over Q[z1,z2,..z_(n-1)]/<m1,m2,..,m_(n-1)>
with(ListTools):
phiminpoly:=proc(m,Xmin1)
local n,M,i,Xmin;
n:=nops(m);
M:=[seq(1..n)];
Xmin:=Reverse(Xmin1);
M[1]:=rpoly(m[1],Xmin[n]);
for i from 2 to n do
Xmin[-i..n];
M[i]:=phirpoly(rpoly(m[i],Xmin[-i..n]),op(M[1..i-1]));
getring(M[i]);
od;
return(M);
end:
with(ListTools):
with(PolynomialTools):
with(ArrayTools):
with(RandomTools):
#phi : L=Q(a1,..an)---->K=Q(gamma) is an isomorhism from Q[x1,..xn]/<m1,..mn> to Q[z]/<m(z)>. If we want to call phi, then flag:=1 else if we want to call phi^(-1) then flag:=-1
#X is the list of minimal polynomials' variables
#If we call phi then the minpoly of the single extension is w.r.t Z:=op(indets(M)
#If we call phi^(-1) then we need to have the variable of the single extension
#Inputs: the polynomial f1, m=[m1,..mn]=list of minimal polynomials from Q[x1,..xn]/<m1,..mn>,Xmin=minimal polynomials' variables, Xpoly=polynomial's variables, flag=1 or -1, gama,M,A,Ainv= The out puts of LAMinpoly
Phi:=proc(f1,m,Xmin,Xpoly,flag,M,A,Ainv)
local n,Z,deg,d,monombasis,F,final,newf,i,basis2,f,fmin,Ringfmin,Ring,R1,R2,R3,R33,Rfinal,R, phimin,XT,Xmin1,Rf;
#The variable of the single extension's minimal polynomial
Z:=op(indets(M));
Ring:=getring(f1);
#The characteristic of the field
R1:=op(1,Ring);
n:=nops(m); #nops(m)=nops(X)
deg:=[seq(degree(m[i],Xmin[i]),i=1..n)];
d:=degree(M); #d:=mul(deg[i],i=1..n);
monombasis:= Basemaker(Xmin,deg); #set it as an input
f:=expand(rpoly(f1));
if flag=1 then #we want to calculate Phi(f) so f is in L and we need Ainv
R2:=[op(Xpoly),Z];
#Get the RECDEN format of the minimal polynomial M
fmin:=phirpoly(rpoly(1,Z),rpoly(M,[Z])); #Make the input M as RECDEN
Ringfmin:=getring(fmin);
#This is the Minpoly of single extension in recden version [[],[],[]]
R33:=[op(1,op(3,Ringfmin))];
#The ring of single extension is Rfinal
Rfinal:=[R1,R2,R33];
F:= Monomlist(f,monombasis,Xmin);
final:=Ainv.F;
#write final as a polyomial w.r.t Z in Q[Z]/<m(Z)>
newf:=add(final[i]*Z^(i-1),i=1...d);
return rpoly(newf,Rfinal); #It is costy to convert
elif flag=-1 then #We want to call phi^(-1)
#Note: Every two bases of a vector space have the same cardinality. Hence, if the base of Q[x,y]/<m1,m2> has d items then the basis of Q[z]/<m(z)> has d items too. Remember these two vector spaces are isomorphic.
#When we call phi^(-1) we have a polynomial over Q[Z]/<m(Z)> and we want to map it to the corresponding polynomial in Q[x,y]/<m1,m2>
# write f w.r.t the basis ord:=[1,Z,Z^2,..,Z^(d-1)] of Q[Z]/<m(Z)> where deg(M(Z))=d
basis2:=Array(1..d);
for i to d do
basis2[i]:=Z^(i-1); #A basis for the single extension Q[Z]/<M(Z)>
od;
F:=Monomlist(f,basis2,Z); #The coordinate of f relative to basis2
final:=A.F;
newf:=expand(RED((add(final[i]*monombasis[i],i=1..d)),m,Xmin)); #it's in Q[x,y]/<m1,m2>
#Now, we are to construct the ring of multiple extensions
phimin:=phiminpoly(m,Xmin);
Xmin1:= [ seq( Xmin[-i], i=1..nops(Xmin) ) ]; #Reverse did not work here:(
XT:=[op(Xpoly),op(Xmin1)];
#R1; #The characteristic of the field
if R1<>0 then #If the characteristic of the field <>0
newf:=phirpoly(rpoly(newf,XT),op(phimin),R1);
else #If the characteristic of the field =0
newf:=phirpoly(rpoly(newf,XT),op(phimin));
fi;
return(newf);
fi;
end:#PGCD
#X is the list of the polynomials' variables, XT is the list of the whole variables
#We did not use minpoly in PGCD So PGCD is the same for single & multiple extensions
PGCD:= proc(A,B,p)
local k,Xk,minmon,least,conta,contb,c,ap,bp,g,interpol,lcinterpol,j,Aj,Bj,Cj,gj,G,lcCj,lcA,lcB,lcAk, lcBk,Randfunction,LM,lm,ring,ring2,count,XT,X,n,m,AL,BL,L,i,GCDL,deginterpol,prod, interpolj,vj,prodj,q;
#Check if A and B are in the same field
checkconsistency(A,B);
ring:=getring(A);
#XT is the whole variables containing the polynomial's variables and the minpolys variables
XT:=getvars(A);
n:=nops(XT);
m:=nops(getalgexts(A));
X:=XT[1..n-m];
#k must be the number of polynomial variables
k:=nops(X);
Xk:=X[1..k-1]; #Polynomial variables= The whole variables - the variable of the minpoly
#Base case
if k = 1 then
G:=gcdrpoly(A,B);
#G:=scarpoly(lcrpoly(G),G); #Monic GCD
return G;
fi;
#Calculate the Content/Primpart of A and B
conta:=contrpoly(A,Xk,'ap');
contb:=contrpoly(B,Xk,'bp');
#c is the gcd of the content of A and B
c:=gcdrpoly(conta,contb);
#g is the gcd of the leading coefficients of pp of A and B w.r.t Xk
g:=gcdrpoly(lcrpoly(ap,Xk),lcrpoly(bp,Xk));
lcA:=lcrpoly(ap,X);
lcB:=lcrpoly(bp,X);
#q is the array of eval points. We need tracking q to avoid repetitious evaluation points
q:=[];
#Randfunction is the function producing the rand eval points
Randfunction:=rand(p);
lcAk:=lcrpoly(ap,Xk);
lcBk:=lcrpoly(bp,Xk);
#count is the counter of the acceptable eval points
count:=0;
while true do
do
j:= Randfunction() mod p;
until not member(j,q) and evalrpoly(lcAk,X[k]=j) mod p<>0 and evalrpoly(lcBk,X[k]=j) mod p<>0;
Aj:=evalrpoly(ap,X[k]=j);
Bj:=evalrpoly(bp,X[k]=j);
Cj:=PGCD(Aj,Bj,p);
#lm=The leading monomial of Cj
lm:=lmrpoly(Cj,Xk,XT,p);
#m=min(lm,n) where n=min(lm(ap),lm(bp))
gj:=evalrpoly(g,X[k]=j);
Cj:=scarpoly(gj,Cj);
#To count the number of acceptable evaluation points
if count=0 then
minmon:=lm;
least:=lm;
count:=count+1;
deginterpol:=0;
interpol:=Cj;
prod:=rpoly(X[k]-j,ring);
q := [j]; # missing !!!
else
least:=tmrpoly(addrpoly(lm,minmon),X,XT,p);
if lm = least and minmon = least then
prodj:=evalrpoly(prod,X[k]=j);
#print(PRODJ(prod,prodj),q,j);
q:=[op(q),j]; #add j to q
if count=1 then
interpolj:=interpol; #When count=1, interpol=Cj and Cj does not have X[k] in it so evalrpoly returns an error.
else
interpolj:=evalrpoly(interpol,X[k]=j);
fi;
vj:=mulrpoly(invrpoly(prodj),subrpoly(Cj,interpolj));
#The ring to which vj, interpol, and prod belong must be the same but at the second iteration when count=1, vj and interpol do not have y as their variable so we must reconstruct the ringto which vj and interpol belong.
vj:=rpoly(rpoly(vj),ring); # adding X[k]
interpol:=rpoly(rpoly(interpol),ring);
interpol:=addrpoly(interpol,mulrpoly(prod,vj));
prod:=mulrpoly(prod, rpoly(X[k]-j,ring));
deginterpol:=degrpoly(interpol,X[k]); #We did not initialize interpol since its ring changes at each iteration and it causes error. If we fixed a ring for interpol then it was more expensive than simply initialize its degree of y =0 at the first iteration
count:=count+1;
elif lm = least and minmon<>least then
printf(" All the evaluation points before j=%d are unlucky",j);
q:=[j];
prod:=rpoly(X[k]-j,ring); # this was missing!!
interpol:=Cj;
minmon:= least;
count:=1;
elif lm<>least then #Recognizing as an unlucky prime!
printf("%d is an unlucky evaluation point\134n",j);
fi;
fi;
if count>deginterpol+1 then #Line 90 is equivalent to line 30 of PGCD algorithm in the paper
G:=pprpoly(interpol,Xk);
do
#Make our algorithm a Monte Carlo algorithm
L:=[seq(X[i]= Randfunction() mod p, i=2..k)];
AL, BL, GCDL:= evalrpoly(A,L), evalrpoly(B,L), evalrpoly(G,L);
until degrpoly(GCDL,X[1])=degrpoly(G,X[1]);
if divrpoly(AL,GCDL) and divrpoly(BL,GCDL) then
G:=mulrpoly(c,G);
return (G);
fi;
fi;
od;
end:
untrace(PGCD);#MGCD
with(ListTools):
#Takes A,B in Q[z]/<m(z)>[x1,..xn] and gives GCD(A,B)
MGCD:= proc(A1,B1,sprime)
local A,B,ringA,X,Xpoly,Minpolys,cc,i,m,n,Xmin,conta,contb,k,ap,bp,LCap,LCbp,lmap,lmbp,minmon,H,prod,p,Ap,Bp,Cp,lcCp,lmCp,least,CRT,H2,test1,test2,summon,irat,M,gama,gamma,phimatrix,invphimatrix,PhiAp,PhiBp, lcmintest,xx,yy,ca,cb,nca,ncb,C,counter,j,s,nfac,fac,ss,l,LA,Rand,ZD,L,AL,BL,GCDL,wt,LCH2,slc,randfunction,minpolys,st,tt,laminpolytime,pgcdtime;
#We use semi of the inputs to eliminate fractions
tt := time();
laminpolytime := 0;
pgcdtime := 0;
A:=semi(A1);
B:=semi(B1);
checkconsistency(A,B,alpha);
ringA:=getring(A);
#get the whole variables
X:=getvars(A);
n:=nops(X);
#get the minimal polynomials
minpolys:=getalgexts(A); #The list of minimal polynomials
m:=nops(minpolys);
#get the polynomials' variables
Xpoly:=X[1..n-m];
#get the minimal polynomials' variables
Xmin:=Reverse(X[n-m+1..n]);
#We need to convert minpolys from rpoly to simple polynomials
Minpolys := map(rpoly,minpolys);
ap, bp:= A,B;
#k is the nops of the variables of polynomials
k:= nops(Xpoly);
H:=[];
#we do not need bound when we are using iratrecon since the loop terminates when iratrecon gives the proper answer
if nargs=2 then p := 2^31; else p:= sprime-1; fi;
#The main loop
while true do
#recognizing lc-bad and det-bad primes
do
p:= nextprime(p);
printf("MGCD:prime=%d\134n",p);
lcmintest:=true;
for i to m do #p must not divide lc of ^M_i
if den(minpolys[i]) mod p=0 then lcmintest:=false; fi;
od;
#p is a bad prime if p|LC(A)
if rpoly(lcrpoly(A,Xpoly)) mod p = 0 then
printf("p=%d is an lc-bad prime\134n",p);
fi;
until rpoly(lcrpoly(A,Xpoly)) mod p <> 0
and
#As we used the semi of inputs as our inputs, we do not need to check if den(A)& den(B) mod p<>0
lcmintest=true;
#If p | lcoeff of A^,M^, where A^ is the semi associate of A and m^ is the semi-associate of the minimal polynomial then p is a lc-bad prime.
Rand:=rand(1..p-1);
C:=[seq(Rand(),i=1..m-1)];
gama := Xmin[1]+add(C[i]*Xmin[i+1],i=1..m-1);
gamma := phirpoly(rpoly(gama,getcofring(A)),p);
st := time();
LA:=LAMinpoly(Minpolys,p,gamma,Xmin,X[n]);
laminpolytime := laminpolytime + time()-st;
if LA=FAIL then
printf("p=%d is a det-bad prime or C=%a is an inappropriate list of constants\134n",p,C);
next;
fi;
M,phimatrix,invphimatrix:=LA;
#printf("gamma:=%a and M(z)=%a\134n",gama,M);
Ap:= phirpoly(ap,p); #Now everything is in Zp
Bp:= phirpoly(bp,p);
PhiAp:=Phi(Ap,Minpolys,Xmin,Xpoly,1,M,phimatrix,invphimatrix); #In Z_p(gamma)
PhiBp:=Phi(Bp,Minpolys,Xmin,Xpoly,1,M,phimatrix,invphimatrix);
st := time();
Cp:=traperror(PGCD(PhiAp,PhiBp,p));
pgcdtime := pgcdtime+time()-st;
#if there is a zero divisor then go to the first loop and change the prime
#print Z-D primes
if Cp=lasterror and nops([Cp])=2 and Cp[1]="inverse does not exist" then
ZD:=Cp[2]; printf("p=%d is a ZD prime ZD=%a\134n",p,ZD);
next;
elif Cp=lasterror then error lasterror;
fi;
#When the modular gcd is constant we are done!
if degrpoly(Cp)=0 then printf("Cp=%a is a constant",Cp); return(rpoly(1,X)) fi;
Cp:=Phi(Cp,Minpolys,Xmin,Xpoly,-1,M,phimatrix,invphimatrix);
getring(Cp);
lcCp:= lcrpoly(Cp,Xpoly[1..k-1]);
lmCp:=rpoly(lmrpoly(Cp,Xpoly[1..k-1],X));
if not assigned(minmon) then
CRT:=Cp;
minmon:=lmCp;
summon:=minmon;
least:=minmon;
H:=[Cp];
else
summon:=lmCp+minmon;
least:=rpoly(tmrpoly(rpoly(summon,X),Xpoly,X));
Cp:= scarpoly(invrpoly(lcrpoly(Cp)),Cp);
#Test for unlucky primes
if lmCp = least and minmon<>least then
if p>sprime then printf("p=%d and All the previous primes were unlucky\134n",p); fi;
H:=[Cp];
CRT:=Cp;
minmon:=least;
elif lmCp = least and minmon= least then
H:=[op(H),Cp];
CRT:=ichremrpoly(map(retextsrpoly,H));
elif lmCp<>least then
printf("%a is an unlucky prime\134n",p);
fi;
fi;
irat:=irrrpoly(CRT);
if irat<> FAIL then
H2:=subsop(1=ringA,irat);#We do not use the bound.When RNR has a correct solution the loop termintes.
#Make our algorithm a Monte Carlo algorithm
randfunction:=rand(p);
do
#Make our algorithm a Monte Carlo algorithm
L:=[seq(Xpoly[i]= randfunction(), i=2..k)];
AL, BL, GCDL:= evalrpoly(A,L), evalrpoly(B,L), evalrpoly(H2,L);
until degrpoly(GCDL,Xpoly[1])=degrpoly(H2,Xpoly[1]);
if divrpoly(AL,GCDL) and divrpoly(BL,GCDL) then
tt := time()-tt;
if timings=true then printf("MGCD: time=%.3f LAminpoly=%.3f PGCD=%.3f\134n",tt,laminpolytime,pgcdtime); fi;
return (H2);
fi;
fi;
od;
end:#This MGCD does not reduce the multiple extension to a single extension#MGCD using multiple extensions
with(ListTools):
#Takes A,B in Q[z]/<m(z)>[x1,..xn] and gives GCD(A,B)
MGCD2:= proc(A1,B1,sprime)
local A,B,ringA,X,Xpoly,Minpolys,minpolys,cc,i,m,n,Xmin,conta,contb,k,ap,bp,LCap,LCbp,lmap,lmbp,minmon,H,prod,p,Ap,Bp,Cp,lcCp,lmCp,least,CRT,H2,test1,test2,summon,irat,M,gama,phimatrix,invphimatrix,PhiAp,PhiBp, lcmintest,xx,yy,ca,cb,nca,ncb,C,counter,j,s,nfac,fac,ss,l,LA,Rand,ZD,L,AL,BL,GCDL,wt,LCH2,slc,randfunction,st,tt,pt;
#We use semi of the inputs to eliminate fractions
tt := time();
pt := 0.0;
A:=semi(A1);
B:=semi(B1);
checkconsistency(A,B,alpha);
ringA:=getring(A);
#get the whole variables
X:=getvars(A);
n:=nops(X);
#get the minimal polynomials
minpolys:=getalgexts(A); #The list of minimal polynomials
m:=nops(minpolys);
#get the polynomials' variables
Xpoly:=X[1..n-m];
#get the minimal polynomials' variables
Xmin:=Reverse(X[n-m+1..n]);
#We need to convert minpolys from rpoly to simple polynomials
Minpolys := map(rpoly,minpolys);
ap, bp:= A,B;
#k is the nops of the variables of polynomials
k:= nops(Xpoly);
H:=[];
#we do not need bound when we are using iratrecon since the loop terminates when iratrecon gives the proper answer
if nargs=2 then p := 2^31 else p:= sprime-1 fi;
#The main loop
while true do
#recognizing lc-bad and det-bad primes
do
p:= nextprime(p);
printf("MGCD2:prime=%d\134n",p);
lcmintest:=true;
for i to m do #p must not divide lc of ^M_i
if den(minpolys[i]) mod p=0 then
lcmintest:=false;
fi;
od;
#p is a bad prime if p|LC(A)
if rpoly(lcrpoly(A,Xpoly)) mod p = 0 then
printf("p=%d is an lc-bad prime\134n",p);
fi;
until rpoly(lcrpoly(A,Xpoly)) mod p <> 0
and
#As we used the semi of inputs as our inputs, we do not need to check if den(A)& den(B) mod p<>0
lcmintest=true; #If p | lcoeff of A^,M^, where A^ is the semi associate of A and m^ is the semi-associate of the minimal polynomial then p is a lc-bad prime.
Ap:= phirpoly(ap,p); #Now everything is in Zp
Bp:= phirpoly(bp,p);
st := time();
Cp:=traperror(PGCD(Ap,Bp,p));
pt := pt+time()-st;
#if there is a zero divisor then go to the first loop and change the prime
#print Z-D primes
if Cp=lasterror and nops([Cp])=2 and Cp[1]="inverse does not exist" then
ZD:=Cp[2]; printf("p=%d is a ZD prime ZD=%a\134n",p,ZD);
next;
fi;
#When the modular gcd is constant we are done!
if degrpoly(Cp)=0 then printf("Cp=%a is a constant",Cp); return(rpoly(1,X)) fi;
lcCp:= lcrpoly(Cp,Xpoly[1..k-1]);
lmCp:=rpoly(lmrpoly(Cp,Xpoly[1..k-1],X));
if not assigned(minmon) then
CRT:=Cp;
minmon:=lmCp;
summon:=minmon;
least:=minmon;
H:=[Cp];
else
summon:=lmCp+minmon;
least:=rpoly(tmrpoly(rpoly(summon,X),Xpoly,X));
Cp:= scarpoly(invrpoly(lcrpoly(Cp)),Cp);
#Test for unlucky primes
if lmCp = least and minmon<>least then
if p>sprime then printf("p=%d and All the previous primes were unlucky\134n",p); fi;
H:=[Cp];
CRT:=Cp;
minmon:=least;
elif lmCp = least and minmon= least then
H:=[op(H),Cp];
CRT:=ichremrpoly(map(retextsrpoly,H));
elif lmCp<>least then
printf("%a is an unlucky prime\134n",p);
fi;
fi;
irat:=irrrpoly(CRT);
if irat<> FAIL then
H2:=subsop(1=ringA,irat);#We do not use the bound.When RNR has a correct solution the loop termintes.
#Make our algorithm a Monte Carlo algorithm
randfunction:=rand(p);
do
#Make our algorithm a Monte Carlo algorithm
L:=[seq(Xpoly[i]= randfunction(), i=2..k)];
AL, BL, GCDL:= evalrpoly(A,L), evalrpoly(B,L), evalrpoly(H2,L);
until degrpoly(GCDL,Xpoly[1])=degrpoly(H2,Xpoly[1]);
if divrpoly(AL,GCDL) and divrpoly(BL,GCDL) then
tt := time()-tt;
if timings=true then printf("MGCD2: time=%.3f PGCD=%.3f\134n",tt,pt); fi;
return (H2);
fi;
fi;
od;
end:
untrace(PGCD):
mini1:=z^2-2;
mini2:=w^2-3;
M1:=rpoly(mini1,z);
M2:=phirpoly(rpoly(mini2,[w,z]),M1);
f1:=phirpoly(rpoly((w+5)*(x+w+y)*(14*x+2*w+z*s),[x,y,s,w,z]),M1,M2);
f2:=phirpoly(rpoly((x+w+y)*(x+2*w+z+s),[x,y,s,w,z]),M1,M2);
timings := false:
Mahsa1:=MGCD(f1,f2,5);MGCD2(f1,f2,5);M1:=z1^2-2;
M2:=z2^2-3;
M3:=z3^2-5;
M4:=z4^2-7;
M5:=z5^2-11;
g:=x^2+z1*z2*x*y+53124/15321*z3*z4*x-z5*y^2-z1*z3*y+z2*z4+12345/67891:
g:=rpoly(g,[x,y,z5,z4,z3,z2,z1],[M1,M2,M3,M4,M5]);
a:=rpoly(x^2+z1*x+z2*z3*y+z4*z5*y^2+1,[x,y,z5,z4,z3,z2,z1],[M1,M2,M3,M4,M5]);
b:=rpoly(x^2+z5*x+z1*z3*y+z4*z3*y^2-1,[x,y,z5,z4,z3,z2,z1],[M1,M2,M3,M4,M5]);k:=1;
timings := true;
T:=[seq(0,k=1..10)];
for k to 5 do
f1:=mulrpoly(g,powrpoly(a,k));
f2:=mulrpoly(g,powrpoly(b,k));
printf("k=%d deg(f1,x)=%d deg(f1,y)=%d deg(f2,x)=%d deg(f2,y)=%d deg(g,x)=%d deg(g,y)=%d\134n",
k,degrpoly(f1,x),degrpoly(f1,y),degrpoly(f2,x),degrpoly(f2,y),degrpoly(g,x),degrpoly(g,y));
t:=time():
Mahsa1:=MGCD(f1,f2,2^31);
tMGCD1:=time()-t;
tt:=time():
Mahsa2:=MGCD2(f1,f2,2^31);
tMGCD2:=time()-tt;
printf("Check=%a\134n",evalb(Mahsa1=Mahsa2));
T[k]:=tMGCD2/tMGCD1;
od:
T;A := powrpoly(a,5):
B := powrpoly(b,5):
for k from 6 to 10 do
A:=mulrpoly(a,A);
B:=mulrpoly(b,B):
f1:=mulrpoly(g,A);
f2:=mulrpoly(g,B);
printf("k=%d deg(f1,x)=%d deg(f1,y)=%d deg(f2,x)=%d deg(f2,y)=%d deg(g,x)=%d deg(g,y)=%d\134n",
k,degrpoly(f1,x),degrpoly(f1,y),degrpoly(f2,x),degrpoly(f2,y),degrpoly(g,x),degrpoly(g,y));: MGCD(f1,f2,2^31);
od: