def gcd(a,b):
if b==0: return a
return gcd(b,a%b)
gcd(12,18)
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 :
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)
(1, -3, 2)
xgcd(26,17)
(1, 2, -3)
Si $a=bq+r$ et $d=a\wedge b = b\wedge r$ vérifie $d = u'b+v'r$, alors $$ d = u'b + v'(a-bq) = v'a +(u'-qv')b. $$
def inversemod(a,n):
d,u,v = xgcd(a,n)
if d==1: return u%n
else: return None
p = 2**127-1
print (p)
print (inversemod(666,p))
170141183460469231731687303715884105727 100909560761089108909934662113775107751
print (inversemod(8,26))
None
Pour trouver les nombres premiers entre 2 et n, on part d'une liste pp contenant seulement 2, et d'une liste mm
contenant tous les nombres entre 2 et n qui ne sont pas multiple de 2. mm[0] est donc 3, qui est forcément premier
puisque pas multiple des nombres premiers qui lui sont inférieurs. On fait donc passer 3 à la fin de pp et on élimine
de mm tous les multiples de 3. On recommence, 5 passe dans pp et les multiples de 5 sont éliminés. A chaque
étape, mm[0] est un nombre premier $p$, on l'ajoute à pp et on élimine ses multiples de mm, jusqu'à ce
que mm soit vide.
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[1:] if m%p]
return pp
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]
ll = primes(100000)
len(ll)
9592
A chaque étape, m%2 est un bit de l'exposant $m$. Le suivant est obtenu en remplaçant m par m//2 = m>>1.
9>>1
4
def powermod(a,m,n):
res = 1; x = a
while m:
if m & 1: res = (res*x) % n
x = (x*x) % n
m >>=1 # m //= 2
return res
powermod(3,2**107-2,2**107-1)
1
Un entier $<30$ qui n'est pas premier est forcément divisible par 2, 3, ou 5. On repart de 7.
Partant de cette observation, on pourrait améliorer l'algorithme pour faire moins de tests (wheel factorization).
Exercice : le faire ...
Solution : fichier ent3.py disponible la semaine prochaine.
def trial_division(n):
if n%2 == 0: return 2
p = 3
while p*p <= n:
if n%p == 0: return p
p +=2
return n
trial_division(2**32+1)
641
ll = [2,3]
p = 5
i = 2
while i<200000:
if trial_division(p) == p:
ll.append(p)
i+=1
p+=2
print (ll[-10:])
[2749919, 2749921, 2749991, 2750021, 2750029, 2750053, 2750071, 2750123, 2750131, 2750159]
from random import randrange
def is_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 is_pseudo_prime(a,n): return False
return True
mm = [2,3]
p = 5
i = 2
while i<200000:
if fermat(p):
mm.append(p)
i+=1
p+=2
L = set(ll); M=set(mm)
failed = M.difference(L)
print (failed)
{2113921, 1193221, 294409, 1909001, 950797, 1152271, 552721, 62745, 399001, 114589, 29341, 2465, 876961, 736291, 488881, 1569457, 252601, 1461241, 1615681, 46657, 2433601, 314821, 334153, 512461, 670033, 340561, 1857241, 15841, 35425, 75361, 162401, 63973, 525673, 2628073, 852841, 2508013, 91001}
len(failed)
37
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
nn = [2,3]
p = 5
i = 2
while i<200000:
if miller_rabin(p):
nn.append(p)
i+=1
p+=2
N = set(nn)
failed = L.difference(N)
print (failed)
print (len(failed))
set() 0
from functools import reduce
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([3,4,5],[17,11,6])
785
Lorsque $n$ (impair) n'est pas premier, $\left(\frac{m}{n}\right)$ est appelé symbole de Jacobi. Son calcul est basé sur la loi de réciprocité quadratique de Gauss.
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)
jacobi(1001,9907)
-1
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
n = 2**(2**5)+1
test_sol_stra(3,n)
False
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
aa = [2,3]
p = 5
i = 2
while i<200000:
if solovay_strassen(p):
aa.append(p)
i+=1
p+=2
A = set(aa)
failed = L.difference(A)
print (failed)
print (len(failed))
set() 0