{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def shanks(p,g,h): # Solves g^x=h in F_p using the Shanks algorithm. \n",
" n=floor(p**0.5)+1\n",
" list_a=[]\n",
" i=0\n",
" temp=1\n",
" print 'Producing the list of powers of g...',\n",
" while i<=n-1:\n",
" list_a.append( [temp,i] )\n",
" temp=(temp*g)%p\n",
" i=i+1\n",
" print 'done.' # list_a now contains [g^0,0],[g^1,1],..,[g^{n-1},n-1]\n",
" # We are keeping track of both the power AND the exponent to make our life\n",
" # easier later\n",
" print 'Producing the list of h*g^{-jn}...',\n",
" list_b=[]\n",
" ginv=fastpower(g,p-2,p)\n",
" u=fastpower(ginv,n,p)\n",
" i=0\n",
" temp=h\n",
" while i<=n-1:\n",
" list_b.append( temp)\n",
" temp=(temp*u)%p\n",
" i=i+1\n",
" print 'done.' # list_b now contains hg^0,h*g^{-n},h*g^{-2n},...,h*g^{-(n-1)n}\n",
" print 'Starting collision algorithm for two lists of length', len(list_a)\n",
" \n",
" l=collision(list_a,list_b) # Do the collision algorithm. l[0]=i will be the index in the first list\n",
" # and l[1]=j will be the index in the second list\n",
" if l=='F':\n",
" print \"\\n There is no solution\"\n",
" return\n",
" return l[0]+n*l[1] # The solution is x=i+n*j.\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def fastpower(x,e,p): # Compute x^e in Z/p using the fast powering algorithm\n",
" temp=1\n",
" list=[]\n",
" if e==0:\n",
" return 1\n",
" while e>0:\n",
" if e%2==0:\n",
" list.append(2)\n",
" e=e/2\n",
" continue\n",
" list.append(1)\n",
" e=e-1\n",
" continue\n",
" a=list.pop()\n",
" temp=x\n",
" while len(list)>0:\n",
" a=list.pop()\n",
" if a==1:\n",
" temp=(temp*x)%p\n",
" continue\n",
" temp=(temp*temp)%p\n",
" return temp\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1, 2, 2, 2, 1]\n"
]
},
{
"data": {
"text/plain": [
"512"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fastpower(2,9,1024)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def collision(l,m): #Given two lists l=[(a0,b0),(a1,b1),...] and m=[x0,x1,x2,...]\n",
" #having the same length, see if some ai matches an xj.\n",
" #If no, return 'F' for 'Fail'.\n",
" i=0 #If yes, return the pair [i,j] giving the indices for the matching elements\n",
" print 'Sorting the first list...',\n",
"\n",
" l.sort(key=sortFirst) #First sort the list l so that a0<=a1<=a2<=...\n",
" print 'done.' #This can take a while for very long lists.\n",
" print 'Testing elements of the second list...', \n",
" i=0 #Now test elements of m one by one to see if they match one of the ai terms.\n",
" n=len(m)\n",
" while i<=n-1:\n",
" if i%100==0:\n",
" print '.',\n",
" x=m[i]\n",
" u=is_in(l,x) #The function is_in decides if x matches one of the ai\n",
" if u=='F':\n",
" i=i+1\n",
" continue\n",
" print 'done'\n",
" return [u,i]\n",
" print 'done.'\n",
" return 'F'\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def is_in(l,x): #Given l=[[a0,b0],[a1,b1],...] and x, decide if x equals one of the ai.\n",
" #Use the method of \"compare to middle term\" and divide the problem in half.\n",
" if len(l)==0: #If no, return 'F'. If yes, return the bi. \n",
" return 'F' #Assumes that a0<=a1<=a2<=... from the beginning.\n",
" first=0 \n",
" last=len(l)-1 #Set up pointers \"first\" and \"last\", originally at the start and end of the list l.\n",
" while first<last:\n",
" middle=first+((last-first)//2) #Go to roughly the middle term\n",
" u=l[middle][0] #and pull out the a-value\n",
" if u==x: #If it matches our x, return the corresponding b-value.\n",
" return l[middle][1]\n",
" if u<x: #If x is bigger, reset the pointers to look at the last half of the list\n",
" first=middle+1 #and repeat\n",
" continue\n",
" last=middle-1 #If x is smaller, reset the pointers to look at the first half of the list.\n",
" #If we come out of the loop, we never found our x.\n",
" if first>last: #Either we exhausted all elements, or we narrowed it down to exactly one.\n",
" return 'F'\n",
" u=l[first]\n",
" if u[0]==x: #Do one last check to see if we have the x.\n",
" return u[1]\n",
" return 'F'\n"
]
},
{
"cell_type": "code",
"execution_count": 245,
"metadata": {},
"outputs": [],
"source": [
"def findprime(lower,upper,tries):\n",
" i=1\n",
" while i<=tries:\n",
" a=randint(lower,upper)\n",
" if is_prime(a):\n",
" return a\n",
" i=i+1\n",
" print 'None found'\n",
" return 0"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def order(x,p):\n",
" i=1\n",
" temp=x\n",
" while i<=p:\n",
"# if i%100==0:\n",
"# print '\\r testing ',i,\n",
" if temp==1:\n",
" return i\n",
" temp=(temp*x)%p\n",
" i=i+1\n",
" print 'Error'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"35295289"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def find_prim_gen(p): # Find a primitive generator for F_p\n",
" x=2 # Warning: For big primes (9 or more digits) this can take a while.\n",
" while x<=p-1:\n",
" print 'Testing ',x,'...',\n",
" u=order(x,p)\n",
" print 'Order is ',u,'...',\n",
" if u==p-1:\n",
" print 'So found a generator'\n",
" return x\n",
" print 'So not a generator'\n",
" x=x+1\n",
" print 'Error'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Testing 2 ... Order is 18198982 ... So not a generator\n",
"Testing 3 ... Order is 72795928 ... So found a generator\n"
]
},
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"find_prim_gen(a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def dlp_naive(p,g,h): # Solve g^x=h in Z/p, given g and h. Uses the brute force algorithm. \n",
" i=0 # i keeps track of the exponent we are trying\n",
" temp=1 #Start with g^0\n",
" while i<=p:\n",
" if temp==h: #if we found it, return the exponent\n",
" return i\n",
" temp=(temp*g)%p # Otherwise multiply by one more copy of g and increase i.\n",
" i=i+1\n",
" print 'Error: no solution found'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"83562159"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dlp_naive(a,2,55)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dlp_naive(23,5,19)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def ph_estimate(n):\n",
" i=1\n",
" save=0\n",
" while i<=n:\n",
" b=findprime(10^20,10^21,35)\n",
" if b==0:\n",
" continue\n",
" l=factor(b-1)\n",
" count=0\n",
" for a in l:\n",
" count=count+a[1]*a[0]\n",
" save=save+(b+0.0)/count\n",
" i=i+1\n",
" return (save+0.0)/n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def miller_rabin_witness(n,a): # Determines if a is a Miller-Rabin witness for n\n",
" k=0 # Returns True or False, accordingly\n",
" temp=n-1\n",
" while (temp%2)==0: # First divide n-1 by 2 repeatedly, until we can't \n",
" k=k+1 # do so anymore.\n",
" temp=temp/2\n",
" q=(n-1)/(2^k) # Have n-1 = 2^k * q where q is odd\n",
" A=fastpower(a,q,n) # A=a^q in Z/n\n",
" if A==1: # If a^q=1, we DON'T have a M-R witness\n",
" return False\n",
" i=0\n",
" while i<=k-1:\n",
" if A==n-1: # Now test to see if A=-1 in Z/n\n",
" return False\n",
" A=(A**2)%n # If not, square A and continue the loop\n",
" i=i+1\n",
" # If we get here then none of a^q,a^{2q},...,a^{2^{k-1}q}\n",
" return True # equal -1 in Z/n, so we DO have a M-R witness"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"def count_mrw(n): # This counts the number of Miller-Rabin witnesses for n\n",
" i=1\n",
" count=0\n",
" while i<=n-1:\n",
" if miller_rabin_witness(n,i):\n",
" count=count+1\n",
" i=i+1\n",
" return count"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"count_mrw(9)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def miller_rabin_test(n,k): # Randomly choose a number between 2 and n-1\n",
" i=1 # and check if it is a M-R witness for n.\n",
" while i<=k: # If yes, return it\n",
" a=randint(2,n-1) # If no, do it again...but stop after k tries.\n",
" if miller_rabin_witness(n,a):\n",
" print 'The number is composite: a Miller-Rabin witness is',a\n",
" return\n",
" i=i+1\n",
" print 'The number is probably prime'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def is_prime_mr(n,k): #Test if n is prime using Miller-Rabin, k times\n",
" if n==2:\n",
" return True\n",
" if n==1: \n",
" return False\n",
" j=1\n",
" while j<=k:\n",
" a=randint(2,n-1)\n",
" if miller_rabin_witness(n,a)==True:\n",
" return False\n",
" j=j+1\n",
"# print \"Probably prime\"\n",
" return True"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def findaprime_mr(lower,upper,trials):\n",
" i=1\n",
" while i<=trials:\n",
" a=randint(lower,upper)\n",
" if is_prime_mr(a,10):\n",
" return a\n",
" i=i+1\n",
" print 'No prime found'\n",
" return 0"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def pollard(N,bound): # Perform Pollard's p-1 algorithm\n",
" a=2 # to factor N, using \"bound\" as the number of steps in the loop.\n",
" j=2 # See page 139 of the text. \n",
" while j<=bound:\n",
" a=(a^j)%N\n",
" d=gcd(a-1,N)\n",
" if 1<d and d<N:\n",
" return d # Have found a nontrivial factor, so return it\n",
" j=j+1\n",
" print 'No factor found'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"elliptic_p=13\n",
"elliptic_A=3\n",
"elliptic_B=8"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [],
"source": [
"def is_nonsing():\n",
" N=(4*elliptic_A^3+27*elliptic_B^2)%elliptic_p\n",
" return (N!=0)"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"is_nonsing()"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [],
"source": [
"# Given x and y, returns [gcd,A,B] where gcd(x,y)=Ax+By. Finds A and B with the Euclidean Algorithm.\n",
"\n",
"def euc_alg(x,y):\n",
" N=y\n",
" d=x\n",
" A1=0 \n",
" B1=1\n",
" A2=1\n",
" B2=0\n",
" while not(d==0):\n",
" quotient=N//d # Divide N by d, get the quotient\n",
" remainder=N%d # and the remainder\n",
" if remainder==0: # This means we found the gcd at previous step \n",
" l=[d,A2,B2] # so just return \n",
" return l\n",
" Atemp=A1-quotient*A2\n",
" Btemp=B1-quotient*B2\n",
" A1=A2\n",
" B1=B2\n",
" A2=Atemp\n",
" B2=Btemp\n",
" N=d\n",
" d=remainder\n",
" return 0"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[3, 2, -1]"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"euc_alg(15,27)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [],
"source": [
"# Finds the inverse of a in Z/q, if it exists, even when q is not prime.\n",
"# Returns 0 if the inverse does not exist\n",
"\n",
"def find_inverse(a,q):\n",
" l=euc_alg(a,q)\n",
" if l[0]==1:\n",
" return l[1]%q\n",
" if l[0]==-1:\n",
" return (0-l[1])%q\n",
" return 0\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 220,
"metadata": {},
"outputs": [],
"source": [
"# This code adds two points l1=[x1,y1] and l2=[x2,y2] in an elliptic curve.\n",
"# It assumes that the global variables elliptic_p, elliptic_A, and elliptic_B have been defined already\n",
"# Also, the point at infinity is represented by ['I','I'], and we use 'E' to represent an error that occurred.\n",
"#\n",
"def e_add(l1,l2):\n",
" x1=l1[0]\n",
" x2=l2[0]\n",
" y1=l1[1]\n",
" y2=l2[1]\n",
" if x1=='E': #If l1 is an \"error\" point, just return it again as an \"error.\"\n",
" return l1\n",
" if y1=='E':\n",
" return l2\n",
" if x1=='I': #If l1 is the point at infinity, return l2\n",
" return l2\n",
" if x2=='I': #If l2 is the point at infinite, return l1\n",
" return l1\n",
" if x1==x2 and ((y1+y2)%elliptic_p==0): #If l1 and l2 are mirror images of each other, return the point at infinity\n",
" return ['I','I']\n",
" if x1==x2 and y1==y2: #If the points are the same, use the tangent method\n",
" temp=find_inverse(2*y1,elliptic_p)\n",
" if temp==0:\n",
" return ['E',2*y1]\n",
" temp=temp*(3*x1^2+elliptic_A)\n",
" lambd=temp%elliptic_p # Note: \"lambda\" has a special meaning in Python, so can't use it\n",
" # for a variable. We are using \"lambd\" instead.\n",
" if x1==x2 and y1!=y2: # Should never encounter this case;\n",
" return ['E','E'] # if we do, return an \"error\"\n",
" if not(x1==x2):\n",
" temp=find_inverse(x2-x1,elliptic_p)\n",
" if temp==0:\n",
" return ['E',x2-x1]\n",
" lambd=(temp*(y2-y1))%elliptic_p\n",
" x3=lambd^2-x1-x2 # These are the formulas from Theorem 6.6 in the book.\n",
" x3=x3%elliptic_p\n",
" y3=lambd*(x1-x3)-y1\n",
" y3=y3%elliptic_p\n",
" l=[x3,y3]\n",
" return l"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3*67"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [],
"source": [
"# Finds all the points on the given elliptic curve. Assumes elliptic_p, elliptic_A, and elliptic_B have been defined\n",
"# as global variables.\n",
"\n",
"def e_points():\n",
" l=[]\n",
" squares=[]\n",
" x=0\n",
" while x<elliptic_p:\n",
" temp=(x^2)%elliptic_p\n",
" squares.append(temp)\n",
" x=x+1\n",
" x=0\n",
" while x<elliptic_p:\n",
" y2=(x^3+elliptic_A*x+elliptic_B)%elliptic_p\n",
" if squares.count(y2)>0:\n",
" i=squares.index(y2)\n",
" Q=[x,i]\n",
" l.append(Q)\n",
" if not(i==0):\n",
" Q=[x,elliptic_p-i]\n",
" l.append(Q) \n",
" x=x+1\n",
" l.insert(0,['I','I'])\n",
" return l"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [],
"source": [
"def e_power(P,n): #returns nP, using the fast powering algorithm\n",
" if n==1:\n",
" return P\n",
" if n%2==0:\n",
" temp=e_power(P,n/2)\n",
" temp=e_add(temp,temp)\n",
" return temp\n",
" temp=e_power(P,(n-1)/2)\n",
" temp=e_add(temp,temp)\n",
" temp=e_add(temp,P)\n",
" return temp"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# returns the multiplicative order of a in Z/N\n",
"# i.e., the smallest value of r such that a^r=1\n",
"#\n",
"def mult_order(a,N):\n",
" exp=1\n",
" temp=a\n",
" while exp<=N:\n",
" if temp==1:\n",
" return exp\n",
" temp=(temp*a)%N\n",
" exp=exp+1\n",
" print a,'is not a unit'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [],
"source": [
"# returns the additive order of P in the elliptic curve \n",
"# i.e., returns the smallest value of n such that nP=0\n",
"#\n",
"def e_order(P):\n",
" exp=1\n",
" temp=P\n",
" N=elliptic_p^2\n",
" while exp<=N:\n",
" if temp==['I','I']:\n",
" return exp\n",
" temp=e_add(temp,P)\n",
" exp=exp+1\n",
" print 'Error'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": 271,
"metadata": {},
"outputs": [],
"source": [
"# Performs the Lenstra algorithm for integer factorization using the point P\n",
"# (the integer to be factored should be set globally to be the value of elliptic_p )\n",
"def lenstra(P):\n",
" n=1\n",
" temp=P\n",
" while n<=50000:\n",
" temp=e_power(temp,n)\n",
" if temp[0]=='E':\n",
" print 'The Lenstra algorithm stopped at',n,'!'\n",
" print 'The number that does not have an inverse is ',temp[1]\n",
" return\n",
" n=n+1\n",
" print 'The algorithm got all the way to 50000! without finding a problem.'\n",
" return"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 8.7",
"language": "",
"name": "sagemath"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.15"
}
},
"nbformat": 4,
"nbformat_minor": 2
}