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, 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 et d=ab=br vérifie d=ub+vr, alors

d=ub+v(abq)=va+(uqv)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 (impair) n'est pas premier, (mn) 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 []: