San Vu Ngoc - IRMAR
%matplotlib inline
#notebook
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.cm as cm
import matplotlib.image as img
from scipy import signal
import cv2
plt.rcParams['figure.figsize'] = (10.0, 6.0)
Ératosthène est un savant Grec qui a vécu aux alentours de 200 ans avant J.-C.
def crible (n): # crible d'Ératosthène (sans optimisation)
P = np.full(n, True, dtype=bool)
P[0]=False # 0 n'est pas premier
P[1]=False # 1 n'est pas premier
prem = 2
while (prem < n):
if P[prem]:
print (str(prem) + ' est premier !')
for j in range(2*prem, n, prem):
if P[j]:
print (" on supprime:" + str(j))
P[j] = False
prem += 1
return P
P=crible(100)
P
def affiche_premier (T):
for i in range(np.size(T)):
if T[i]:
print i,
affiche_premier(P)
Image10 = np.zeros((1,100))
Image10[0,:] = P
plt.imshow(Image10, interpolation='none', cmap=cm.gray)
# Les nombres premiers apparaissent en blanc (True=1=blanc) et les autres en noir (False=0=noir)
Image10 = np.reshape(P,(10,10))
plt.imshow(Image10, interpolation='none', cmap=cm.gray)
for i in range(10):
for j in range(10):
plt.annotate(str(10*j+i), xy=(i-0.2, j+0.1), color='r')
(version un peu optimisée)
def crible2 (n): # crible d'Erathosthène (avec borne optimisée)
P = np.full(n, True, dtype=bool)
P[0] = False # 0 n'est pas premier
P[1] = False # 1 n'est pas premier
prem = 2
while (prem < n):
if P[prem]:
for j in range(prem*prem, n, prem):
if P[j]:
P[j] = False
prem += 1
return P
n=100
T=crible2(n*n)
affiche_premier(T)
Il y a une infinité de nombres premiers...
Mais est-ce qu'il y a autant de nombres premiers "autour de 1000" et "autour de 10^45" ?
Image100 = np.reshape(T,(n,n))
plt.imshow(Image100, interpolation='none', cmap=cm.gray)
T1000 = crible2(1000000)
Image1000 = np.reshape(T1000,(1000,1000))
C'est-à-dire le nombre de premiers dans une fenêtre donnée autour de chaque (x,y). On choisit une fenêtre carrée 3x3 ou 11x11.
Fenetre3=np.ones((3,3))
Fenetre11=np.ones((11,11))
Voisins10=signal.convolve2d(Image10.astype(int),Fenetre3, mode='same')
plt.imshow(Voisins10, interpolation='none', cmap=cm.gray)
for i in range(10):
for j in range(10):
plt.annotate(str(10*j+i), xy=(i-0.2, j+0.1), color='r')
plt.colorbar()
Voisins100=signal.convolve2d(Image100.astype(int),Fenetre11, mode='same')
plt.imshow(Voisins100, interpolation='none', cmap=cm.gray)
plt.colorbar()
plt.imshow(Voisins100, interpolation='none', cmap=cm.spectral)
plt.colorbar()
# un peu plus rapide:
#Voisins2=cv2.filter2D(Image.astype(int),-1,Fenetre)
Voisins1000=signal.convolve2d(Image1000.astype(int),Fenetre11, mode='same')
plt.imshow(Voisins1000, interpolation='lanczos', cmap=cm.spectral)
plt.colorbar()
#img.imsave('Voisins1000x1000.png', Voisins1000, cmap=cm.spectral)
def comptage (P):
n = np.size(P)
C = np.zeros(n)
x = 0 # variable de comptage
for i in range (n):
if P[i]:
x += 1
C[i] = x # nombre de premiers inférieurs ou égaux à i
return C
comptage(crible2(100))
n = 1000
x = np.arange(0, n, 1)
P = crible2(n)
y = comptage(P)
Image_compt = plt.plot(x, y)
#plt.axis([0, 6, 0, 20])
Est-ce que c'est le cas ??
x2 = np.arange(2, n, 1)
y2 = x2/np.log(x2)
fig, ax = plt.subplots()
ax.plot(x,y, color='b')
ax.plot(x2,y2, color='r')
plt.show()
C'est presque le cas, mais pas tout-à-fait. On montre que $\pi(x)$ se comporte en première approximation comme $x/\ln x$. (On "perd" donc un facteur $1/\ln(x)$)
fig, ax = plt.subplots()
ax.plot(x2,y2/y[2:], color='g')
ax.plot(x2,np.ones(np.size(x2)))
ax.set_ylim([0,1.5])
plt.show()
La convergence vers 1 est très lente !
Pour $x=10^{23}$, on est seulement à 1 moins 2%.
On connaît des équivalents qui convergent plus vite. Par exemple:
$${\mathrm{Li}}(x) = \int_2^x \frac{dt}{\ln t}.$$
x3 = np.arange(2,1000,0.1)
y3 = 1/np.log(x3)
plt.plot(x3,y3)
np.trapz(y3,x3)
def f(x):
X = np.arange(2,x,0.1)
y = 1/np.log(X)
return(np.trapz(y,X))
f(3)
Y3 = [f(i) for i in x3]
fig, ax = plt.subplots()
ax.plot(x3,Y3, color='g')
ax.plot(x2,y2, color='r')
ax.plot(x,y, color='b')
plt.show()