# This version was used to get the timings for the benchmark

restart;
with(PolynomialTools):
currentdir("/Users/mmonagan/Desktop");
read "recden";
#The recden codes are in the following Code Edit Region for convenience

GCD := Algebraic:-RecursiveDensePolynomials:-gcdrpoly;
MUL := Algebraic:-RecursiveDensePolynomials:-mulrpoly;

# This uses implementation of gcdrpoly and mulrpoly in the Maple kernel to try to
# get decent timings.  Michael Monagan

#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]:=MUL(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:=GCD(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:=GCD(conta,contb);
   #g is the gcd of the leading coefficients of pp of A and B w.r.t Xk
   g:=GCD(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\n",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\n",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\n",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\n",p,C); 
            next;
         fi;

 M,phimatrix,invphimatrix:=LA;
#printf("gamma:=%a and M(z)=%a\n",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\n",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\n",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\n",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\n",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\n",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\n",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\n",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\n",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\n",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\n",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..40)];
A := a;
B := b;
for k to 40 do
  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\n",
     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\n",evalb(Mahsa1=Mahsa2));
  T[k]:=tMGCD2/tMGCD1;
  A := mulrpoly(a,A);
  B := mulrpoly(b,B);
od:
T;