discrete_log := proc # Demonstration of the Pollard rho method for discrete log # Copyright: Erich Kaltofen 2006, 2016 ( a # residue ,g # residue of order ord (initially primitive root) ,p # prime modulus ,alpha0 ,beta0 ,ord # g, g^2,..., g^ord equiv 1 contains the residue a ) local randomf ,equl ,mu_lam; randomf := proc(x) # x is a list: [x[1],a,g,p,alpha,beta] where x[1] = (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; expo_a := 13; if x_i < p/3 then g_to_expo:=g &^ expo_g mod p; return([x_i*g_to_expo mod p, a,g,p, x[5],(x[6]+expo_g) mod (ord) ]); fi; if x_i < 2*p/3 then return ([x_i^2 mod p, a,g,p, (2*x[5]) mod (ord), (2*x[6]) mod (ord) ]); fi; a_to_expo:=a &^ expo_a mod p; return([x_i*a_to_expo mod p, a,g,p, (x[5]+expo_a) mod (ord), x[6] ]); end; # randomf equl := proc(x,y) # returns false or # [true, dlog] or [true, false], the latter when the # collisions did not produce a result local alpha1,alpha2,beta1,beta2 ,d ,ord_over_d ,a1_minus_a2_over_d ,dlog_a, dlog_a_to_d, dlog_a_adj ,G # new base g^(ord/d) [see below] ,G_to_i ,a_adj ,i ; if x[1] = y[1] then if printlevel >= 0 then print("Collision: ", x,y); fi; alpha1 := x[5]; alpha2 := y[5]; beta1 := x[6]; beta2 := y[6]; d := igcd(alpha1 - alpha2, ord); if d = 1 # check if lucky then # Eureka! # a^(alpha1-alpha2) equiv g^(beta2 - beta1) # and (alpha1 - alpha2) has an inverse mod ord dlog_a := (beta2-beta1)*(alpha1-alpha2)^(-1) mod (ord); if printlevel >= 0 then print("Success: " ,dlog_a ,ifactor(alpha1 - alpha2) ); fi; return([true,dlog_a]); fi; # if d = 1 # case d > 1 if printlevel >= 0 then print("Failure:" ,ifactor(alpha1 - alpha2) ,"ord=",ifactor(ord), "gcd = ", d); fi; if d = ord then if printlevel >= 0 then print("No order reduction possible"); fi; return([true, false]); # indicate failure fi; # no order reduction possible # reduce the order ord_over_d := ord/d; a1_minus_a2_over_d := (alpha1 - alpha2)/d; # (a^d)^( (alpha1 - alpha2)/d ) equiv g^(beta1 - beta1) dlog_a_to_d := ((beta2-beta1) *( a1_minus_a2_over_d^(-1) mod ord_over_d ) ) mod ord; # a^d equiv g ^ dlog_a_to_d, with dlog_a_to_d divisible by d # (all d-th powers have discrete log divisible by d, which # is a divisor of the order of g) a_adj := a * g &^ (-iquo(dlog_a_to_d, d)) mod p; # Since ( a * g^(-dlog_a_to_d/d) )^d equiv 1, it has order # a divisor of d, so is in the set { g^(i*ord/d) } G := g &^ ord_over_d mod p; # avoid recursive call for small d if d <= 20 then # check each exponent print("Iterative continuation with" ,a_adj ,G ,p ,"ord =",d ); G_to_i := 1; for i from 0 to d-1 do if G_to_i = a_adj then dlog_a_adj := i; break; fi; G_to_i := G * G_to_i mod p; od; else # call discrete log recursively print("Recursive continuation with" ,a_adj ,G ,p ,"init = ", alpha0, beta0 ,"ord =",d ); dlog_a_adj := discrete_log(a_adj ,g &^ ord_over_d mod p ,p,alpha0,beta0,d); if dlog_a_adj = false # failure in recursive call then # try additional calls with different starting alpha0, beta0 for i from 2 to 5 while dlog_a_adj = false do if printlevel >= 1 then print("Try different initial values in recursive call"); fi; dlog_a_adj := discrete_log(a * g&^(-iquo(dlog_a_to_d, d)) mod p ,g &^ ord_over_d mod p ,p ,(1+alpha0 &^ 3) mod d ,(1+beta0 &^ 3) mod d ,d); od; if dlog_a_adj = false # failure in all other calls recursive call then return [true,false]; fi; fi; # dlog_a_adj = false # Eureka again: recursive call succeeded fi; # d <= 20 return [true, (iquo(dlog_a_to_d, d) + dlog_a_adj * ord_over_d) mod ord # note: a_adj equiv ( g^(ord/d) )^dlog_a_adj # so a * g^(-dlog_a_to_d/d) equiv g^(dloga_a_adj * ord/d) ]; else return false; # no collision fi; # x[1] = y[1] end; mu_lam := cycle_indices([a^alpha0 * g^beta0 mod p ,a, g, p ,alpha0, beta0 ] ,randomf, equl ,'return_mu_lambda'=false ,'additional_data_var'='DLOG' ); return (DLOG); end: # discrete log