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:

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 RatreconTime := time() - time();
local lambdaTime := time() - time();
local BenOrTime := time() - time();

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:
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:
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;
     do
        alphak := ra();
     until not evalb(alphak in alpha);
     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);

   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:

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:

KY := proc(BA,BB,X,T3,da1,db1,dg1,p)

   local i,j,k,ra,n,beta,sigma;
   local T,P,alpha,alphak,BR,x,T1,T2,y,E,M,F,G,C,D,v,ce,lambda,c,RR,Gbas,t,mm,V,a,gp,st;

   ra := rand(1..p-1);
   n := nops(X);
   x := X[1];

   beta := [seq(ra(),i=1..n)];
   sigma := [seq(1,i=1..n)]; T := T3;
   P := [seq(0,i=1..n)];
   P[1] := 2;
   for i from 2 to n do P[i] := nextprime(P[i-1]); od; #print(P);

   alpha := Array(1..da1+db1-2*dg1+1);
   for i from 1 to da1+db1-2*dg1+1 do
      alphak := ra();
      while evalb(alphak in alpha) do
         alphak := ra();
      od;
      alpha[i] := alphak;
   od;

   BR := proc(Y)
         local a,b;

         a := BA(Y,p);
         b := BB(Y,p);

         if b=0 then return FAIL; fi;

         return a/b mod p;

      end proc;


   v := Array(1..T);print(T);

   j := 1;

   do
      v := Concatenate(2,v,Array(1..T)); #print(v);
      #v := [op(v),seq(0,i=1..T)];
      #if T=2 then T1 := 0; else T1 := T; fi;
      #T2 := 2*T-1; print(T1); print(T2);
      while j <= 2*T do

      #if j mod 100 = 0 then print(j); fi;

         y := Array(1..da1+db1-2*dg1+1);
         st := time();
         for k from 1 to da1+db1-2*dg1+1 do
            E := [alpha[k],seq(beta[i]*(alpha[k]-sigma[1])+sigma[i] mod p,i=2..n)];
            y[k] := BR(E);
         od;
         evalTime := evalTime + time() - st;

         st := time();
         M := mul(x-alpha[k],k=1..da1+db1-2*dg1+1);
         F := Interp(alpha,y,x) mod p;
         G := Ratrecon(F,M,x) mod p;
         C := numer(G); D := denom(G); #print(G);
         ce := Eval(C,x=sigma[1]) mod p; #print(ce);
         v[j] := BA(sigma,p)/ce mod p;

         RatreconTime := RatreconTime + time() - st;
         #print("after v");
         for k from 1 to n do sigma[k] := sigma[k]*P[k] mod p; od;
         j := j+1;
      od;
      #print(v);
      st := time();
      lambda := BMEA(convert(v,list),p,z); #print(lambda);
      lambdaTime := lambdaTime + time() - st;
      #print(degree(lambda));
      T := 2*T; #print(T);
   until degree(lambda,z) < T/2;

   T := T/2;

   #print(lambda);

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

     st := time();

    #print("start of Ben-Or Tiwari");
      RR := Roots(lambda) mod p; #print(RR);
      k := nops(RR); #print(RR);
      Gbas := Array(1..k,fill=1);
      for i from 1 to k do
         t := RR[i][1]; #print(t);
         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);
               #print(Gbas[i]);
            od;
         od;
      od;
   #print("after roots");

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

      #print("after Vandermonde"); print(a); print(Gbas);

      #assemble the gcd
      gp := add(a[i]*Gbas[i],i=1..k);

      BenOrTime := BenOrTime + time() - st;
      #print("end of Ben-Or Tiwari");
      #print(gp);

      return gp,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,dg,alpha2,tdega,tdegb,tdegg;


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

   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)]; dg := [seq(0,i=1..n)]; q := nextprime(2^60); ra := rand(1..q-1);
   alpha := [seq(ra(),i=1..n)];

   st := time();
    for i from 1 to n do
       a := UnivInterpNew(BA,X,alpha,i,q); #print(Factor(Expand(a) mod q) mod q);
       b := UnivInterpNew(BB,X,alpha,i,q); #print(Factor(Expand(b) mod q) mod q);
       da[i] := degree(a);
       db[i] := degree(b);
       dg[i] := degree( Gcd(a,b) mod q);
   od;
   uniDegTime := uniDegTime + time() - st;

  # print(da,db);

   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);
   tdega := degree(a);
   tdegb := degree(b);
   tdegg := degree(Gcd(a,b) mod q);
   totDegTime := totDegTime + time() - st;

  #print(tdega,tdegb,tdegg);

   LM := 0;
   T := 2;

   while true do

      p := nextprime(p); #if p > 1030 then return -1 fi;
      st := time();
      gp,T := KY(BA,BB,X,T, tdega,tdegb,tdegg,p);
      gcdPTime := gcdPTime + time() - st;

      #print(Factor(gp) mod p);

      #make the leading coefficient monic
      gp := gp/lcoeff(gp,X) mod p;


      if gp <> FAIL then

         compare := 0;
         if LM <> 0 then
            lcoeff(gp,X,'LM2'); #print(LM2);
            if LM <> LM2 then
               #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:

printTimings := proc();

print("modular timings");
printf("Eval time=%8.2fs\n",(evalTime));
printf("Ratrecon time=%8.2fs\n",(RatreconTime));
printf("lambda time=%8.2fs\n",(lambdaTime));
printf("BenOr time=%8.2fs\n",(BenOrTime));

print("main timings");
printf("uni deg time=%8.2fs\n",(uniDegTime));
printf("tot deg time=%8.2fs\n",(totDegTime));
printf("gcd mod p time=%8.2fs\n",(gcdPTime));
printf("chrem + iratrecon time=%8.2fs\n",(chremTime));
printf("check time=%8.2fs\n",(checkTime));


end proc;