localB:=proc(F,x,y,M,z)
local dx,dy,w,nF,R,B0,B1,B2;
    dy,dx,w:=degree(F,y),degree(F,x),degree(M,z);
    nF:=eval(numer(F),RootOf(M)=z);
    R:=(2*dy-1)!*dy^dy*((w+1)*(dx+1))^(2*dy-1)*norm(nF,infinity)^(2*dy-1);
    B0:=R*(norm(M,infinity)+1)^((w-1)(2*dy-2));
    B1:=(w+1)^((2*w-1)/2)*norm(M,infinity)^(w-1)*B0^w;
    B2:=w^w*(w+1)^((2*w-1)/2)*norm(M,infinity)^(2*w-1);
    max(denom(F),B1,B2);
end;

globalB:=proc(F,x,y,M,z)
local dx,dy,w,delta,nF,R1,B0,R2,R3,B3,B4;
    dx,dy,w:=degree(F,x),degree(F,y),degree(M,z);
    delta:=dx*(2*dy-1);
    nF:=eval(numer(F),RootOf(M)=z);
    R1:=(2*dy-1)!*dy^dy*((w+1)*(dx+1))^(2*dy-1)*norm(nF,infinity)^(2*dy-1);
    B0:=R1*(norm(M,infinity)+1)^((w-1)(2*dy-2));
    R2:=2^(delta+w)*(delta+1)^(1/2)*(w+1)^(7*w/2)*B0^(delta)*norm(M,infinity)^(4*delta);
    R3:=(2*delta-1)!*delta^delta*(w+1)^(2*delta-1)*R2^(2*delta-1);
    B3:=R3*(norm(M,infinity)+1)^((w-1)*(2*delta-2));
    B4:=(w+1)^((2*w-1)/2)*norm(M,infinity)^(w-1)*B3^w;
    max(denom(F),B4);
end;

DGP:=proc(F,x,y,M,z)
# This algorithm computes a bound B such that all prime number greater
# than B induce a local good reduction
    nextprime(localB(F,x,y,M,z));
end;

GDGP:=proc(F,x,y,M,z)
# This algorithm computes a bound B such that all prime number greater
# than B induce a global good reduction
    nextprime(globalB(F,x,y,M,z));
end;

FindP:=proc(B,d,eps)
# the function Find-p described in the paper
local K,C;
    if B<3 then return nextprime(d) fi;
    K:=2*ln(B)/(eps*ln(ln(B)))+2*d/ln(d);
    C:=ceil(evalf(max(2*d,K*ln(K)^2)));
    nextprime(rand(d+1..C)())
end;

MCGP:=proc(F,x,y,M,z,eps)
# the function Monte Carlo Good Prime described in the paper
local B;
    B:=evalf(localB(F,x,y,M,z));
    FindP(B,degree(F,y),eps/3)
end;

LVGP:=proc(F,x,y,M,z,tcRF)
# the function Las Vegas Good Prime described in the paper
# tcRF is the trailing coefficient of resultant_y(F,F_y) if we want a
# local reduction, and the resultant_x(S,S_x) where S is the squarefree part of
# resultant_y(F,F_y) if we want a global reduction
local dy,R,N1,N2,B,pi,p,i;
    dy:=degree(F,y);
    R:=numer(tcRF);
    N1:=abs(evala(Norm(R)));
    N2:=abs(discrim(M,z));
    B:=max(denom(F),N1,N2);
    pi:=denom(F)*N1*N2;
    p:=nextprime(dy);
    for i from 1 to 5 while pi mod p = 0 do
        p:=nextprime(p)
    od;
    while pi mod p = 0 do
        p:=FindP(B,dy,1/6)
    od;
    p;
end;


GMCGP:=proc(F,x,y,M,z,eps)
# the function Global Monte Carlo Good Prime described in the paper
local B;
    B:=evalf(globalB(F,x,y,M,z));
    FindP(B,degree(F,y),eps/4)
end;