Python v2.7

[modification le 20/3/2012: refonte complète de la page.]

En informatique, générer du hasard, en l’occurrence des nombres au hasard, ce n'est pas facile. La plupart des algorithmes connus ont des défauts.

Un bon générateur doit avoir au moins ces qualités:

• que tous les nombres aient la même chance d'être tirés,
• que chaque tirage soit indépendant des précédents (donc, même en connaissant une partie de la série, on est incapable de prévoir la suite)
• et que la période soit assez grande, c'est à dire que la séquence de départ ne se retrouve pas à l'identique trop rapidement.

On va coder ici l'un des “bons” algorithmes: le “Blum Blum Shub”. Les explications sont chez wikipedia: http://fr.wikipedia.org/wiki/Blum_Blum_Shub.

Même s'il est considéré comme lent, cet algorithme est suffisamment bon pour être reconnu par les cryptographes.

Pour comprendre la base de l'algorithme, voilà un petit code qui est la traduction fidèle de l'explication de wikipedia:

gr = 4575 # => graine
np1 = 24375763 # => 1er nombre premier
np2 = 28972763 # => 2e nombre premier
pnp = np1*np2 # => produit des 2 nombres premiers
 
def bitalea():
    """retourne un bit (0 ou 1) tiré au hasard """
    global gr, pnp
    gr = (gr * gr) % pnp
    return gr & 1 # on retourne seulement le bit le moins significatif
 
ch = ""
for i in range(80):
    ch += str(bitalea())
print ch

Ce qui affiche la liste des 80 bits demandés:

11001111110001111100001000010110001000010010111001010011000010101001011101110100

Dans cette chaine, tous les bits de gauche à droite sont tirés au hasard. A noter que, de ce fait, ils sont aussi au hasard de droite à gauche!

On voit que l'algorithme dépend de 3 nombres: gr = la “graine”, np1 = 1er nombre premier et np2 = 2e nombre premier. Mais c'est le produit de ces 2 nombres premiers qu'on utilise dans la fonction.

Chacun des 2 nombres premiers doivent en plus être “congru à 3 modulo 4”, c'est à dire qu'on doit avoir: (np-3)%4=0

La graine peut être quelconque, mais certains articles (pas tous), sur ce générateur indiquent que l'on améliore la qualité de l'algorithme, et en particulier la longueur de la période de la série, en choisissant une graine telle que: PGCD(gr, pnp)=1. Comme c'est facile à faire, c'est ce que nous ferons ici.

Le code plus loin dans cette page donne des solutions pour calculer des nombres qui vont bien: nombres premiers et graine.

Utilisé en cryptographie, donc pour empêcher quelqu'un de non autorisé de déchiffrer un message crypté, il faut utiliser des nombres premiers très grands, de sorte qu'il soit quasi impossible, dans l'état actuel de la technique, de retrouver les 2 nombres premiers à partir de leur produit (factorisation des grands nombres). On verra plus loin comment construire des nombres premiers aussi grands que ceux utilisés dans le cryptage RSA.

Point important: si on recommence à lancer le même algorithme avec les mêmes trois nombres, on obtient la même suite! Cela peut être intéressant pour certaines applications (la cryptographie justement), mais pas toujours.

Si on veut du hasard depuis le lancement de la fonction, on peut calculer la graine avec des éléments matériels qui évoluent en permanence (ex: le nombre d'heures depuis le 1/1/1970, nombre d'octets occupés sur le disque dur système, etc…).

Par exemple en utilisant le nombre de seconde depuis l'epoch (=1/1/1970):

import time
graine = int(time.time())
print graine
1332223339

Pour accélérer les calculs, on peut extraire plusieurs bits à chaque fois, mais dans la limite de log(log(pnp)) afin de maintenir la qualité du hasard. Voyons à quoi ça correspond en gros:

pnp à    2  chiffres => 1 bit à extraire
pnp à    4  chiffres => 2 bits à extraire
pnp à    9  chiffres => 3 bits à extraire
pnp à   24  chiffres => 4 bits à extraire
pnp à   65  chiffres => 5 bits à extraire
pnp à  176  chiffres => 6 bits à extraire
pnp à  477  chiffres => 7 bits à extraire
pnp à 1295  chiffres => 8 bits à extraire

Nous utiliserons ces valeurs en tant que simplification pour éviter le calcul de log(log(pnp))

Comme on ne veut pas seulement des bits au hasard mais des nombres entiers ou flottants, et que pour être plus rapides, il va falloir exploiter la possibilité d'utiliser plusieurs bits à chaque calcul, on va créer plusieurs fonctions construites comme des itérateurs.

Reprenons d'abord la fonction précédente qui fournit un bit à la fois:

def genbit_bbs(graine, pnp):
    """itérateur retourne un bit au hasard (entier = 0 ou 1)"""
    while True:
        graine = pow(graine, 2, pnp)
        yield int(graine & 1) # retourne le bit de poids faible

On peut utiliser cet itérateur de 2 façons différentes, par exemple pour obtenir une liste de 50 bits au hasard:

gr = 4575 # => graine
np1 = 24375763 # => 1er nombre premier
np2 = 28972763 # => 2e nombre premier
pnp = np1*np2 # => produit des 2 nombres premiers
 
# 1ère utilisation possible de l'itérateur
gen = genbit_bbs(gr, pnp)
n = 50
L1 = []
for i in xrange(0, n):
    L1.append(gen.next())
print L1
 
# 2ème utilisation possible de l'itérateur
n = 50
L2 = []
for i, r in enumerate(genbit_bbs(gr, pnp)):
    L2.append(r)
    if i>=n:
        break
print L2

On va ici obtenir des nombres entiers au hasard de 0 à 7 (bornes comprises), qui pourront nous servir pour faire des rotations d'octets (cryptographie).

Et comme on va extraire 3 bits d'un coup, il nous faudra un produit pnp aussi grand ou plus grand que 9 chiffres (cf. condition du log(log(pnp)) ).

def gen3bits_bbs(graine, pnp):
    """itérateur retourne un nombre entier au hasard de 0 à 7 (bornes comprises)
       version rapide qui nécessite:  pnp à au moins 9 chiffres
    """
    while True:
        graine = pow(graine, 2, pnp)
        yield int(graine & 7) # retourne les 3 bits de poids faible 

Exemple d'utilisation: on va calculer une liste de 50 nombres de 0 à 7 au hasard:

gr = 4575 # => graine
np1 = 24375763 # => 1er nombre premier
np2 = 28972763 # => 2e nombre premier
pnp = np1*np2 # => produit des 2 nombres premiers
 
gen = gen3bits_bbs(gr, pnp)
n = 50
L = []
for i in xrange(0, n):
    L.append(gen.next())
print L
[1, 1, 6, 2, 7, 1, 1, 3, 7, 1, 0, 2, 2, 1, 1, 7, 1, 7, 0, 2, 6, 4, 1, 2, 4, 2, 6, 5, 6, 1, 3, 0, 4, 4, 5, 6, 6, 4, 4, 3, 4, 6, 7, 4, 3, 3, 5, 0, 2, 1]

Comme on a souvent besoin d'un générateur d'octets, on va en fabriquer un qui soit en même temps simple et rapide.

Sa seule contrainte est qu'il nécessitera que pnp, le produit des 2 nombres premiers, ait au moins 24 chiffres décimaux: cela permettra de créer l'octet avec seulement 2 calculs puis qu’alors on pourra extraire 4 bits en même temps à chacun des calculs, tout en respectant le log(log(pnp)).

def gen8bits_bbs(graine, pnp):
    """itérateur retourne un octet au hasard selon le générateur Blum Blum Shub
       version rapide qui nécessite:  pnp à au moins 24 chiffres décimaux
    """
    while True:
        graine = pow(graine, 2, pnp)
        r = graine & 15 # récup des 4 bits inférieurs
        graine = pow(graine, 2, pnp)
        r |= (graine & 15)<<4 # récup et ajout des 4 bits supérieurs
        yield int(r) # retourne l'octet 

Ce code est assez rapide: avec 2 nombres premiers de 512 bits, chaque nouvel octet ne demande que 6/100000 secondes!

On va créer ici un code assez général dont on se servira plus loin: générer des nombres entiers composés d'un nombre donné de bits. L’intérêt, c'est qu'on pourra toujours extraire le maximum de bits à chaque fois, connaissant le produit pnp:

def genalea_bbs(graine, pnp, nbits=8):
    """itérateur retourne un nombre entier de nbits (défaut=8) au hasard 
       version rapide qui extrait le nb de bits maxi selon la taille de pnp
    """
    # calcul du nombre de bits qu'on peut extraire d'un seul coup
    E = ((4,2), (9,3), (24,4), (65,5), (176,6), (477,7), (1295,8))
    v0, ok = 1, False
    for k, v in E:
        if pnp < 10**k:
            nbext, ok = v0, True
            break
        else:
            v0 = v
    if not ok:
        nbext = 8 # pnp pas trouvé dans la liste: 8 bits maxi extractibles        
    # calcul des masques
    masque1 = (1<<nbext)-1 # exemple: masque=7 pour ex=3 bits
    q, r = divmod(nbits, nbext) # on extrait q fois nbext bits + 1 fois r bits
    masque2 =  (1<<r)-1 # exemple: masque=3 pour r=2 bits
    # itérateur
    while True:
        n = 0
        for i in xrange(q):
            graine = pow(graine, 2, pnp)
            n = n | ((graine & masque1)<<(i*nbext))
        graine = pow(graine, 2, pnp)
        n = n | ((graine & masque2)<<(q*nbext))
        yield int(n)

Exemple d'utilisation: obtenir une liste de 10 nombres au hasard composés de 32 bits.

gr = 4575 # => graine
np1 = 24375763 # => 1er nombre premier
np2 = 28972763 # => 2e nombre premier
pnp = np1*np2 # => produit des 2 nombres premiers
 
gen = genalea_bbs(gr, pnp, 32)
n = 10
L = []
for i in xrange(0, n):
    L.append(gen.next())
print L
[258274697, 840757842, 54191377, 1049513388, 4142440795L, 3127086357L, 3178412812L, 660394131, 1557182431, 1857775838]

On va utiliser la fonction précédente pour générer des nombres entiers quelconques entre 2 valeurs n1 et n2:

def genentier_bbs(graine, pnp, n1, n2):
    """itérateur retourne des nb entiers au hasard entre n1 et n2 (bornes comp.) """    
    dn = n2-n1
    nbits = len(bin(dn))-2 # nb de bits de dn
    gen = genalea_bbs(graine, pnp, nbits)
    while True:
        n = gen.next()
        if n <= dn:
            yield n1+n

Exemple d'utilisation: liste de 20 nombres au hasard compris entre 1000 et 9999 (bornes comprises).

gr = 4575 # => graine
np1 = 24375763 # => 1er nombre premier
np2 = 28972763 # => 2e nombre premier
pnp = np1*np2 # => produit des 2 nombres premiers
 
gen = genentier_bbs(gr, pnp, 1000, 9999)
n = 10
L = []
for i in xrange(0, n):
    L.append(gen.next())
print L
[8881, 5752, 9655, 2126, 8538, 7403, 6163, 4349, 7242, 6931, 1735, 8086, 8846, 1087, 3717, 1373, 9006, 5937, 6413, 7424]

On va générer ici des nombres décimaux au hasard entre 0.0 et 1.0, et dont le nombre de chiffres après la virgule correspond à la précision des flottants Python (donc des “double” du C).

Le principe est simple: on va générer des nombres entiers compris entre 0 et 10**17, qu'on va diviser à la fin par 10**17.

def genflottant_bbs(graine, pnp, dn=100000000000000000):
    """itérateur retourne des nb flottants au hasard entre 0 et 1 (bornes comp.)"""    
    nbits =  len(bin(dn))-2 # nb de bits de dn (=précision demandée pour le nb)
    gen = genalea_bbs(graine, pnp, nbits)
    while True:
        n = gen.next()
        if n<=dn:
            yield n/dn

Exemple d'utilisation: générer une liste de 10 nombres flottants au hasard.

    gr = 4575 # => graine
    np1 = 24375763 # => 1er nombre premier
    np2 = 28972763 # => 2e nombre premier
    pnp = np1*np2 # => produit des 2 nombres premiers
 
    gen = genflottant_bbs(gr, pnp)
    n = 10
    L = []
    for i in xrange(0, n):
        L.append(gen.next())
    print L
[0.1629546695655156, 0.8298921698743204, 0.43617398308304073, 0.9818278003550986, 0.13144280334580596, 0.4636071208964701, 0.9566706890174403, 0.00977197244581526, 0.34704090744976773, 0.8498426522964183]

Les nombres premiers qui “conviennent” sont:

• des nombres premiers. Ils seront assez grands, surtout en cryptographie, puisque la sécurité dépendra directement de la difficulté à factoriser
• qui sont en plus “congrus à 3 modulo 4”.

Dans les essais que j'ai fait, j'ai trouvé que des nombres premiers de 24 bits permettaient déjà de trouver des séries aléatoires d'un million de digits (0 à 9) sans rencontrer de période à l'intérieur.

Pour trouver de tels nombres premiers de 24 bits:

• on tire au hasard avec randint (du module random) un nombre compris entre 8388608 (=int(“100000000000000000000000”, 2)) et 16777215 (=int(“111111111111111111111111”, 2)) bornes comprises
• on le teste avec nos 2 conditions: être congru à 3 modulo 4 et être premier
• et on boucle jusqu'à avoir trouvé un nombre qui passe le test.

Pour le test dit “de primalité” (dire si un nombre est premier), j'ai mis au point une fonction qui fait appel aux bons algorithmes selon la taille du nombre:

• <1024 ⇒ recherche dans une liste
• < 100000000 ⇒ méthode des divisions
• au dessus: méthode probabiliste de Miller-Rabin (étudié sur une autre page de ce site).

Voilà un code qui fait ça:

#############################################################################
# Test de primalité
# Python v2.7
 
from random import randint
 
#############################################################################
def estpetitpremier(n):
    """dit si un nombre n<=1024 est premier"""
    return n in [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,1009,1013,1019,1021]
 
#############################################################################
def lsqrt(x):
    """Racine carrée entière d'un nombre entier x (méth. Héron d'Alexandrie)"""
    # trouver une valeur approchée de la racine pour aller plus vite
    i = 0
    while (x>>i)>0: # trouver l'indice du bit le plus fort
        i += 1 
    r1 = 1 << (((i+1)>>1)-1)
    # calcul de la racine en partant de la racine approchée r1
    while True:
        r2 = (r1+x//r1)>>1 
        if abs(r1-r2) < 2:
            if r1*r1 <= x and (r1+1)*(r1+1) > x:
                return r1
        r1 = r2
 
#############################################################################
def estpremierdiv(n):
    """dit si un nombre est premier par la méthode des divisions """
    # si n==1 ou n==pair => ne renvoyer True que pour n==2
    if n==1 or (n & 1)==0:
        return n==2
    # autre cas: on teste tous les nombres impairs à partir de 3
    k=3
    r=lsqrt(n)
    while k<=r:
        if n % k == 0:
            return False
        k+=2
    return True
 
#############################################################################
def _millerRabin(a, n):
    """fonction utilitaire. Ne pas appeler directement => appeler millerRabin"""
    # trouver s et d pour transformer n-1 en (2**s)*d
    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1
    # calculer l'exponentiation modulaire (a**d)%n
    apow = pow(a,d,n) # =(a**d)%n
    # si (a**d) % n ==1 => n est probablement 1er
    if apow == 1:
        return True
    for r in xrange(0,s):
        # si a**(d*(2**r)) % n == (n-1) => n est probablement 1er
        if pow(a,d,n) == n-1:
            return True
        d *= 2
    return False
 
#############################################################################
def millerRabin(n, k=20):
    """Test de primalité probabiliste de Miller-Rabin
       à utiliser de préférence quand n>100 millions
       k = nombre de tests à faire pour limiter les "faux positifs"
       utilise la fonction utilitaire: "_millerRabin"
    """
    # si n==1 ou n==pair => ne renvoyer True que pour n==2
    if n==1 or (n & 1)==0:
        return n==2
    # recommencer le test k fois: seul les nb ayant réussi k fois seront True
    for repete in xrange(0, k):
        # trouver un nombre au hasard entre 1 et n-1 (bornes inclues)
        a = randint(1, n-1)
        # si le test echoue une seule fois => n est composé
        if not _millerRabin(a, n):
            return False
    # n a réussi les k tests => il est "probablement premier"
    return True
 
#############################################################################
def estpremier(n):
    """dit si n est premier, en utilisant la méthode adaptée à la taille de n"""
    # si n<1024 => recherche dans une liste
    if n<1024:
        return estpetitpremier(n)
    # si n<100 millions => recherche par la méthode des divisions
    if n<100000000:
        return estpremierdiv(n)
    # sinon => recherche par le test probabiliste de Miller-Rabin
    return millerRabin(n, k=20)
 
#############################################################################
def genpremier_bbs(nb):
    """génère un nombre premier de nb bits, congru à 3 modulo 4"""
    p1 = 1 << (nb-1) # ex pour nb=8 bits: int('10000000', 2) = 128
    p2 = (p1 << 1)-1 # ex pour nb=8 bits: int('11111111', 2) = 255
    while True:
        p = randint(p1, p2)
        if (p-3)%4==0 and estpremier(p):
            return p

Utilisation: calculer 2 nombres premiers de 24 bits (en moins de 1/1000 de seconde!):

np1 = genpremier_bbs(24)
21439307
np2 = genpremier_bbs(24)
19013711

Grâce au test de primalité probabiliste de Miller-Rabin, on peut faire beaucoup mieux: Voilà 2 nombres premiers de …512 bits, calculés en quelques secondes (typiquement moins de 5 secondes):

np1 = genpremier_bbs(512)
14889955028678998228128809491538085572595354855300108148273587724116803436124655102599750654713732760112845637394204141228033860781535350246536112810743291
np2 = genpremier_bbs(512)
19323592535440214565041225886718274424436684585050206777038241984399868553805819858163825122973836817287441334977521615639643536088579097250021033947376771

Et donc, leur produit m = np1*np2 sera impressionnant:

print np1*np2
287727423845221976148781429438006511662739674238561051464833561044370017465704962067561076398062137747948585913210676375204119705097371076208207350543791619895824275118754318225416172636176745811393141141711884050017011766526096895591212541694770174764265839983999580130947994713008751252234603173316137493361

Et rien n'empêche d'aller plus loin avec des nombres premiers de 1024 bits (obtenus en typiquement moins de 30 secondes) dont le produit de 2048 bits sera quasi impossible à factoriser avec les techniques actuelles:

np1 = genpremier_bbs(1024)
126372531954418782670747930963559908367043074176552978854728812339634981881297738186484473148885678679817250057381314722686293252256406164004406590313221431594658364724234232601401479859842407134408936812192579280935450252285976584619444811530638558312046094061538487405595898712527059751636380445254298891031
np2 = genpremier_bbs(1024)
123555046844417096837492760317210404666514985905974312915561461451032289102818991294313133675920913246356864032351746814305842372838661180375179088442987931118501244894316663138810505256806261996615220224557941959117257174831629245837169465651034746618434770124678300829334605342426353004782298329345860198811
pnp = np1*np2
15613964105475809149101856910076623309061665516829150898003657641548758372119130773220742570326648288440464650696363563962218616699785594529796973690739204928012833486461687431221425900688370093408218430569798106175136011852868770344375989591408285286072458551713907700373387026966011961633646082159782644472215672024144524954674028008051764350113725668783612255520181138771971248581193864024826475777003802153835548515397602250896961506283420930928960373446251318305080664906123225112150284604088594565795943116506544346028045261052234630606963476566870359641285441058986974252799623555884333020188903103398484764141

Il faut cependant avoir conscience qu'on manipule dans ce dernier cas des nombres de plus de 600 chiffres! La sécurité obtenue se paiera par un plus grand temps de calcul.

Certains articles, mais pas tous, sur ce générateur indiquent que l'on améliore sa qualité, et en particulier la longueur de sa période, en choisissant une graine telle que: PGCD(gr, pnp)=1.

Voilà un code qui fera ça très bien. Il nécessite l'importation de randint du module random, mais on pourrait aussi utiliser genentier_bbs() avec pour graine le nombre de secondes depuis l'epoch (voir plus haut):

from random import randint
 
#############################################################################
def pgcd(a,b):
    """pgcd(a,b): 'Plus Grand Commun Diviseur' de a et b (nb entiers)"""
    while b<>0:
        a,b=b,a%b
    return a
 
#############################################################################
def gengraine_bbs(pnp, g1=1, g2=None):
    """retourne une graine au hasard entre g1 et g2, telle que pgcd(g,pnp)==1"""
    if g2==None:
        g2 = pnp
    while True:
        g = randint(g1, g2)
        if pgcd(g, pnp)==1:
            return g

L'utilisation est facile:

np1 = 24375763
np2 = 28972763
gr = gengraine_bbs(np1*np2)
print gr
471089100952622

Par défaut, gr est choisi entre 1 et pnp=np1*np2. Si on veut une graine plus petite, on peut le préciser à l'appel en fixant x1 et x2:

np1 = 24375763
np2 = 28972763
gr = gengraine_bbs(np1*np2, 10000, 99999)
print gr
16249

Pour cette évaluation, on va travailler sur une liste d'1 million d'octets, dont on va évaluer la qualité de la répartition des 8 millions de bits.

On va de plus travailler avec des nombres premiers de 512 bits et une graine à 6 chiffres:

np1 = genpremier_bbs(512)
# 14889955028678998228128809491538085572595354855300108148273587724116803436124655102599750654713732760112845637394204141228033860781535350246536112810743291
np2 = genpremier_bbs(512)
# 19323592535440214565041225886718274424436684585050206777038241984399868553805819858163825122973836817287441334977521615639643536088579097250021033947376771
pnp = np1*np2
graine = gengraine_bbs(pnp, 100000, 999999)
# 792819

On va aussi voir ce que ça donne avec des nombres premiers plus petits comme:

np1 = genpremier_bbs(64)
# 9331871689144683103
np2 = genpremier_bbs(64)
# 16946680743821830423
pnp = np1*np2
# 158144250258244299340119084951739442569
graine = gengraine_bbs(pnp, 100000, 999999)
# 125309

1- Cet algorithme est considéré comme “lent”. Mais lent à quel point?

Avec la fonction genalea_bbs(graine, pnp) ci-dessus et les valeurs choisies, on obtient cette liste de un million de nombres au hasard (de 0 à 255) en à peu près 1 minute. Cela fait tout de même 0.00006 seconde par nombre au hasard!

Quand on prend des nombres premiers plus petit, ça s'améliore! Par exemple, avec des nombres premiers de 64 bits, ce qui n'est déjà pas mal pour une utilisation courante, on obtient cette même liste d'un million de nombres au hasard en une dizaine de secondes, soit 0.000011 seconde par nombre au hasard, chacun de ces nombres étant composés de 8 bits.

En résumé, on pourra obtenir un nombre au hasard de 0 à 255 en moins de 1/10000 seconde: pas mal, non pour un algorithme “lent”?

2- Est-ce que dans cette liste de un million de nombres, tous les nombres ont une chance égale d'être tirés?

Comme le phénomène générateur est le tirage de bits, on va travailler directement sur les 8 millions de bits.

Quand on compte le nombre de bits à '1' et ceux à '0', on constate qu'il y en a le même nombre à 0.1% près!

On pourrait confirmer avec un test statistique (Test du χ²) pour vérifier que l'écart qui reste est bien du au hasard et non à un défaut d'algorithme.

3- Est-ce que chaque chiffre tiré peut être suivi de n'importe quel chiffre?

Il s'agit de vérifier la résistance de la série à la prévision: si on connait un nombre, il ne faut pas qu'on puisse “deviner” le suivant.

On va travailler ici aussi sur les bits: on va vérifier:

• que chaque '0' peut être suivi d'un '0' ou d'un '1' et que leurs fréquences d'apparition sont égales
• que chaque '1' peut être suivi d'un '0' ou d'un '1' et que leurs fréquences d'apparition sont égales

Quand on fait ça, on trouve des écarts de fréquence < 0.1%: c'est très bon!

4- Est-ce qu'on a une période assez longue?

Comme les séries pseudo-aléatoires sont le plus souvent périodiques, on vérifie que la période est assez longue.

En alignant tous les 8 millions de bits dans une même chaine de caractères, on vérifie qu'on ne retrouve pas les 50 premiers caractères dans la suite de la chaine.

En recommençant de nombreuses fois, même avec les nombres premiers de 64 bits, je n'ai pas retrouvé une seule fois cette séquence de départ. La période d'extraction des bits dans ces conditions est donc supérieure à 8 millions!

Finalement, il n'est pas mal ce générateur!

On va créer un cas concret: crypter et décrypter un fichier contenant n'importe quoi (texte, image, musique, etc…) en faisant une rotation de chaque octet de k bits, k étant donné par une série aléatoire donnée par la fonction gen3bits_bbs décrite plus haut. Cette série donnera une liste de valeurs entières comprises entre 0 et 7 (bornes comprises).

Puisqu'on est en crypto, on va calculer 2 nombres premiers très grands: 512 bits, ce qui donnera un produit de …1024 bits (environ 300 chiffres!). C'est à dire un nombre très difficilement factorisable avec les techniques actuelles. On prendra une graine à 6 chiffres:

    np1 = 14889955028678998228128809491538085572595354855300108148273587724116803436124655102599750654713732760112845637394204141228033860781535350246536112810743291
    np2 = 19323592535440214565041225886718274424436684585050206777038241984399868553805819858163825122973836817287441334977521615639643536088579097250021033947376771
    pnp = np1*np2
    graine = 792819

pour crypter et décrypter, voilà les codes proposés:

#############################################################################
# rotation du mot b (de taille n bits), de k bits à droite (k de 0 à 7)
rotd = lambda b, k=1, n=8: ((b>>k)|(b<<(n-k)))&((1<<n)-1)
 
#############################################################################
def encrypte_bbs(srce, dest, graine, pnp, lgbuf=8192):
    """crypte les octets du fichier srce par rotation de k bits à droite
       k (0 à 7) est tiré au hasard selon le générateur GenaleaBBS(graine, pnp)
       et enregistre le résultat dans le fichier dest
       les fichiers srce et dest sont donnés avec leurs chemins
       lgbuf = taille du buffer de lecture (=8192 octets par défaut)
     """
    gen = gen3bits_bbs(graine, pnp)
    # traitement du fichier srce et enregistrement dans le fichier dest    
    fs = open(srce, 'rb')
    fd = open(dest, 'wb')
    while True:
        bufs = fs.read(lgbuf)
        if bufs=="":
            break # =fin de lecture du fichier source
        bufd = []
        for oct in bufs:
            k = gen.next() # valeur suivante de rotation au hasard
            bufd.append(chr(rotd(ord(oct),k))) # rotation de k bits à droite
        fd.write(''.join(bufd))
    fd.close()
    fs.close()
 
#############################################################################
# rotation du mot b (de taille n bits), de k bits à gauche (k de 0 à 7)
rotg = lambda b, k=1, n=8: ((b<<k)|(b>>(n-k)))&((1<<n)-1)
 
#############################################################################
def decrypte_bbs(srce, dest, graine, pnp, lgbuf=8192):
    """decrypte les octets du fichier srce par rotation de k bits à gauche
       k (0 à 7) est tiré au hasard selon le générateur GenaleaBBS(graine, pnp)
       et enregistre le résultat dans le fichier dest
       les fichiers srce et dest sont donnés avec leurs chemins
       lgbuf = taille du buffer de lecture (=8192 octets par défaut)
     """
    # initialise le générateur de nombres pseudo-aléatoires
    gen = gen3bits_bbs(graine, pnp)
    # traitement et enregistrement    
    fs = open(srce, 'rb')
    fd = open(dest, 'wb')
    while True:
        bufs = fs.read(lgbuf)
        if bufs=="":
            break # =fin de lecture du fichier source
        bufd = []
        for oct in bufs:
            k = gen.next() # valeur suivante de rotation au hasard
            bufd.append(chr(rotg(ord(oct),k))) # rotation de k bits à gauche
        fd.write(''.join(bufd))
    fd.close()
    fs.close()

Et voilà mon code d'essai:

np1 = 14889955028678998228128809491538085572595354855300108148273587724116803436124655102599750654713732760112845637394204141228033860781535350246536112810743291
np2 = 19323592535440214565041225886718274424436684585050206777038241984399868553805819858163825122973836817287441334977521615639643536088579097250021033947376771
pnp = np1*np2
graine = 792819
 
fichier1 = "DSC_0258.JPG" # image de 3 944 490 octets
fichier2 = "DSC_0258_crypt.JPG"
fichier3 = "DSC_0258_decrypt.JPG"
 
print( "crypte")
t = clock()
encrypte_bbs(fichier1, fichier2, graine, pnp)
t = clock()-t
print( "fin de crypte (temps =", t, "sec.)")
 
print( "decrypte")
t = clock()
decrypte_bbs(fichier2, fichier3, graine, pnp)
t = clock()-t
print( "fin de decrypte (temps =", t, "sec.)")

Comme on le voit dans le code, chaque octet subit une rotation à droite pour le cryptage et à gauche pour le décryptage. La valeur de la rotation [0..7] est la même dans les 2 cas, grâce au fait que le générateur donne exactement la même série en démarrant à partir des mêmes 2 nombres (graine et pnp).

Résultat: appliqué à une photo de 3853 ko (3944490 octets), le cryptage puis le décryptage ont produit les effets attendus: le fichier crypté est illisible, et le fichier décrypté reproduit fidèlement la photo initiale.

Chacune des opérations de cryptage et de décryptage n'a demandé que 2 minutes. Il faut voir la taille des opérations: on manipule tout de même des nombres de 300 chiffres pour chaque extraction de bits! Si on avait plusieurs fichiers à crypter/décrypter de la même façon, on pourrait calculer d'avance la liste des rotations aléatoires.

Avec des nombres premiers de 1024 bits (⇒ produit 2048 bits!), le temps de cryptage et de décryptage monte à environ 7 minutes pour traiter ce fichier de 3853 ko, ce qui n'est pas mal quand on manipule des nombres de 600 chiffres….

Avec des nombres premiers de 64 bits, le temps de cryptage et de décryptage descend à environ 20 secondes pour traiter ce fichier de 3853 ko.

Comment ça marche dans les échanges?

• l'émetteur et le destinataire se mettent d'accord sur l'algorithme (ici, celui de “Blum Blum Shub”)
• L'émetteur envoie au destinataire les 2 nombres (la graine et le produit des 2 nombres premiers) par un moyen sûr (transmission cryptée ou courrier).
• Il crypte ensuite les fichiers souhaités et les envoie au destinataire
• le destinataire décrypte les fichiers

Est-ce que c'est un cryptage “solide”? On a vu plus haut que l'algorithme BBS produit de très bonnes séries aléatoires, et qu'il est peu vraisemblable que quelqu'un puisse retrouver les nombres utilisés pour produire la série (dont un nombre de 2048 bits = environ 600 chiffres). Donc, oui, c'est du solide. Mais on a vu des mathématiciens tellement astucieux sur ce genre de sujet…

Amusez-vous bien!