with(ArrayTools);

CCodeDir := "./";

EVAL64 := define_external('eval64s',
    FF::integer[8],
    XX::integer[8],
    AA::integer[8],
    pp::integer[8], LIB=cat(CCodeDir,"evalGP1.so"),
    RETURN::integer[8]):

pointer := proc(f) option inline; addressof(f)-4*(2^63-1); end:

Eval64 := proc(f,X,alpha,p)
    if type(f,integer) then modp(f,p); else EVAL64(pointer(f),pointer(X),pointer(alpha),p); fi;
end:

#Interpolate a univariate polynomial from black box A (pre computed)
#X - evaluation points
#x - Interpolate in variable x
#j - index for particular variable
#p - prime
#degBound - if you know the degree bound for a particular polynomial
ProbeInterpUni := proc(A,X,xj,j,p,degBound := 3)

   local i, m, ra, Y, g,g2,f,xPoints,yPoints, check;

   ra := rand(1..p-1);
   m := nops(X);

   check := 0;
   g := -1;
   Y := X;

   Y[j] := ra();
   xPoints := [Y[j]]; yPoints := [A(Y,p)];
   g2 := Interp(xPoints,yPoints,xj) mod p;

   while g <> g2 or check < degBound do
      #print("in here");
      g := g2;
      while evalb(Y[j] in xPoints) do
         Y[j] := ra();
      od;
      xPoints := [op(xPoints),Y[j]]; yPoints := [op(yPoints),A(Y,p)];
      g2 := Interp(xPoints,yPoints,xj) mod p;
      if g = g2 then check := check + 1;
      else check := 0; fi;
   od;

   return g;

end proc:
ProbeInterpDegUni := proc(BB,X,x,j,deg,alpha,p)

   local i,xpoints,ypoints,t,ra,f,beta;

   beta := alpha;

   if p-1 < deg+1 then
      printf("p is too small to perform this interpolation");
      return FAIL;
   fi;


   ra := rand(1..p-1);

   xpoints := [seq(0,i=1..deg+1)]; ypoints := [seq(0,i=1..deg+1)];
   #print(xpoints); print(ypoints);

   #fill xpoints with unique elements of Zp
   i := 1;
   while i <= deg+1 do
      t := ra();
      if not evalb(t in xpoints) then
         xpoints[i] := t;
         i := i+1;
      fi;
   od;

#print(xpoints);

   for i from 1 to deg+1 do
      beta[j] := xpoints[i];
      ypoints[i] := BB(beta,p);
   od;

   f := Interp(xpoints,ypoints,x) mod p;

   return f;

end proc:
TdegEstimate := proc(BB,X,p)

   local i,j,k,xx,yy,ra,n,beta,fest1,fest2,count,alpha,gamma;

   n := nops(X);

   ra := rand(1..p-1);
   beta := [seq(ra(),i=1..n)];
   fest1 := 0; fest2 := 1; count := 0;

   xx := []; yy := [];

   while fest1 <> fest2 or count < 3 do
      alpha := ra();
     while evalb(alpha in xx) do alpha := ra(); od;
      xx := [op(xx),alpha];
      yy := [op(yy),BB(alpha*beta,p)];
      fest1 := Interp(xx,yy,x) mod p;
      if fest1 = fest2 then count := count+1;
      else  count := 0; fi;
      fest2 := fest1;
   od;

   return degree(fest1,x);

end proc:


UnivInterpNew := proc(B,X,alpha2,j,p)

   local n,ra,beta,i,M,b,alpha,k,alphak,yk,ba,Ma,vk;

   n := nops(X);

   if j < 1 or j > n then
      return "FAIL";
   fi;

   ra := rand(1..p-1);
   beta := alpha2;
   alpha := [];


   M := 1; b[-1] := 0; k := -1;
   do
     k := k+1;
     alphak := ra();
     if evalb(alphak in alpha) then
        alphak := ra();
     fi;
     alpha := [op(alpha),alphak];
     beta[j] := alphak;
     yk := B(beta,p); #print(yk);
     ba := Eval(b[k-1],x=alphak) mod p;
     Ma := Eval(M,x=alphak) mod p;
     vk := (yk - ba)/Ma mod p;
     b[k] := b[k-1] + vk*M mod p; #print(b[k]);
     M := M*(x-alphak);
     #print(b[k]);

   until vk = 0;


   return b[k];

end proc:
TotalDegNew := proc(BA,X,BETA1,BETA2,p)

   local n,ra,beta,i,j,M,bbb,alpha,k,alphak,yk,ba,Ma,vk,beta2,gamma;

   n := nops(X);

   ra := rand(1..p-1);
   alpha := [];

   M := 1; bbb[-1] := 0; k := -1;
   do
     k := k+1;
     alphak := ra();
     if evalb(alphak in alpha) then
        alphak := ra();
     fi;
     alpha := [op(alpha),alphak]; #print(alpha);
     gamma := [seq(BETA1[j]*alphak + BETA2[j] mod p,j=1..n)];
     yk := BA(gamma,p);
     ba := Eval(bbb[k-1],x=alphak) mod p;
     Ma := Eval(M,x=alphak) mod p;
     vk := (yk - ba)/Ma mod p;
     bbb[k] := bbb[k-1] + vk*M mod p;
     M := M*(x-alphak);
   until vk = 0;


   return bbb[k];

end proc:

uniquePoints := proc(ra,n)

   local i,j,t,flag,Y;

   Y := [seq(0,i=1..n)];
   Y[1] := ra();
   i := 2;
   while i <= n do
      t := ra();
      flag := 0;
      for j from 1 to i-1 do
         if t = Y[j] then flag := 1; fi;
      od;
      if flag = 0 then Y[i] := t; i := i+1; fi;
   od;

   return Y;

end proc:
BMEA := proc(a,p,x) local n,m,R0,R1,V0,V1,i;
  n := iquo( nops(a), 2 );
  m := 2*n-1;
  R0 := x^(2*n);
  R1 := add( a[m+1-i]*x^i, i=0..m );
#lprint("R0=",R0);
#lprint("R1=",R1);
  V0 := 0;
  V1 := 1;
  while n <= degree(R1,x) do
     R0,R1 := R1,Rem(R0,R1,x,'Q') mod p;
#lprint("R0=",R0);
#lprint("R1=",R1);
     V0,V1 := V1,Expand(V0-Q*V1) mod p;
#lprint("V0=",V0);
#lprint("V1=",V1);
  od;
  i := 1/lcoeff(V1,x) mod p;
  i*V1 mod p;
end:
VandermondeSolve := proc( v::{Vector,list}, m::{Vector,list}, p::prime, shift::integer:=0 )
local t,i,j,M,x,a,q,r,s,temp;
   t := numelems(v);
   if numelems(m) <> t then error "v and m must be the same size"; fi;
   printf("Maple code Vandermonde solver: t=%d  p=%d\n",t,p);
   M := 1;
   for r in m do
       M := Expand( M*(x-r) ) mod p;
   od;
   a := Vector(t);
   for j to t do
       q := Quo(M,x-m[j],x) mod p;
       r := 1/Eval(q,x=m[j]) mod p;
       s := 0;
       for i to t do s := s + v[i]*coeff(q,x,i-1); od;
       a[j] := r*s mod p;
       if shift=0 then next fi;
       r := 1/m[j] mod p;
       r := r &^ shift mod p;
       a[j] := r*a[j] mod p;
   od;
   if type(v,list) then a := convert(a,list); fi;
   return a;
end:
Interp1var := proc( B::procedure, X::list, alpha::Array, deg::list, N::integer, j::posint, p::prime )
    # Input: The black box B, X=[x_1,...,x_N], alpha=[alpha_1,...,alpha_N], an evaluation point;
    #        deg is the list of degrees of all variables; N is the no.of variables;
    #        j is the chosen variable index; p is a prime.
    # Output: Det(A(alpha_1,...,x_j,...,alpha_N)) mod p
    local dj, xx, i, F, beta;
    if nops(X)<>N then
        printf("Interp1var: no.of variables not equal to N\n"); return 'Interp1var_FAIL';
    fi;
    if j > N or j < 1 then
        printf("Interp1var: chosen variable index out of range\n"); return 'Interp1var_FAIL';
    fi;
    dj := deg[j]; xx := Array(1..dj+1,i->i);
    beta := Array(alpha); F := Array(1..dj+1);
    for i to dj+1 do beta[j] := i; F[i] := B( beta, p ); od;
    return Interp( xx, F, X[j] ) mod p;
end:
MakeUniqueList := proc(d,p)

   local ra,list,alpha,i;

   ra := rand(1..p-1);
#print(d);
   #list := [seq(0,i=1..d)];
   list := Array(1..d);

   list[1] := ra();
   for i from 2 to d do
      alpha := ra();
      while evalb(alpha in list) do
        alpha := ra();
      od;
      list[i] := alpha;
    od;

    return list;

end proc:
Interp1varKD := proc(B, X, P, a, b, e, d, p)

   local i, j, k, n,xPoints, yPoints, E ;

   n := nops(X);

   #get random unique points for interpolation
   xPoints := MakeUniqueList(d+1,p);
   #yPoints := [seq(0,i=0..d)];
   yPoints := Array(1..d+1);

   for i from 1 to d+1 do
         E := [xPoints[i],seq(e*(P[k+1]-a[k]*P[1]-b[k]) + a[k]*xPoints[i] + b[k] mod p,k=1..n-1)];
         yPoints[i] := B(E,p);
   od;

   return Interp(xPoints,yPoints,X[1]) mod p;

end proc:
CreateBBGcd := proc(F,X,degs,p)

   local i,j,k,n,ra,a,b,c,xPoints,yPoints,alpha,E,Interps,f0bar,f1bar,gamma0bar,BG;
   local r;

   r := nops(F);
   ra := rand(1..p-1);
   n := nops(X);

   a := [seq(ra(),i=2..n)]; b := [seq(ra(),i=2..n)]; c := [seq(ra(),i=3..r)];

  # print(a); print(b); print(c);

   #Interpolate polys in 1 variable
   for i from 1 to r do
      Interps[i] := Interp1varKD(F[i],X, [seq(0,i=1..n)],a,b,0,degs[i],p);
   od;
 #  print(Interps[1]); print(Interps[2]);

   f0bar := Interps[2] + add(c[i-2]*Interps[i],i=3..r) mod p;
   f1bar := Interps[1];

   #print(f0bar); print(f1bar);

   gamma0bar := Gcd(f0bar,f1bar) mod p; #print(gamma0bar);

   #Create the black-box to be returned
   BG := proc(P,q)

      local FF,A,B,C,D,Y,Gamma0bar,delta,Fbar1,F0bar1,F1bar1,Gammabar1,II,Gammas,e,F0bar,F1bar,Fbare;
      local ii,jj,kk,rr,nn,G,co;

      rr := r; nn := n; Y := X;
      A := a; B := b; C := c; D := [seq(degree(Interps[ii]),ii=1..rr)];
      FF := F;
      Gamma0bar := gamma0bar;
      delta := degree(Gamma0bar);

      #if the cardinality of Zq is not large enough to interpolate correctly
      if q < D[1]*D[2] + delta then return FAIL; fi;

      #do initial evaluation test for i1= 1
      for ii from 1 to rr do
        Fbar1[ii] := Interp1varKD(FF[ii],Y,P,A,B,1,D[ii],q);
      od;

      F0bar1 := Fbar1[2] + add(C[ii-2]*Fbar1[ii],ii=3..rr);
      F1bar1 := Fbar1[1];

      Gammabar1 := Gcd(  F0bar1, F1bar1 ) mod q;
      #print(Gammabar1);

      if degree(Gammabar1) = delta then
         return Eval(Gammabar1,X[1]=P[1]) mod q;
      elif degree(Gammabar1) > delta then

         II := [seq(0,ii=0..delta+1)]; Gammas := [Gamma0bar,seq(0,ii=1..delta+1)];

         for ii from 2 to delta+1 do

            do

              e := ra();
              while evalb(e in II) do
                 II[ii] := e;
              od;

              #perform bivariate interpolation for each black-box
              for jj from 1 to rr do
                 Fbare[jj] := Interp1varKD(FF[ii], X, P, A, B, e, D[jj], q);
              od;

              F0bar := Fbare[2] + add(C[jj-2]*Fbare[jj],jj=3..rr);
              F1bar := Fbare[1];

              Gammas[ii] := Gcd( F0bar, F1bar ) mod q;

              if degree(Gammas[ii]) < delta then return FAIL; fi;

            until degree(Gammas[ii]) = delta;


         od;

         G := x^delta;
         for jj from 0 to delta-1 do
            co := Interp(II,[seq(coeff(Gammas[kk],jj),kk=1..delta+1)],y) mod q;
            G := G + x^jj*co;
         od;

         return Eval(G,[x=P[1],y=1]) mod q;

      fi;

   end proc;


   return BG;

end proc:
with(ArrayTools):

local uniDegTime := time() - time();
local totDegTime := time() - time();
local checkTime := time() - time();
local gcdPTime := time() - time();
local chremTime := time() - time();
local evalTime := time() - time();
local createBBTime := time() - time();
local lambdaTime := time() - time();
local BenOrTime := time() - time();

KD := proc(BA,BB,X,T2,da,db,tdeg,pi)

   with(ArrayTools):

   local T,Tprev,i,n,primes,ra,beta,sigma,gamma,P,j,ai,bi,alpha,BR;
   local G,x,y,C,M,F,k,LM,Term,Terms,c,LM2;
   local r,D,ce,RR,Gbas,t,mm,v,gi,gprime,q,ra2,a2,q2,g2,xx,lambda,a,st,BG;

      n := nops(X); #print(X);
      x := X[1];

      ra := rand(1..pi-1);

      T := T2;
      beta := [seq(ra(),i=1..n)];
      gamma := [seq(ra(),i=1..n)];
      sigma := gamma;
      #print(sigma);
      G := Array(1..T);
      #G := [seq(0,i=0..2*T-1)];

      P := [2,seq(0,i=2..n)];
      for j from 2 to n do P[j] := nextprime(P[j-1]); od;

#print("before a");
      st := time();
      BG := CreateBBGcd([BA,BB],X,tdeg,pi);
      createBBTime := createBBTime + time() - st;

      j := 1;
      do
         #print("infinite loop");

         G := Concatenate(2,G,Array(1..T)); #print(G);
         st := time();
         while j <= 2*T do
            if j mod 100 = 0 then print(j); fi;
           #print("before");
            G[j] := BG(sigma,pi);
            #print("after");
            for k from 1 to n do sigma[k] := sigma[k]*P[k] mod pi; od;
            j := j+1;
         od;
         evalTime := evalTime + time() - st;

         st := time();
         lambda := BMEA([seq(G[i],i=1..2*T)],pi,z);
         lambdaTime := lambdaTime + time() - st;
         #print("In the T loop");
         T := 2*T;
      until degree(lambda,z) < T/2:
      T := T/2;
     #print("made it out"); #print(lambda);

     c := 0;
     while coeff(lambda,z,c) = 0 do
        c := c+1;
     od;
     lambda := Quo(lambda,z^c,z) mod pi; T := T-c;

     st := time();
     #print("start of Ben-Or Tiwari");
      RR := Roots(lambda) mod pi; #print(RR);
      k := nops(RR); #printf("k=%lld\n",k);
      Gbas := Array(1..k,fill=1);
      for i from 1 to k do
         t := RR[i][1];

         for j from 1 to n do
            while irem(t,P[j]) = 0 do
               Gbas[i] := Gbas[i]*X[j];
               t := t/P[j]; #print(t);
            od;
         od;
      od;
   #print("after roots");
BenOrTime := BenOrTime + time() - st;


      #solve for the coefficients of the gcd using a Vandermonde matrix
      mm := [seq(RR[i][1],i=1..k)];
      v := <seq(G[i],i=1..k)>;
      a := VandermondeSolve(v,mm,pi);

    # print("after Vandermonde");

      #assemble the gcd
      #gi := Array(1..k);
      #for i to k do gi[i] := a[i]*Gbas[i]; od;
      gi := add(a[i]*Gbas[i],i=1..k); #print("end of Ben-Or Tiwari");
#print(gi);
      #account for random gamma vector
      #print(gi);
      print(k);
      print(nops(gi));
      if k > 1 then
         #lcoeff(gi,X,'LM');
         #print(nops(gi)); print(k); print(gi);
         #Terms := [op(gi)]; #print(nops(Terms));
         Terms := Array(1..k);
         for j from 1 to k do Terms[j] := op(gi)[j]; od; #print(Terms);
         for j from 1 to k do
            lcoeff(Terms[j],X,'LM2'); #print(LM2);
            c := Eval(LM2,{seq(X[i]=gamma[i],i=1..n)}) mod pi;
            Terms[j] := Terms[j]/c mod pi;
         od; #print(Terms);
         gi := add(Terms);
      fi;

      #make the gcd monic in lex order
      gi := gi/lcoeff(gi,X) mod pi;
      #print(gi);

      return gi,k+1;

end proc:

BBMGCD:= proc(BA,BB,X)

   local i,j,k,gp,g,G,Gprime,gprime,p,q,M,n,da,db,a,b,ra,x1,alpha;
   local compare,LM,ra2,gamma,ge,aq,bq,T,st, alpha2,ta,tb,tdegg;


   M,G := 0,0;
   n := nops(X); p := nextprime(2^61);

   if n<1 then
      return "You need at least 1 variable to use this algorithm";
   fi;

   x1 := X[1];

   da := [seq(0,i=1..n)]; db := [seq(0,i=1..n)]; q := nextprime(10000000); ra := rand(1..q-1);
   alpha := [seq(ra(),i=1..n)];
   st := time();
   for i from 1 to 1 do
      a := UnivInterpNew(BA,X,alpha,i,q);
      b := UnivInterpNew(BB,X,alpha,i,q);
      da[i] := degree(a);
      db[i] := degree(b);
   od;
   uniDegTime := uniDegTime + time() - st;

   alpha := [seq(ra(),i=1..n)];
   alpha2 := [seq(ra(),i=1..n)];
   st := time();
      a := TotalDegNew(BA,X,alpha,alpha2,q);
      b := TotalDegNew(BB,X,alpha,alpha2,q);
   ta := degree(a);
   tb := degree(b);
   tdegg := degree(Gcd(a,b) mod q);
   totDegTime := totDegTime + time() - st;


   LM := 0;
   T := 2;
   #print(da[1]); print(db[1]);

   while true do

      p := nextprime(p); #if p > 1030 then return -1 fi;
      st := time();
      gp,T := KD(BA,BB,X,T,da[1],db[1],[ta,tb],p);
     gcdPTime := gcdPTime + time() - st;

      #return gp;

      #return -1;

      #gp := Gcd(A,B) mod p;
      #print (gp);
      if gp <> FAIL then


         compare := 0;
         if LM <> 0 then
            lcoeff(gp,X,'LM2'); #print(LM2);
            if LM <> LM2 then
               print("should terminate here");
               return FAIL;
               #for i from 1 to n do
                for i from n to 1 by -1 do
                  k := degree(LM2,X[i]) - degree(LM,X[i]);
                  #if degree(LM2,X[i])>0 and degree(LM,X[i])=0 or degree(LM2,X[i]) < degree(LM,X[i]) then
                   if k < 0 then
                     compare := 1; break;
                  #elif degree(LM2,X[i]) > degree(LM,X[i]) then
                   elif k > 0 then
                     compare := 2; break;
                  fi;
               od;
            fi;
         fi;

         #First prime
         if G=0 then
            gprime := iratrecon(gp,p);
            if gprime <> FAIL then gprime := denom(gprime)*gprime; fi;
            M,G,Gprime := p,gp,gprime;
            lcoeff(gp,X,'LM');

         elif compare = 1 then
            print("starting over");
            gprime := iratrecon(gp,p);
            if gprime <> FAIL then gprime := denom(gprime)*gprime; fi;
            M,G,Gprime := p,gp,gprime;
            lcoeff(gp,X,'LM');

         elif compare = 0 then

            print("normal case");

            st := time();
            g := chrem([G,gp],[M,p]); #print(g);
            M := M*p;
            gprime := iratrecon(g,M);
            if gprime <> FAIL then gprime := denom(gprime)*gprime; fi;
            chremTime := chremTime + time() - st;

           # print(gprime);

            if gprime <> FAIL and Gprime <> FAIL and gprime = Gprime then

                  q := nextprime(123456789); ra2 := rand(1..q-1);
                  do
                      gamma := [seq(ra2(),i=2..n)];
                      ge := Eval(gprime,{seq(X[i]=gamma[i-1],i=2..n)}) mod q;
                 until degree(ge,X[1]) = degree(gprime,X[1]);

                 st := time();
                 aq := ProbeInterpUni(BA,[0,op(gamma)],X[1],1,q);
                 bq := ProbeInterpUni(BB,[0,op(gamma)],X[1],1,q);
                  if Divide(aq,ge) mod q and Divide(bq,ge) mod q then
                    checkTime := checkTime + time() - st;
                    return gprime;
                 fi;
            fi;

            return FAIL;

            G,Gprime := g,gprime;

         elif compare = 2 then
            #print("ignore");
            #unlucky prime. do nothing

         fi;


      fi;

   od;


end proc: