# This gcd method assumes that x and y are positive def mygcd x,y g,r = y,x%y g,r = r,g%r while r>0 return g end # Verbose gcd method, again assumes that x and y are positive # Returns gcd, number of divisions, worst-case estimate for number of divisions def vgcd x,y g,r = y,x%y count = 1 print g, "\n" while r>0 g,r = r,g%r count = count + 1 print g, "\n" end return g, count, Math.log(Math.sqrt(5) * y)/Math.log((1+Math.sqrt(5))/2)+1 end # Extended gcd also returns c,d such that cx+dy=g def extgcd x,y g,q,r,a,b,c,d = y,x/y,x%y,1,0,0,1 while r>0 a,b,c,d = c,d,a-q*c,b-q*d q,s = g/r,g%r g,r = r,s end return g,c,d,c*x+d*y end # Exponentiate b^e % n def raise b,e,n x,t,y = b,e,1 while t>0 if (t%2==0) x,t=(x*x)%n,t/2 else y,t=(x*y)%n,t-1 end end return y end # Solve the equation ax = b mod n # Assumes n is positive def axbmodn a,b,n a,b = -a,-b if a<0 g = mygcd a,n if b%g != 0 return false, nil else at = a/g bt = b/g nt = n/g c = extgcd(at,nt)[1] x0 = (c*bt)%nt ans = Array.new(g) {|i| x0 + i*nt} return true, ans end end # To prepare for Miller-Rabin, factor n-1 def factornminusone n s,m=0,n-1 s,m=s+1,m>>1 while m&1==0 return s,m end # Miller-Rabin # Input: base b, candidate prime n such that n-1 = 2^s m # Output: true if n could be prime, false if n is definitely composite def mrprime b,n,s,m b=raise b,m,n return true if b==1 (s-1).times {return true if b==n-1 b=(b*b)%n } return true if b==n-1 return false end # Find first integer that Miller-Rabin thinks is prime # for base 2 and 3 but not base 5. def findfirst235 nmax=10**6 nmax.times {|n| nn=2*n+13 print "Starting a new block of 10000\n" if nn%10000==1 s,m=factornminusone nn if (mrprime 2,nn,s,m) and (mrprime 3,nn,s,m) and not (mrprime 5,nn,s,m) print nn, " passes the Miller-Rabin test for b=2 and b=3 but not b=5.\n" return nn end } end # Check the probabilistic primality of n: # First do trial division by some small primes # Second do the Fermat test for some bases # Third do the Miller-Rabin test for some bases def checkprimality n smallprimestable=[2,3,5,7,11,13,17,19] fbasestable=[2,3] #,5,7,11,13,17,19,23,29,31,37,41,43,47] mrbasestable=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47] smallprimestable.each {|p| return false if n%p==0} fbasestable.each {|b| return false if raise(b,n-1,n)!=1} s,m=factornminusone n mrbasestable.each {|b| return false if mrprime(b,n,s,m)==false} return true end # Find a big probable prime at least n def findbigprobprime n print "Searching for a probable prime at least ", n, ".\n" n += 1 if n%2==0 while checkprimality(n)==false print "~" n +=2 end print "\nThe first one found is ", n, ".\n" return n end # Lucas-Pocklington-Lehmer test # Use the base b to investigate the primality of n=pu+1 # where p is prime and 1\le u
1 print "Found a factor: ", g, ".\n" return g end pollardprimestable.each {|p| ell=(Math.log(n)/Math.log(p)).ceil print "Computed l = ", ell, ".\n" b=raise(b,p**ell,n) print "Computed b = ", b, ".\n" g=mygcd(b-1,n) if g==n print "The algorithm has failed.\n" return false elsif g>1 print "Found a factor: ", g, ".\n" return g end } print "The algorithm has not found a factor.\n" return false end # p is prime, g generates (Z/pZ)^*, b=a^l. # Find l, the discrete log of b to base a. def babygiant a, b, p m = Math.sqrt(p).ceil am = raise a,m,p # Make a list with entries [a^(mj) % p, j] for j=0,...,m-1, # and then sort it by the first elements of its entries list = [[1,0]] (m-1).times {|j| list<<[(list[j][0]*am)%p,j+1] } print "babygiant has made the first list" list = list.sort_by {|x| x[0]} print "babygiant has sorted the first list" # Now find a pair [(ba^(-i)) % p, i] whose first entry matches an [a^(mj) % p, j] in the list. # The code here does not exploit the fact that the list has been sorted. ainv = extgcd(a,p)[1] candidate = [b,0] candidate = [(candidate[0]*ainv)%p,candidate[1]+1] while list.assoc(candidate[0])==nil # So now b=a^(i+mj). print "babygiant has found candidate" return candidate[1]+m*list.assoc(candidate[0])[1] end # Auxiliary function for discrete logarithm Pollard rho method def f x,a,b,ea,eb,p if x%3==0 return (a*x)%p,(ea+1)%(p-1),eb elsif x%3==1 return (b*x)%p,ea,(eb+1)%(p-1) else return (x*x)%p,(2*ea)%(p-1),(2*eb)%(p-1) end end def pollardrho2 a, b, p x,eax,ebx = 1,0,0 y,eay,eby = f(x,a,b,eax,ebx,p) result = [false, nil] while result == [false,nil] while x != y x,eax,ebx = f(x,a,b,eax,ebx,p) y,eay,eby = f(y,a,b,eay,eby,p) y,eay,eby = f(y,a,b,eay,eby,p) end print "Found (x,ea,eb) = (", x, ", ", eax, ", ", ebx, ")\n" print "Found (y,ea,eb) = (", y, ", ", eay, ", ", eby, ")\n" # So now a^(eax + l ebx) = a^(eay + l eby), # i.e., l(ebx-eby) = eay-eax mod p-1 print "Need to solve the equation ", ebx-eby, " l = ", eay-eax, " mod ", p-1, "\n" result = axbmodn(ebx-eby,eay-eax,p-1) # print "Result: ", result, ".\n" if(result==[false,nil]) x,eax,ebx = f(x,a,b,eax,ebx,p) y,eay,eby = f(y,a,b,eay,eby,p) y,eay,eby = f(y,a,b,eay,eby,p) end end return result[1][0] end # Elliptic curve routines # The curve is y^2 = x^3 + ax + b mod n # Assuming x,y>=0 to make sure the gcd routine behaves def ecadd x1,y1,x2,y2,a,b,n if ((y1*y1 - x1*x1*x1 - a*x1 - b)%n != 0) or ((y2*y2 - x2*x2*x2 - a*x2 - b)%n != 0) print "e.c. add routine: a point does not lie on the curve.\n" return false end g=extgcd((x2-x1)%n,n)[0..1] print "g = ", g[0], ", ", g[1], ".\n" if g[0] != 1 print "gcd(x2-x1,n) = ", g[0], ".\n" return false else lambda = ((y2-y1)*(g[1]))%n newx = (lambda*lambda-x1-x2)%n newy = (-lambda*(newx-x1)-y1)%n print "Sum is ", newx, ", ", newy, ".\n" return true end end def ecdouble x,y,a,b,n if (y*y - x*x*x - a*x - b)%n != 0 print "e.c. double routine: the point does not lie on the curve.\n" return false end g=extgcd(2*y,n)[0..1] print "g = ", g[0], ", ", g[1], ".\n" if g[0] != 1 print "gcd(2y,n) = ", g[0], ".\n" return false else lambda = ((3*x*x+a)*(g[1]))%n newx = (lambda*lambda-2*x)%n newy = (-lambda*(newx-x)-y)%n print "Double is ", newx, ", ", newy, ".\n" return true end end