discrete_log := proc(a,g,p,alpha0,beta0)
# Demonstration of the Pollard rho method for discrete log
# Copyright: Erich Kaltofen 2006
   local randomf, equl, pl, mu_lam;
      randomf := proc(x)
      # x is a list: [x_i,a,g,p,alpha,beta] where x_i = a^alpha g^beta mod p
         local x_i
              ,p
              ,g
              ,a
              ,expo_g,g_to_expo
              ,expo_a,a_to_expo
              ;
         x_i := x[1]; a := x[2]; g := x[3]; p := x[4];
         expo_g:=11; g_to_expo:=g &^ expo_g mod p;
         expo_a := 13; a_to_expo:=a &^ expo_a mod p;
         if x_i < p/3 then return([x_i*g_to_expo mod p,
                                   a,g,p,
                                   x[5],(x[6]+expo_g) mod (p-1)
                                  ]);
         fi;
         if x_i < 2*p/3 then return ([x_i^2 mod p,
                                     a,g,p,
                                     (2*x[5]) mod (p-1),
                                     (2*x[6]) mod (p-1)
                                     ]);
         fi;
         return([x_i*a_to_expo mod p,
                 a,g,p,
                 (x[5]+expo_a) mod (p-1), x[6]
                ]);
      end;

      equl := proc(x,y)
                local alpha1,alpha2,beta1,beta2;
                if x[1] = y[1] then
                   print("Collision: ", x,y);
                   alpha1 := x[5]; alpha2 := y[5];
                   beta1 := x[6]; beta2 := y[6];
                   if igcd(alpha1 - alpha2, p-1) = 1
                      then print("Success: ",
                                 (beta2-beta1)*(alpha1-alpha2)^(-1) mod (p-1));
                   fi;
                   return true
                else return false;
                fi;
              end;

      pl := printlevel; # printlevel := 3;
      mu_lam := cycle_indices([a^alpha0 * g^beta0 mod p, a,g,p, alpha0,beta0], randomf, equl);
      printlevel := pl;
      return (mu_lam);
end: