#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:
#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:

#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:
#ures computes the resultant of two univariate polynomials using MEA
ures:=proc(a, b)
local q,r,s,t,k,n,m,h,R,c,re,Ring,v,prod,i,l;
r[0]:=a;
r[1]:=b;
m[0]:=r[0];
n[0]:=degrpoly(a);
n[1]:=degrpoly(b);
k:=1;
re:=[a];
v:=0;
Ring:=getring(a);
prod:=rpoly(1,Ring);
while r[k]<>rpoly(0,Ring) do
   m[k]:=monic(r[k]);
   re:=[op(re),r[k]];
   r[k+1]:=remrpoly(m[k-1],m[k]);
   n[k+1]:=degrpoly(r[k+1]);
   #If gcd(f,g)<>constant then res(f,g)=0
   if r[k+1]=rpoly(0,Ring) and n[k]<>0  then return(0); fi;
   prod:=scarpoly(powrpoly(lcrpoly(r[k]),n[k-1]),prod);
   v:=v+n[k]*n[k-1];
   k:=k+1;
od;
#v:=add(n[i]*n[i+1],i=1..k-2);
prod:=scarpoly((-1)^v,prod);
return(prod);

end:
removexrpoly:=proc(f,x) 
#remove the main variable from the ring of f and return f with the new ring
  local Rf,X,n,Y,i,R,nr,j,f1,ff,XX,M,E;
   M,XX,E:=op(getring(f));
   if Search(x,XX)=0 then return(f); fi;
 
  f1:=makemainrpoly(f,x);
  Rf:=getring(f1);
  X:=op(2,Rf);
  n:=nops(X);
  Y:=[seq(0,i=1..n-1)]; #we use this in PRS where there is just one extention and
  for i to n-1 do
     Y[i]:=X[i+1];
  od;
  nr:=nops(Rf);
  R:=[seq(0,i=1..nr)];
  R[1]:=op(1,Rf);
  R[2]:=Y;
  for j from 3 to nr do 
   R[j]:=op(j,Rf);
  od;
  ff:=rpoly(f);
  return( rpoly(ff,R));
end:
#PRS
#X is the list of the polynomials' variables, XT is the list of the whole variables
PRS:= proc(AA,BB,x,p)
local A,B,H,k,Xk,conta,contb,c,ap,bp,g,interpol,j,Aj,Bj,Cj,gj,R,lcAk, lcBk,Randfunction,ring,count,XT,X,n,m,i,q,vj,interpolj,prodj,prod,dya,dyb,Bound1,Bound2;
#To compute the resultant w.r.t x_i,we must make x_i the main variable
   A:=makemainrpoly(AA,x);
   B:=makemainrpoly(BB,x);
#Check if A and B are in the same field
   checkconsistency(A,B);
   ring:=getring(A);
#XT is the list of all the 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. In MRS we call RSP after phi_{gamma} which means we just have one extension
   k:=nops(X);
   Xk:=X[1..k-1];
#Xk := X[2..k];  #Polynomial variables= The whole variables - the variable of the minpoly
   #Base case
   if k = 1 then  
      R:=ures(A,B,x); #we use MEA to compute the resultant of two univariate polynomials
      return(R);    
  fi;
#Calculate the Content/Primpart of  A and B w.r.t Xk to use it for identifying lc-bad evaluation points
   conta:=contrpoly(A,Xk,'ap');
   contb:=contrpoly(B,Xk,'bp');
   lcAk:=lcrpoly(ap,Xk);  
   lcBk:=lcrpoly(bp,Xk);
#q is the array of eval points. We need tracking q to avoid repetitious evaluation points and H is the array of interpolated polynomials
   H:=[];
   q:=[];
   #Randfunction is the function producing the rand eval points 
   Randfunction:=rand(p);
   #count is the counter of the acceptable eval points
   count:=0;   #At the end of the loop, we need to compare the new resulatant with the previous resultant which we initialize it to be A
   dya,dyb:=degrpoly(A,x),degrpoly(B,x);
  Bound1:=dyb*degrpoly(A,X[k])+dya*degrpoly(B,X[k]);
  #Bound2:=degree(rpoly(A),{x,X[k]})+degree(rpoly(B),{x,X[k]});
  for i from 0 to Bound1 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(A,X[k]=j);
      Bj:=evalrpoly(B,X[k]=j);
      Cj:=PRS(Aj,Bj,x,p);
      if count=0 then 
         q:= [j]; 
         H:=[Cj]; 
         interpol:=Cj;
         prod:=rpoly(X[k]-j,ring);
         count:=count+1;
      else 
         q:=[op(q),j];  #add j to q
         H:=[op(H),Cj]; #add Rj to H
         prodj:=evalrpoly(prod,X[k]=j);
         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(removexrpoly(invrpoly(prodj),x),subrpoly(removexrpoly(Cj,x),removexrpoly(interpolj,x)));
          getring(vj);
getring(prod);
             #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)); #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;
      fi;   
   od;
     R:=interpol;
     return (R);
end:
#MRS with primitive element
with(ListTools): 
#Input: A,B in Q[z1,..zn]/<m1,..mn>[x1,..xk] and variable x_i
#Output: resultant(A,B,x)  
MRS:= proc(A1,B1,x,sprime)
local count, previousresult,A,B,ringA,X,Xpoly,Minpolys,i,m,n,Xmin,k,H,prod,p,Ap,Bp,Cp,CRT,H2,test1,test2,irat,M,gama,gamma,phimatrix,invphimatrix,PhiAp,PhiBp, lcmintest,ca,cb,nca,ncb,C,counter,j,s,LA,Rand,ZD,wt,randfunction,minpolys,st,tt,laminpolytime,pgcdtime,PPresult,Nump;
#We use semi of the inputs to eliminate fractions
   previousresult:=0; #first result test
   PPresult:=0; #second result test
   count:=0;
   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);
   Minpolys :=map(rpoly,map(semi,minpolys));
   #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 returns a repeatious answer
   if nargs=3 then p := 2^31; else p:= sprime-1; fi;
   #The main loop
   Nump:=0;
   while true do
   #recognizing lc-bad and det-bad primes
      do
         p:= nextprime(p);
        # printf("MRS: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;  printf("p=%d is an lc-bad prime minpoly case\n",p);  fi;
         od;
         #p is a bad prime if p|LC(A)  or p|LC(A)
         if rpoly(lcrpoly(A,Xpoly)) mod p = 0 or  rpoly(lcrpoly(B,Xpoly)) mod p = 0 then 
            printf("p=%d is an lc-bad prime polynomial case\n",p); 
         fi;
      until rpoly(lcrpoly(A,Xpoly)) mod p <> 0 and rpoly(lcrpoly(B,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;
   
   Ap:= phirpoly(A,p); #Now everything is in Zp[z]/<M(z)>
   Bp:= phirpoly(B,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(PRS(PhiAp,PhiBp,x,p));
   #if verify(Cp,resrpoly(PhiAp,PhiBp,x)) then printf("PRS is working well\n"); fi;
   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;
   #phi_gamma inverse
   Cp:=Phi(Cp,Minpolys,Xmin,Xpoly,-1,M,phimatrix,invphimatrix);
   if count=0 then 
      CRT:=Cp;
      H:=[Cp];
      count:=count+1;
   else
      H:=[op(H),Cp];
      CRT:=ichremrpoly(map(retextsrpoly,H));
      count:=count+1;     
   fi;
   Nump:=Nump+1; #Number of good primes
   irat:=irrrpoly(CRT);
   if  irat<> FAIL then 
#We do not use the bound.When the result of RNR does not change the loop termintes.This strategy makes our algorithm  Monte Carlo 
     if irat=previousresult and irat=PPresult then
         tt := time()-tt;
         if timings=true then printf("MRS: time=%.3f  LAminpoly=%.3f  PGCD=%.3f Nump=%d\n",tt,laminpolytime,pgcdtime, Nump); fi;
          irat:=subsop(1=ringA,irat);
         return (removexrpoly(irat,x),count); 
      else  PPresult:=previousresult; previousresult:=irat;  
      fi;
    fi;

od;
end: