#Author := Ayoola Jinadu
# If you find a bug in the code please email me at ajinadu@sfu.ca

SameDeg := proc(A,s)local k,i,n;
  k := degree(A[1],s);
  n := nops(A);
  for i to n do
     if IsEqual(k,degree(A[i],s))=false then return (false,i); fi;
  od;
return k,-1;
end:


Same := proc(A)
local k,i;
 k := A[1];
 for i to numelems(A) do
     if k <> A[i] then return false; fi;
 od;
 return true;
end:

KronInv := proc(m,r,X) local i,k,A;
  A := btfactor(m,r);
  k := 1;
  for i to numelems(X) do
    if A[i] <> 0 then k := X[i]^A[i]*k; fi;
  od;
return k;
end:

btfactor := proc(m::posint,W::list(posint))
  local n,P,x,i,k,q;
  n := numelems(W);
  x := m;
  for i from 1 to n  do
     P[i] := irem(x,W[i],'q');
     x := q;
  od;
  [seq( P[i],i=1..n),x];
end:

BMEA := proc(a,p,x) local m,n,R0,R1,V0,V1,i;
  n := iquo(numelems(a) , 2 );
  m := 2*n-1;
  R0 := x^(2*n);
  R1 := add(a[m+1-i]*x^i,i=0..m);
  V0 := 0;
  V1 := 1;
  while n <= degree(R1,x) do
     R0,R1 := R1,Rem(R0,R1,x,'Q') mod p;
     V0,V1 := V1,Expand(V0-Q*V1) mod p;
  od;
  i := 1/lcoeff(V1,x) mod p;
  i*V1 mod p;
end:

VSolve := proc(v,u,p)
local t,i,j,M,x,a,q,r,s,Q,temp;
  t := numelems(v);
  M := modp1(ConvertIn(1,x),p);
  for i to t do
    temp := modp1(ConvertIn(x-u[i],x),p);
    M := modp1(:-Multiply(temp,M),p);
  od;
  a := Vector(t);
  for j to t do
    Q := modp1(Quo(M,ConvertIn(x-u[j],x)),p);
    r := 1/modp1(Eval(Q,u[j]),p) mod p;
    Q := modp1(ConvertOut(Q),p); #dense list#
    s := add(v[i]*Q[i],i=1..t) mod p;
    a[j] := r*s mod p;
  od;
  a := convert(a,list);
  return a;
end:

Invert := proc(A,r,X,A1,p,Mb,u);
 local T,G,f1,ij,e,Y,j,k,i,t,M1,Ma;
  if u <> 0 then
     Ma :=  [seq(1/Mb[i] mod p, i=1..nops(Mb)) ];
     M1 := [seq(i&^u mod p, i in Ma ) ];
  fi;
 f1 := 0;
 for ij to numelems(A) do
   e := KronInv(A[ij],r,X);
   if u <> 0 then f1 := f1+e*(A1[ij]*M1[ij]) mod p ;
     else f1 := f1+e*(A1[ij]) mod p
   fi;
 od;
return f1;
end:



InvBot := proc(CN,Rx,X,p,u,LamA,ta)
  local Lam1,f1,R1,M1,r1,E1,i,Mon1,A1,n;
  Lam1 := LamA; n := nops(X)-1;
  if Same(CN) = true and CN[1] <> 1 then f1 := CN[1];
      elif Same(CN) = true and CN[1] = 1 then f1 := 1;
        else
             Lam1 := modp1(ConvertIn(Lam1,s),p);
             R1 := modp1(Roots(Lam1), p);
             M1 := [ seq(r1[1],r1=R1)];
             E1 := [ seq( ModularLog(M1[i],ta,p,method=optimal), i=1..nops(M1)) ];
             Mon1 := [ seq( CN[i], i=1..nops(M1) )];
             A1 := VSolve(Mon1,M1,p);
             if n > 1 then
                f1 := Invert(E1,Rx,X[2..n+1],A1,p,M1,u);
               elif n=1 and u <> 0 then f1 := (A1[1]/M1[1]&^u)*X[2]^E1[1] mod p;
             fi;
      fi;
   return f1;
end: