rnd(a)=subst(a.pol,variable(a.pol),2)%3
rho(V,a,b)=
{
  my(X,u,v,h);
  [X,u,v] = V;
  h = rnd(X);
  if(h==0,[a*X,u+1,v],
     h==1,[b*X,u,v+1],
          [X^2,2*u,2*v]);
}

floyd(a,b)=
{
  my(X1,X2);
  X1=[a^0,0,0]; X2=X1;
  until(X1[1]==X2[1],
    X1=rho(X1,a,b);
    X2=rho(rho(X2,a,b),a,b));
  [X1,X2];
}

dlog(a,b)=
{
  my(X1,X2);
  [X1,X2]=floyd(a,b);
  lift(-(X1[3]-X2[3])/(X1[2]-X2[2])*Mod(1,fforder(b)));
}

t=ffgen([2,19],'t);
b=ffprimroot(t);
e=random(fforder(b));
a=b^e;
e
dlog(a,b)
t=ffgen([2,31],'t);
b=ffprimroot(t);
e=random(fforder(b));
a=b^e;
e
dlog(a,b)