rnd(P)=subst(P[1].pol,variable(P[1].pol),2)%3;

 rho(E,P,Q,V)=
 {
   my(X,a,b,h);
   [X,a,b] = V;
   h = rnd(X);
   if(h==0,[elladd(E,X,P),a+1,b],
      h==1,[elladd(E,X,Q),a,b+1],
           [ellmul(E,X,2),2*a,2*b]);
 }

 floyd(E,P,Q)=
 {
   my(X1,X2);
   X1=[P,1,0]; X2=[P,1,0];
   until(X1[1]==X2[1],
     X1=rho(E,P,Q,X1);
     X2=rho(E,P,Q,rho(E,P,Q,X2)));
   [X1,X2];
 }

 dlog(E,P,Q,o)=
 {
   my(X1,X2);
   [X1,X2]=floyd(E,P,Q);
   -(X1[3]-X2[3])/(X1[2]-X2[2])%o;
 }

 test(N)=
 {
   p =randomprime(N); a =ffgen(p);
   until(isprime(ellcard(E)), E = ellinit([1,random(p)],a));
   Q = ellgenerators(E)[1];
   o = ellgroup(E)[1];
   e = random(ellcard(E));
   P = ellmul(E,Q,e);
 }
 test(2^35);
 dlog(E,P,Q,o)
 ##
 elllog(E,P,Q,o)
 ##