In [1]:
def shanks(p,g,h): # Solves g^x=h in F_p using the Shanks algorithm.
n=floor(p**0.5)+1
list_a=[]
i=0
temp=1
print 'Producing the list of powers of g...',
while i<=n-1:
list_a.append( [temp,i] )
temp=(temp*g)%p
i=i+1
print 'done.' # list_a now contains [g^0,0],[g^1,1],..,[g^{n-1},n-1]
# We are keeping track of both the power AND the exponent to make our life
# easier later
print 'Producing the list of h*g^{-jn}...',
list_b=[]
ginv=fastpower(g,p-2,p)
u=fastpower(ginv,n,p)
i=0
temp=h
while i<=n-1:
list_b.append( temp)
temp=(temp*u)%p
i=i+1
print 'done.' # list_b now contains hg^0,h*g^{-n},h*g^{-2n},...,h*g^{-(n-1)n}
print 'Starting collision algorithm for two lists of length', len(list_a)
l=collision(list_a,list_b) # Do the collision algorithm. l[0]=i will be the index in the first list
# and l[1]=j will be the index in the second list
if l=='F':
print "\n There is no solution"
return
return l[0]+n*l[1] # The solution is x=i+n*j.
In [2]:
def fastpower(x,e,p): # Compute x^e in Z/p using the fast powering algorithm
temp=1
list=[]
if e==0:
return 1
while e>0:
if e%2==0:
list.append(2)
e=e/2
continue
list.append(1)
e=e-1
continue
a=list.pop()
temp=x
while len(list)>0:
a=list.pop()
if a==1:
temp=(temp*x)%p
continue
temp=(temp*temp)%p
return temp
In [3]:
fastpower(2,9,1024)
Out[3]:
In [4]:
def collision(l,m): #Given two lists l=[(a0,b0),(a1,b1),...] and m=[x0,x1,x2,...]
#having the same length, see if some ai matches an xj.
#If no, return 'F' for 'Fail'.
i=0 #If yes, return the pair [i,j] giving the indices for the matching elements
print 'Sorting the first list...',
l.sort(key=sortFirst) #First sort the list l so that a0<=a1<=a2<=...
print 'done.' #This can take a while for very long lists.
print 'Testing elements of the second list...',
i=0 #Now test elements of m one by one to see if they match one of the ai terms.
n=len(m)
while i<=n-1:
if i%100==0:
print '.',
x=m[i]
u=is_in(l,x) #The function is_in decides if x matches one of the ai
if u=='F':
i=i+1
continue
print 'done'
return [u,i]
print 'done.'
return 'F'
In [5]:
def is_in(l,x): #Given l=[[a0,b0],[a1,b1],...] and x, decide if x equals one of the ai.
#Use the method of "compare to middle term" and divide the problem in half.
if len(l)==0: #If no, return 'F'. If yes, return the bi.
return 'F' #Assumes that a0<=a1<=a2<=... from the beginning.
first=0
last=len(l)-1 #Set up pointers "first" and "last", originally at the start and end of the list l.
while first<last:
middle=first+((last-first)//2) #Go to roughly the middle term
u=l[middle][0] #and pull out the a-value
if u==x: #If it matches our x, return the corresponding b-value.
return l[middle][1]
if u<x: #If x is bigger, reset the pointers to look at the last half of the list
first=middle+1 #and repeat
continue
last=middle-1 #If x is smaller, reset the pointers to look at the first half of the list.
#If we come out of the loop, we never found our x.
if first>last: #Either we exhausted all elements, or we narrowed it down to exactly one.
return 'F'
u=l[first]
if u[0]==x: #Do one last check to see if we have the x.
return u[1]
return 'F'
In [6]:
def findprime(lower,upper,tries):
i=1
while i<=tries:
a=randint(lower,upper)
if is_prime(a):
return a
i=i+1
print 'None found'
return 0
In [7]:
def order(x,p):
i=1
temp=x
while i<=p:
# if i%100==0:
# print '\r testing ',i,
if temp==1:
return i
temp=(temp*x)%p
i=i+1
print 'Error'
return
In [22]:
a
Out[22]:
In [8]:
def find_prim_gen(p): # Find a primitive generator for F_p
x=2 # Warning: For big primes (9 or more digits) this can take a while.
while x<=p-1:
print 'Testing ',x,'...',
u=order(x,p)
print 'Order is ',u,'...',
if u==p-1:
print 'So found a generator'
return x
print 'So not a generator'
x=x+1
print 'Error'
return
In [9]:
find_prim_gen(a)
In [ ]:
In [ ]:
In [10]:
def dlp_naive(p,g,h): # Solve g^x=h in Z/p, given g and h. Uses the brute force algorithm.
i=0 # i keeps track of the exponent we are trying
temp=1 #Start with g^0
while i<=p:
if temp==h: #if we found it, return the exponent
return i
temp=(temp*g)%p # Otherwise multiply by one more copy of g and increase i.
i=i+1
print 'Error: no solution found'
return
In [ ]:
In [ ]:
In [12]:
dlp_naive(23,5,19)
Out[12]:
In [13]:
def ph_estimate(n):
i=1
save=0
while i<=n:
b=findprime(10^20,10^21,35)
if b==0:
continue
l=factor(b-1)
count=0
for a in l:
count=count+a[1]*a[0]
save=save+(b+0.0)/count
i=i+1
return (save+0.0)/n
In [14]:
def miller_rabin_witness(n,a): # Determines if a is a Miller-Rabin witness for n
k=0 # Returns True or False, accordingly
temp=n-1
while (temp%2)==0: # First divide n-1 by 2 repeatedly, until we can't
k=k+1 # do so anymore.
temp=temp/2
q=(n-1)/(2^k) # Have n-1 = 2^k * q where q is odd
A=fastpower(a,q,n) # A=a^q in Z/n
if A==1: # If a^q=1, we DON'T have a M-R witness
return False
i=0
while i<=k-1:
if A==n-1: # Now test to see if A=-1 in Z/n
return False
A=(A**2)%n # If not, square A and continue the loop
i=i+1
# If we get here then none of a^q,a^{2q},...,a^{2^{k-1}q}
return True # equal -1 in Z/n, so we DO have a M-R witness
In [15]:
def count_mrw(n): # This counts the number of Miller-Rabin witnesses for n
i=1
count=0
while i<=n-1:
if miller_rabin_witness(n,i):
count=count+1
i=i+1
return count
In [16]:
count_mrw(9)
Out[16]:
In [17]:
def miller_rabin_test(n,k): # Randomly choose a number between 2 and n-1
i=1 # and check if it is a M-R witness for n.
while i<=k: # If yes, return it
a=randint(2,n-1) # If no, do it again...but stop after k tries.
if miller_rabin_witness(n,a):
print 'The number is composite: a Miller-Rabin witness is',a
return
i=i+1
print 'The number is probably prime'
return
In [18]:
def is_prime_mr(n,k): #Test if n is prime using Miller-Rabin, k times
if n==2:
return True
if n==1:
return False
j=1
while j<=k:
a=randint(2,n-1)
if miller_rabin_witness(n,a)==True:
return False
j=j+1
# print "Probably prime"
return True
In [19]:
def findaprime_mr(lower,upper,trials):
i=1
while i<=trials:
a=randint(lower,upper)
if is_prime_mr(a,10):
return a
i=i+1
print 'No prime found'
return 0
In [20]:
def pollard(N,bound): # Perform Pollard's p-1 algorithm
a=2 # to factor N, using "bound" as the number of steps in the loop.
j=2 # See page 139 of the text.
while j<=bound:
a=(a^j)%N
d=gcd(a-1,N)
if 1<d and d<N:
return d # Have found a nontrivial factor, so return it
j=j+1
print 'No factor found'
return
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [57]:
def plot2(f,g,a,b):
P=plot(f(x),x,a,b)+plot(g(x),x,a,b,color='red')
show(P)
return
In [ ]:
In [58]:
def plot3(f,g,h,a,b):
P=plot(f(x),x,a,b)+plot(g(x),x,a,b,color='red')+plot(h(x),x,a,b,color='green')
show(P)
return
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [230]:
prime_pi(10^6)
Out[230]:
In [235]:
plot(prime_pi(x),2,10000)+line([(0,0),(10^4,1229)],color='red')
Out[235]:
In [243]:
plot(prime_pi(x)/x,2,10^5)+plot(1/ln(x),2,10^5,color='red')
Out[243]:
In [258]:
plot(prime_pi(x),10^7,10^7+10^5)+plot(x/(ln(x)-1.08366),10^7,10^7+10^5,color='red')
Out[258]:
In [246]:
prime_pi(10^6)
Out[246]:
In [248]:
N(10^6/ln(10^6))
Out[248]:
In [249]:
def D(x):
a=prime_pi(x)
b=x/ln(x)
return N((a-b)/a)
In [253]:
D(10^12)
Out[253]:
In [259]:
def P(x):
return x/ln(x)
In [262]:
plot3(prime_pi,P,Li,10^5,10^5+10^4)
In [265]:
for k in [1,2,..,10]:
x=10^k
print N(prime_pi(x)),' ',N(Li(x)),' ',N(x/ln(x))
In [271]:
plot(Li(x)-prime_pi(x),10^5,10^8)+plot(1/150*x^0.5*ln(x),10^5,10^8,color='red')
Out[271]:
In [267]:
var('x')
Out[267]:
In [ ]: