In [48]:
#Enable LaTeX printing
from IPython.core.display import *
global MathJax
MathJax = True
def MDPL(string): display(Math(string)) if MathJax else display(Latex(string))
import IPython.core.display
IPython.core.display.HTML("""
div.input {
width: 105ex; /* about 80 chars + buffer */
}
div.text_cell {
width: 105ex /* instead of 100%, */
}
div.text_cell_render {
/*font-family: "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;*/
font-family: "Charis SIL", serif; /* Make non-code text serif. */
line-height: 145%; /* added for some line spacing of text. */
width: 105ex; /* instead of 'inherit' for shorter lines */
}
/* Set the size of the headers */
div.text_cell_render h1 {
font-size: 18pt;
}
div.text_cell_render h2 {
font-size: 14pt;
}
.CodeMirror {
font-family: Consolas, monospace;
}
""")
Out[48]:
In [49]:
#(Extended) Euclid Algorithm
#Code Block
#Derived from Algorithm by Robert W. Sebesta
def Euclid(a,b,ans=0,pt=0,retA=1):
"""Finds the greatest common divisor of a and b and x,y
such that a*x+b*y=gcd(a,b)
Three flags vary the output (Any number may be true):
ans: If true will print the equation a*x+b*y=gcd(a,b)
pt: If true will show the table method for finding the solution
retA: If true will return the list [a,x,b,y,gcd(a,b)]
"""
Temp=_Euc([1,0,a],[0,1,b],[[['A1','A2','A3'],['B1','B2','B3']],[[1,0,a],[0,1,b]]])
T=Temp[0]
if ans:
print '\n({0})({1})+({2})({3})={4}'.format(a,T[0],b,T[1],T[2])
if pt:
for x in Temp[1]: print '{0}\t{1}\t{2}\t{3}\t{4}\t{5}'.format(x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2])
print
if retA:
return [a,T[0],b,T[1],T[2]]
def _Euc(a,b,s):
"""Main recursive loop for Extended Euclidean Algorithm"""
if b[2]==0: return [b[0],b[1],a[2]],s
if b[2]==1: return b,s
else:
Q=a[2]/b[2]
t=[a[n]-Q*b[n] for n in range(3)]
a=b
b=t
return _Euc(a,b,s+[[a,b]])
def MI(a,b):
"""Returns the multiplicative inverse of a mod b"""
temp = Euclid(a,b,ans=0,pt=0,retA=1)[1]
if temp<0:temp+=b
return temp
In [50]:
#Chinese Remainder Theorum
#Code Block
pprint=MDPL
def CRT(A,ptable=0):
"""A is a list of lists such that [[a1,a2],[b1,b2],...] is equivalent to
x = a1 mod a2, x = b1 mod b2,... where a,b,... are integers
The flag ptable if true will output the work to find the solution"""
M=reduce(lambda x,y:x*y,[i[1] for i in A])
if ptable: pprint( 'M={0}'.format(M))
count=1;Ans=[];outS='';sumS=0
for i in A:
Ans=Ans+[[count,i[0],i[1],M/i[1],Euclid(M/i[1],i[1],retA=1)[1]%i[1]]]
#print 'a({0})={1} m({0})={2} M/M({0})={3} c({0})={4}'.format(Ans[count-1][0],Ans[count-1][1],Ans[count-1][2],Ans[count-1][3],Ans[count-1][4])
if ptable: pprint('a_{{{0}}}={1},\ m_{{{0}}}={2},\ \\frac{{M}}{{M_{{{0}}}}}={3},\ c_{{{0}}}={4}'.format(Ans[count-1][0],Ans[count-1][1],Ans[count-1][2],Ans[count-1][3],Ans[count-1][4]))
outS=outS+'({0}*{1}*{2})+'.format(Ans[count-1][3],Ans[count-1][4],Ans[count-1][1])
sumS=sumS+(Ans[count-1][3]*Ans[count-1][4]*Ans[count-1][1])
count=count+1
if ptable:
pprint( outS[:-1] + ' = {0}'.format(sumS) )
pprint( '= {0}\ mod\ {1}'.format(sumS%M,M) )
return sumS%M
#ChRem([[5,7],[7,11],[3,13]]);print
#CRT([[5,8],[3,5]],ptable=1);print
#ChRem([[7,11],[9,7],[2,3],[1,2]])
In [51]:
#Basic Encryption/Decryption Algorithms
#Code Block
from binascii import * #Hex conversion
from collections import deque #Queue
from collections import defaultdict as dd #Dynamic Dictionary
from string import Template #Output formatting
from heapq import heappush, heappop #Minheap
import re
def type2(a):
try:
b=type(a)
except:
b='NULL'
return b
def remDups(seq):
seen = set()
seen_add = seen.add
return reduce(lambda x,y:x+y,[ x for x in seq if x not in seen and not seen_add(x)])
class Cipher:
def __init__(self):
pass
def enc(self,s):
print "No Encryption Algorithm Added"
def dec(self,s):
print "No Decryption Algorithm Added"
class MAS(Cipher):
"""Monoalphabetic (Affine) Substitution Cipher"""
def __init__(self,a,b):
self.a=a
self.b=b
self.an=Euclid(a,95,retA=1)[1]
def enc(self,s):
return reduce(lambda x,y:x+y,[chr(((ord(p)-32)*self.a+self.b)%95+32) for p in s])#if (ord(p)>=97 and ord(p)<=122) else p for p in s])
def dec(self,s):
return reduce(lambda x,y:x+y,[chr((((ord(p)-32)-self.b)*self.an)%95+32) for p in s])#if (ord(p)>=97 and ord(p)<=122) else p for p in s])
class PF(Cipher):
"""Playfair Cipher"""
def __init__(self,s):
Temp=remDups(s.replace('q','')+'abcdefghijklmnoprstuvwxyz')
Temp2=[]
for i in range(5):
Temp2=Temp2+[Temp[i*5:(i+1)*5]]
self.Table=Temp2
def printTable(self):
for i in self.Table:
print '{0}'.format(str(i))
In [52]:
#Primality Tests
##Code Block
from random import randint as ri
ppt=MDPL
class PTest:
def __init__(self):
pass
def TryDiv(self,a):
return [(i,a/i) for i in range(2,int(a**.5)) if a%i==0]
def Fermat(self,a):
pass
def Carmic(self,a):
pass
def MilRab(n,rep=5,showwork=0):
if n<=3: return 1
if n%2==0:
if showwork: ppt("{0}\ is\ even".format(n))
return 0
n1=q=(n-1)
k=0
while q%2==0:
q=(q)/2
k+=1
c=pow(2,q,n)
d,s=q,k #Other explanations use this notation
if showwork: ppt("{2}-1 = 2^{0} {1}".format(k,q,n)) #Show n-1 factored
if q==1: return 1
def _MRHelp(a):
if pow(a,d,n)==1: return 0
for i in range(s): #Witness Loop
temp2=pow(a,2**i*d,n)
if temp2==n-1: return 0
return 1
for temp in range(rep):
a=ri(2,n-2)
if _MRHelp(a): return 0
return 1 #Since nothing says otherwise, return probably prime
In [53]:
#Factorization
#Code Block
class Factorizer:
def __init__(self):
pass
def TryDiv(self,a):
return [(i,a/i) for i in range(2,int(a**.5)) if a%i==0]
def Fermat(self,a):
ret=[]
for i in range(1,100):
temp=(a+i**2)**.5
if temp==int(temp):
temp=int(temp)
ret+=[[temp+i,temp-i]]
return ret
def Pollard(self,a):
pass
def QuadSieve(self,a):
pass
In [54]:
#Pohlig Hellman Discrete Log
#Code Block
#Source: http://www-math.ucdenver.edu/~wcherowi/courses/m5410/phexam.html
from sympy import factorint as fc
ppt=MDPL #LaTeX-style print function
def DLog(g,p,y,ShowWork=0):
"""Prompt for a generator g, prime p, and y=g^x value.
Essentially finds x
Compute x using the Pohlig Hellman algorithm."""
pf=fc(p-1) #Factor p-1
if ShowWork:
ppt("Pollig-Hellman\ Discrete\ Log\ Example")
ppt("Given\ p={0},g={1},\ and\ y={2}\ such\ that\ {1}^x\equiv {2}\mod {0} \equiv {3}\mod {0},\ find\ x:".format(p,g,y,y%p))
ppt("$p-1={0}-1={1}$".format(p,"".join([str(i[0])+"^"+str(i[1])+" * " for i in pf.iteritems()]))[:-3])
ppt("Now\ we\ need\ to\ find\ a\ congruency\ for\ each\ factor:")
x=[]
for i in pf.iteritems(): #For each factor
q=i[0] #Factor
r=i[1] #Power
if ShowWork:
print
ppt("Find\ x_{{{0}}} \mod {1}\ [x_{{{0}}}={2}]:".format(q,q**r,"".join(["c_{{{0}}}*({{{1}}})+".format(ii,q**ii) for ii in range(r)])[:-1]))
yy=y #Initially set base to y for each factor
NewCLast=dict((pow(g,ii*(p-1)/q,p),ii%q) for ii in range(1,q+1)) #Dictionary of possible values for c_? to hit
if ShowWork: ppt("y^{{ k \\frac{{p-1}}{{q^{{c_0}}}} }}\ must\ be\ one\ of\ these:\ {0}".format("".join(["{2}^{{{0}*\\frac{{p-1}}{{{3}}}}} \equiv {1},\ ".format(NewCLast[ii],ii,g,q) for ii in NewCLast.keys()])[:-3] ) )
xTemp=0 #Using cumulative sum for finding each c_?
for k in range(r): #For each c_?
for j in range(1,q+1): #Try values until temp is in the NewCLast dictionary
j=j%q #Since it would normally start with 0, making 0 last by wrapping with mod
trypow=j*(p-1)/(q**(k+1)) #Generate power to raise yy by
temp=pow(yy,trypow,p)
if temp in NewCLast.keys():
j=NewCLast[temp] #Find the power of yy that makes temp
xTemp+=(q**k)*j #Add the cumulating sum for x
if ShowWork: ppt("$ y^{{ k \\frac{{p-1}}{{q^{{{9}}}}} }} \equiv {0}^{{ {7} \\frac{{{1}}}{{{2}}} }} \equiv {0}^{{{3}}} \equiv {4} \mod {6},\ so\ c_{{{8}}} = {7} $".format(yy,p-1,q**(k+1),trypow,temp,temp-p,p,j,k,k+1))
#cLast=[j]
if k<r-1: #Skip division if unneccessary
if ShowWork: ppt("We found\ c_{{{1}}},\ so\ divide\ {0}\ by\ {3}^{{{2}*c_{{{1}}} }}\ ({3}^{{ -{5}*{4} }} \equiv {6})".format(yy,k,q**k,g,j,q**k, (MI(g,p)**(q**j))%p ) )
yy=( yy*MI(g,p)**( q**(k)*(j) ) )%p
break
if ShowWork: ppt("So,\ x_{{{0}}}={1} \mod {2}".format(q,xTemp,q**r))
x+=[[xTemp,q**r]] #Add equivalence to list for CRT
if ShowWork:
print
ppt("Now\ use\ the\ Chinese\ Remainder\ Theorem:")
return CRT(x,ptable=ShowWork)
In [55]:
#P3
#Code Block
from sympy import factorint as fc
from sympy import randprime as rp
from random import randint as ri
from binascii import hexlify
from binascii import unhexlify
def isProbablyPrime(a,b=5):
"""Uses Miller-Rabin to return if a is prime using b repititions"""
return MilRab(a,rep=b,showwork=0)
def getProbablePrime(a,b,speedup=0):
"""Return a prime of bit length a using b repititions of Miller-Rabin"""
r=4
if speedup: r=rp(2**(a-1),2**(a))
else:
while not(isProbablyPrime(r,b)):
r=ri(2**(a-2),2**(a-1))*2-1#int(reduce(lambda x,y:x+y,[str(ri(0,1)) for i in range(a)]),2)
return r
def Gcd(a,b):
"""Returns gcd of integers a and b"""
return Euclid(a,b,ans=0,pt=0,retA=1)[4]
def MultiplicativeInverse(a,b):
"""Returns the inverse of a mod b"""
temp = Euclid(a,b,ans=0,pt=0,retA=1)[1]
if temp<0:temp+=b
return temp
def GenerateRSAKeyPair():
"""generate large (probable) primes p and q,
and appropriate e and d and print them to the screen.
It would be helpful to allow interaction with a file
(put the public key information into one file and
private key information into another file for input
into later functions)."""
p=getProbablePrime(1024,5,1)
q=getProbablePrime(1024,5,1)
n=p*q
e=ri(3,2**10+3)
phiN=(p-1)*(q-1)
while Gcd(e,phiN)!=1: e=ri(3,2**10+3)
d=MultiplicativeInverse(e,phiN)
return (e,n),(d,n)
def RSAEncrypt(p,e):
e,n=e
m=int(hexlify(p),16)
m=hex(pow(m,e,n))[2:-1]
if len(m)%2==1: m='0'+m
m=unhexlify(m)
return m
def RSADecrypt(c,d):
return RSAEncrypt(c,d)
def PohligHellmanDiscreteLogarithm(g,p,y):
"""Prompt for a generator g, prime p, and y=g^x value.
Essentially finds x
Compute x using the Pohlig Hellman algorithm."""
return DLog(g,p,y)
In [83]:
TEST=MAS(59,34)
pt="""Gen. 1:27 So God created man in his own image, in the image of God he created him; male and female he created them."""
ct=TEST.enc(pt)
print ct
print TEST.dec(ct)#.replace("[","")
d=Euclid(17,12195156688,retA=1)[1]
n=12195377863
ct=RSAEncrypt("""Hey Alice, isn't RSA cool!? --Bill""",[17,12195377863])
pt=RSADecrypt(ct,[d,12195377863])
print "\nRSA:"
print ct
print pt
ct=[4048208982,7973981835,2232860160,1228506572,2392521799,2353518709,4207703255,3703820767,250037531,6326907769,3298628742]
for i in ct:
print pow(i,d,n)
In [56]:
from math import log
def texttoblocks(text_,blocklen_) :
Map[Sum[#[[i]] 95^(i-1), {i, Length[#]}] & ,Partition[Join[texttonums[text],Table[0,
{j,Mod[-Length[texttonums[text]],blocklen]}]],blocklen]]
def blocklength(M): return int(log(95,M))
def RSAencode(string_,e_,M_) :
return pow(texttoblocks(string_,blocklength(M_)),e_,M_)