# M1 Cryptographie - TD4

## PGCD récursif

In [1]:
def gcd(a,b):
if b==0: return a
return gcd(b,a%b)

gcd(12,18)

Out[1]:
6


Notons que si $b>a$$b>a$, le premier appel récursif sera

gcd(a,b%a)

Il n'y a donc pas besoin de distinguer ce cas. L'algorithme étendu peut aussi se coder de manière récursive :

In [2]:
def xgcd(a,b):
if b==0:
return (a,1,0)
else:
d,u,v = xgcd(b,a % b)
return (d,v,u-(a/b)*v)
xgcd(17,26)

Out[2]:
(1, -3, 2)

In [3]:
xgcd(26,17)

Out[3]:
(1, 2, -3)


### Explication

Si $a=bq+r$$a=bq+r$ et $d=a\wedge b=b\wedge r$$d=a\wedge b = b\wedge r$ vérifie $d={u}^{\prime }b+{v}^{\prime }r$$d = u'b+v'r$, alors

$d={u}^{\prime }b+{v}^{\prime }\left(a-bq\right)={v}^{\prime }a+\left({u}^{\prime }-q{v}^{\prime }\right)b.$

In [4]:
def inversemod(a,n):
d,u,v = xgcd(a,n)
if d==1: return u%n
else: return None

In [26]:
p = 2**127-1
print p
print inversemod(666,p)

170141183460469231731687303715884105727
100909560761089108909934662113775107751


In [28]:
print inversemod(8,26)

None


In [5]:
def primes(n):
if n<2: return []
pp = [2]; mm = [m for m in range(3,n) if m%2]
while mm:
p = mm[0]
pp.append(p)
mm = [m for m in mm if m%p]
return pp

In [6]:
print primes(1000)

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]


In [7]:
ll = primes(100000)
len(ll)

Out[7]:
9592

In [8]:
pp = primes(1000000)
len(pp)

Out[8]:
78498

In [9]:
def powermod(a,m,n):
res = 1; x = a
while m:
if m%2: res = (res*x) % n
x = (x*x) % n
m /= 2
return res

In [10]:
powermod(3,2**107-2,2**107-1)

Out[10]:
1L

In [11]:
def trial_division(n):

if n%2 == 0: return 2
if n%3 == 0: return 3
if n%5 == 0: return 5
if n<10: return n
p = 3
while p*p <= n:
if n%p == 0: return p
p +=2
return n


In [12]:
trial_division(2**32+1)

Out[12]:
641

In [13]:
ll = [2,3]
p = 5
i = 2
while i<200000:
if trial_division(p) == p:
ll.append(p)
i+=1
p+=2

In [14]:
len(ll)

Out[14]:
200000

In [15]:
print ll[:100]

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541]


In [16]:
print ll[-10:]

[2749919, 2749921, 2749991, 2750021, 2750029, 2750053, 2750071, 2750123, 2750131, 2750159]


In [17]:
from random import randrange

def pseudo_prime(a,n):
return powermod(a,n-1,n)==1

def fermat(n,tests=4):
for i in range(tests):
a = randrange(2,n-1)
if not pseudo_prime(a,n): return False
return True

In [18]:
mm = [2,3]
p = 5
i = 2
while i<200000:
if fermat(p):
mm.append(p)
i+=1
p+=2

In [19]:
L = set(ll); M=set(mm)
failed = M.difference(L)

In [20]:
print failed

set([2113921, 530881, 294409, 1024651, 876961, 1152271, 399001, 1909001, 13981, 2704801, 2100901, 488881, 449065, 1314631, 188461, 172081, 251251, 334153, 252601, 1569457, 1729, 2465, 6601, 512461, 115921, 15841, 46657, 410041, 1857241, 2615497, 1615681, 162401, 340561, 2628073, 14701, 852841, 748657, 2508013, 1773289, 91001])


In [21]:
len(failed)

Out[21]:
40

In [22]:
def test_base(a,n):
m = n-1; k = 0
while not m%2:
k+=1; m/=2
b = powermod(a,m,n)
if b==1: return True
for i in range(k):
x = b
b = pow(b,2,n)
if b==1:
if x != n-1: return False
else: return True
return False

def miller_rabin(n, tests=4):
if n in [2,3]: return True
if not n%2: return False
for i in range(tests):
a = randrange(2,n-1)
if not test_base(a,n): return False
return True

In [23]:
nn = [2,3]
p = 5
i = 2
while i<200000:
if miller_rabin(p):
nn.append(p)
i+=1
p+=2

In [24]:
N = set(nn)
failed = L.difference(N)
print failed
print len(failed)

set([])
0


In [25]:
def solve_chinese(yy,mm):
assert len(yy) == len(mm), "Arguments must have same length."
k = len(yy); m = reduce(lambda x,y:x*y, mm); x=0
for i in range(k):
u = reduce(lambda x,y:x*y, [mm[j] for j in range(k) if j!=i])
v = inversemod(u,mm[i])
x = (x + yy[i]*v*u) % m
return x

solve_chinese([1,2,3,4], [11,13,17,19])

Out[25]:
41823


## Test de Solovay-Strassen

Lorsque $n$$n$ (impair) n'est pas premier, $\left(\frac{m}{n}\right)$$\left(\frac{m}{n}\right)$ est appelé symbole de Jacobi. Son calcul est basé sur la loi de réciprocité quadratique de Gauss.

In [76]:
def jacobi(m,n):
assert n%2, "n must be odd"
m = m%n
if m == 0: return 0
if m == 1: return 1

if m%2 == 0:
m /=2
if n%8 in [1,7]: return jacobi(m,n)
else: return -jacobi(m,n)
elif m%4 == 3 and n%4 == 3: return -jacobi(n,m)
else: return jacobi(n,m)

In [77]:
jacobi(1001,9907)

Out[77]:
-1

In [78]:
def test_sol_stra(a,n):
j = jacobi(a,n)
if j==0: return False
else: return (pow(a,(n-1)/2,n)-jacobi(a,n)) % n  == 0

In [79]:
n = 2**(2**5)+1
test_sol_stra(3,n)

Out[79]:
False

In [80]:
def solovay_strassen(n, tests=4):
for i in range(tests):
a = randrange(2,n-1)
if not test_sol_stra(a,n): return False
return True


In [86]:
aa = [2,3]
p = 5
i = 2
while i<200000:
if solovay_strassen(p):
aa.append(p)
i+=1
p+=2

In [87]:
A = set(aa)
failed = L.difference(A)
#print failed
print len(failed)

0


In []: