# This version of BDP uses modp routines in Zp[x] to speedup the implementation
BDP := proc(A,B,C,x,y,dy,p)
local dx,r,Sig,Tau,sig,tau,M,dm,beta,Ab,Bb,G,S,T,Cb,failed,dsig,dtau,SigT,TauT,k;
# Given A,B,C in Zp[x,y] with deg(A,x)+deg(B,x)>deg(C,x) and gcd(A,B)=1,
# solve sigma A + tau B = C in Zp[x,y] for sigma and tau in Zp[x,y] such that
# deg(sigma,y)+deg(B,y)<=dy, deg(tau,y)+deg(B,y)<=dy and deg(C,y)<=dy, if they exist.
# Algorithm 1 on input of $a_j,f_{j-1},g_{j-1}$ solves the MDP
# $\sigma g_{j-1} + \tau f_{j-1} = c_i$.
# Here if $d_2 = deg(a_j,x_2)$ we have $\deg(sigma g_{j-1}) \le d_2$,
# and $\deg(tau f_{j-1},x_2) \le d_2$ and $\deg(c_i,x_2) \le d_2$.
# So here for BDP we may use dy = d2 = $\deg(a,x_2)$.
  dx := degree(A,x)+degree(B,x);
  if degree(C,x)>=dx then error "C is wrong"; fi;
  r := rand(0..p);
  failed := false;
  dsig := degree(B,x)-1;
  dtau := degree(A,x)-1;
  ET := 0;
  GT := 0;
  A2 := modp2(ConvertIn(A,x,y),p); A2 := [op(A2)];
  B2 := modp2(ConvertIn(B,x,y),p); B2 := [op(B2)];
  C2 := modp2(ConvertIn(C,x,y),p); C2 := [op(C2)];
  one := modp1(ConvertIn(1,x),p);
  Sig := Array(1..dy+1);
  Tau := Array(1..dy+1);
  M := Array(1..dy+1);
  for k to dy+1 do
    do beta := r(); until not member(beta,M);
st := time():
    #Ab := Eval(A,y=beta) mod p;
    Ab := modp1(Eval(A2,beta),p);
    #Bb := Eval(B,y=beta) mod p;
    Bb := modp1(Eval(B2,beta),p);
    #Cb := Eval(C,y=beta) mod p;
    Cb := modp1(Eval(C2,beta),p);
ET := ET+time()-st;
st := time():
    #G := Gcdex( Ab, Bb, x, 'S', 'T' ) mod p;
    G := modp1( Gcdex(Ab,Bb,'S','T'),p);
GT := GT+time()-st;
    if G <> one then # beta is unlucky
        if failed then return FAIL else failed := true fi;
        print(UNLUCKY);
        k := k-1;
        next; # try one more time
    fi;
    # Solve sig Bb + tau Ab = Cb in Zp[x] with deg(sig)<deg(Bb) and deg(tau)<deg(Ab)
    #sig := Rem( Cb*S, Bb, x ) mod p;
    sig := modp1( Rem( Multiply(Cb,S), Bb ), p );
    #tau := Quo( Cb-sig*Ab, Bb, x ) mod p;
    tmp := modp1( Subtract( Cb, Multiply(sig,Ab)), p );
    tau := modp1( Quo(tmp,Bb), p );
    #sig := modp1(ConvertOut(sig,x),p);
    Sig[k] := sig;
    #tau := modp1(ConvertOut(tau,x),p);
    Tau[k] := tau;
    M[k] := beta;
  od;
st := time();
  #Sig := Interp( M, [Sig], y ) mod p;
  M := [seq(M[i],i=1..dy+1)];
  SigT := [seq( [seq(modp1(Coeff(Sig[i],j),p), i=1..dy+1)], j=0..dsig )];
  SigT := [seq( modp1(Interp(M,sig,y),p), sig=SigT )];
  SigT := [seq( modp1(ConvertOut(sig,y),p), sig=SigT )];
  Sig := add( SigT[i+1]*x^i, i=0..dsig );
  #Tau := Interp( M, [Tau], y ) mod p;
  TauT := [seq( [seq(modp1(Coeff(Tau[i],j),p), i=1..dy+1)], j=0..dtau )];
  TauT := [seq( modp1(Interp(M,tau,y),p), tau=TauT )];
  TauT := [seq( modp1(ConvertOut(tau,y),p), tau=TauT )];
  Tau := add( TauT[i+1]*x^i, i=0..dtau );
  if degree(Sig,y)+degree(A,y) > dy or degree(Tau,y)+degree(B,y) > dy then return FAIL fi;
  # We have M(y)|(C-Sig A - Tau B) and deg(M,y) > deg(C-Sig A - Tau B,y) ==> Sig A + Tau B = C
print(EVAL=ET,INTERP=time()-st,GCDEX=GT);
  return Sig,Tau;
end:

A := x^3+y^2*x+1;
B := x^2+y+3;
sig := y*x+3*y;
tau := x^2+y^2+3;
C := expand( sig*A+tau*B );
p := 101;
BDP( A,B,C,x,y,3,p);
BDP( A,B,C,x,y,2,p);

A := x^5+randpoly([x,y],dense,degree=5) mod p:
B := x^5+randpoly([x,y],dense,degree=5) mod p:
sig := randpoly([x,y],dense,degree=4) mod p:
tau := randpoly([x,y],dense,degree=4) mod p:
C := expand(sig*A+tau*B):
Sig,Tau := BDP( A,B,C,x,y,9,p);
ZERO( Expand(Sig-sig) mod p );
ZERO( Expand(Tau-tau) mod p );
BDP( A,B,C,x,y,8,p);

# random inputs will likely not have a polynomial solution
BDP( x^2+y^2,x^2+y*x+1,y*x+y,x,y,5,p);


d := 200;
p := 2^31-1;
A := x^d+randpoly([x,y],dense,degree=d-1) mod p:
B := x^d+randpoly([x,y],dense,degree=d-1) mod p:
sig := randpoly([x,y],dense,degree=d-1) mod p:
tau := randpoly([x,y],dense,degree=d-1) mod p:
C := expand(sig*A+tau*B) mod p:
print(degree(A,x),degree(sig,x),degree(B,x),degree(tau,x));
st := time():
Sig,Tau := BDP( A,B,C,x,y,2*d-1,p):
time()-st;
ZERO( Expand(Sig-sig) mod p );
ZERO( Expand(Tau-tau) mod p );
#BDP( A,B,C,x,y,2*d-3,p);