prdl:=function(a,b,p)
 q:= p-1;
 p1:= p div 3;
 p2:=2*p div 3;
 x:=<1,0,0>;
 y:=<b,1,0>;
 while x[1] ne y[1] do
   if x[1] le  p1 then
     x:=<b*x[1] mod p, (x[2]+1) mod q, x[3]>;
   elif x[1] le p2 then
     x:=<x[1]^2 mod p, 2*x[2] mod q, 2*x[3] mod q>;
   else
     x:=<a*x[1] mod p, x[2], (x[3]+1) mod q>;
   end if;
   if y[1] le  p1 then
     y:=<b*y[1] mod p, (y[2]+1) mod q, y[3]>;
   elif y[1] le p2 then
     y:=<y[1]^2 mod p, 2*y[2] mod q, 2*y[3] mod q>;
   else
     y:=<a*y[1] mod p, y[2], (y[3]+1) mod q>;
   end if;
   if y[1] le  p1 then
     y:=<b*y[1] mod p, (y[2]+1) mod q, y[3]>;
   elif y[1] le p2 then
     y:=<y[1]^2 mod p, 2*y[2] mod q, 2*y[3] mod q>;
   else
     y:=<a*y[1] mod p, y[2], (y[3]+1) mod q>;
   end if;
   x,y;
 end while;
 if IsDivisibleBy(2*(x[3]-y[3]),q) then
   return "failure";
 else
   "The discrete log of",a,"modulo",p,"relative to the base",b,"is";
   if IsEven(x[3]-y[3]) then
     ans:=((y[2]-x[2]) div 2)*InverseMod((x[3]-y[3]) div 2,q div 2) mod (q div 2);
     if IsEven(ans) then
      return (((LegendreSymbol(a,p)-1) div 2)*((p-1) div 2)+ans) mod (p-1);
     else
      return (((LegendreSymbol(a,p)+1) div 2)*((p-1) div 2)+ans) mod (p-1);
     end if;
  else
    return (y[2]-x[2])*InverseMod(x[3]-y[3],q) mod q;
  end if;
 end if;
end function;

choosesafeprime:=function(n)
repeat
p:=RandomPrime(n);
until IsProbablyPrime((p-1) div 2);
b:=1;
repeat
b:=b+1;
until Modexp(b,(p-1) div 2,p) ne 1;
return p,b;
end function;