Discrete logarithm for small numbers

[Ismael Jimenez Calvo <ltqij04@pinar2.csic.es>]

	I miss a function for the discrete logarith in Z_p in the pari
package. I propose and include in this mail --as attachement-- a
function of complexity O( \sqrt p log( \sqrt p ) ) which can find in some
seconds a discrete logarithm in a group of order of about a million.
Somebody knows a faster function usefull in groups of moderate size?
Thank you very much,
		Ismael Jimenez Calvo

/* Discrete logarithm by "baby step big step" method.
   If y=logdiscr(x,b,p) then b^y \equiv x (mod p) */

{
logdiscr(x,b,p, pb=[[1,0]],i,j,jj,m,a,pi,pf,pp,z,fin)= 
if(gcd(p,b)!=1,return(-1),);
m=sqrtint(p-1)+1;a=1;pb=vector(m,X,[1,0]);
for(i=1,m-1,a=a*b%p;pb[i+1]=[a,i]);
pb=veclexsort(pb); /* en pb (b^i,i) */

a=lift(Mod(b,p)^(-m));z=x;
for(j=0,m-1,

 /*binary search*/
 if(j!=0,pi=1,pi=2);pf=m;i=-1;fin=0;
 if(pb[pf][1]==z,i=pb[pf][2];jj=j;break(),);
 while(!fin,pp=(pi+pf)\2;if(pp==pi,fin=1,);
  if(pb[pp][1]>z,pf=pp,
    if(pb[pp][1]<z,pi=pp,
      while(pp>2&&pb[pp-1][1]==pb[pp][1],pp=pp-1);i=pb[pp][2];jj=j;break(2)))
 );
 z=z*a%p;
);

if(i==-1,return(-1),z=(i+jj*m)%p;return(z));
}

/* Aplicacion del metodo de Pohlig-Hellman */
{
logdiscr1(x,b,p, n,ni,f,pi,bi,xi,y,z,zi)=
n=p-1;f=factor(n)~;
for(i=1,length(f),pi=f[1,i]^f[2,i];ni=(p-1)/pi;bi=lift(Mod(b,p)^ni);
 xi=lift(Mod(x,p)^ni);y=logdiscr(xi,bi,p);
 if(y==-1,return(-1),);zi=Mod(y,pi);
 if(i==1,z=zi,z=chinese(z,zi));
);
return(lift(z));
}


/* Metodo de Pollard-rho */
read("rand");
{
logdiscr2(x,b,p,s=[1,0,0,1,0,0],t,r,u,n,y,i )=
until(s[1]==s[4],t=rand(3);if(s[1]==1 && t==2,t=3,);
 if(t==0,s[1]=s[1]*x%p;s[2]=s[2]+1,
    if(t==1,s[1]=s[1]*s[1]%p;s[2]=2*s[2]%(p-1);s[3]=2*s[3]%(p-1),
      s[1]=s[1]*b%p;s[3]=s[3]+1));

 if(t==0,s[4]=s[4]*x%p;s[5]=s[5]+1,
    if(t==1,s[4]=s[4]*s[4]%p;s[5]=2*s[5]%(p-1);s[6]=2*s[6]%(p-1),
      s[4]=s[4]*b%p;s[6]=s[6]+1));

 if(s[1]==1 && t==2,t=3,); 
 if(t==0,s[4]=s[4]*x%p;s[5]=s[5]+1,
    if(t==1,s[4]=s[4]*s[4]%p;s[5]=2*s[5]%(p-1);s[6]=2*s[6]%(p-1),
      s[4]=s[4]*b%p;s[6]=s[6]+1));
);
r=(s[2]-s[5])%(p-1);u=(s[6]-s[3])%(p-1);d=gcd(r,p-1);r=r/d;u=u/d;n=(p-1)/d;
print([r,u,d]);
y=lift(Mod(u,n)*Mod(r,n)^(-1));
print("y_0=",y," n=",n);
for(i=1,d,if(lift(Mod(b,p)^y)==x,break(),y=y+n));print(y);return(y);
}