Progetto Eulero: Problema 3

1. si provano in ordine i numeri primi finché non si trova un divisore del numero da fattorizzare;

   13195 | 2 3 5
   

2. questo è uno dei fattori primi. Si effettua dunque una divisione del numero da fattorizzare per il divisore trovato;

   13195 | 2 3 5
    2639 |
   

3. si ricomincia l'algoritmo ripartendo dal risultato della divisione.

   13195 | 2 3 5
    2639 | 5 7
     377 | 7 11 13
      29 | 13 17 19 23 29
       1 |
   

Questo algoritmo estremamente semplice necessita però della conoscenza di quali siano i numeri primi, un problema antico ma attualissimo. Cogliamo quindi lo spunto per parlare dei numeri primi e della loro ricerca.

I numeri primi

I numeri primi sono il sottoinsieme dei numeri naturali caratterizzati dalla proprietà di essere divisibili solo per sé stessi e per l'unità. Per definizione, inoltre, 1 non è un numero primo. I primi numeri primi sono dunque:

L'ultimo numero primo è... nessuno. Non esiste il massimo numero primo. La dimostrazione, dovuta ad Euclide, è molto semplice ma affascinante: supponiamo per assurdo che i numeri primi siano finiti. Possiamo dunque elencarli tutti:

Ma allora possiamo calcolare un numero p ottenuto moltiplicando tutti i numeri della lista e aggiungendo 1:

questo numero è sicuramente non incluso nella lista, perché p > n_m. Inoltre non è divisibile per nessun numero primo, infatti, per qualunque numero n della lista, p mod n = 1. Dunque n è primo, ma non incluso nella lista, il che rende assurda l'ipotesi di averli elencati tutti.

Con questa osservazione si potrebbe pensare di aver trovato un metodo per generare i numeri primi: se li si conoscono tutti fino a un certo punto, basterebbe moltiplicarli tutti ed aggiungere uno... et voila! Purtroppo questo non funziona: sappiamo solo che p è primo, ma è probabile che non sia il primo tra quelli non ancora generati. Lo si può verificare osservando i primi elementi della lista: 1 + (2 × 3) = 7: il 5 non è stato generato.

Uno degli algoritmi per generare numeri primi è forse uno degli algoritmi più antichi che si conoscano: il Crivello di Eratostene.

Il Crivello di Eratostene

L'algoritmo risale al cirenaico Eratostene: matematico, astronomo e tanto altro, come lo erano i geni dell'antichità. Tra le sue opere fu notevole anche la misura della circonferenza della terra nel III secolo A.C.

L'idea del crivello è quella di eliminare i multipli dei numeri primi trovati: i numeri che restano sono anche essi primi. Diciamo ad esempio di voler conoscere la primalità dei numeri fino al 50: si parte dalla condizione seguente:

    2  3  4  5  6  7  8  9 10
11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30
31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50

1 è escluso per definizione. 2 è primo, dunque cancello dal crivello (con un piolo, il massimo della potenza di calcolo disponibile nel 250 A.C...) tutti i suoi multipli:

    2  3  .  5  .  7  .  9  .
11  . 13  . 15  . 17  . 19  .
21  . 23  . 25  . 27  . 29  .
31  . 33  . 35  . 37  . 39  .
41  . 43  . 45  . 47  . 49  .

Il numero successivo rimasto è 3, ciò vuol dire che è un numero primo: cancello anche i suoi multipli giungendo alla situazione:

    2  3  .  5  .  7  .  .  .
11  . 13  .  .  . 17  . 19  .
 .  . 23  . 25  .  .  . 29  .
31  .  .  . 35  . 37  .  .  .
41  . 43  .  .  . 47  . 49  .

che mi dice che 5 è il successivo numero primo (che mi permette di eliminare il 25 e il 35). Giunti al 7 ci rendiamo conto che esso elimina solo il 49: tutti i multipli precedenti sono già stati eliminati da un precedente fattore primo, per esempio 21 = 3 × 7 è stato già eliminato dal 3. 49 = 7 × 7 ha il 7 come minimo (ed unico) fattore primo: in generale, per trovare tutti i numeri primi fino ad n è sufficiente crivellare i multipli dei primi inferiori alla radice quadrata di n. Lo schema è ora completo:

    2  3  .  5  .  7  .  .  .
11  . 13  .  .  . 17  . 19  .
 .  . 23  .  .  .  .  . 29  .
31  .  .  .  .  . 37  .  .  .
41  . 43  .  .  . 47  .  .  .

La seguente è una semplice implementazione in Python del crivello, non orientata alle prestazioni ma alla chiarezza:

from math import sqrt, ceil

def crivello(n):
    """Restituisci la lista dei numeri primi < n."""

    c = [ True ] * n    # crivello: i numeri cancellati saranno False
    primi = []          # valori da restituire

    # Prova i numeri fino alla radice di n.
    for i in range(2, int(ceil(sqrt(n)))):
        if c[i]:
            # Hai trovato un primo: aggiungilo alla lista...
            primi.append(i)

            # ...e cancella tutti i suoi multipli.
            for j in range(i * i, n, i):
                c[j] = False

    # I numeri rimasti sono tutti primi.
    for i in range(i + 1, n):
        if c[i]:
            primi.append(i)

    return primi

Questa funzione restituisce una lista di numeri primi inferiori ad un massimo scelto:

>>> print crivello(50)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

il che può andar bene... se ci volessimo accontentare. Ma se non si conosce a priori il numero di primi desiderato? Se dopo averne generati un milione ne volessimo uno in più, dovremmo ripetere la generazione dall'inizio. Se invece fossero bastati i primi 20 avremmo sprecato un sacco di tempo. Non sarebbe elegante un generatore di numeri primi su richiesta, senza limiti prefissati?

I generatori

Un generatore è un oggetto in grado di produrre dati su richiesta. I risultati vengono generati e restituiti man mano che servono al cliente: il generatore deve pensare a mantenere lo stato della propria computazione.

Senza l'aiuto da parte del linguaggio, fare entrambe le cose contemporaneamente (mantenere lo stato, produrre risultati su richiesta) può essere complesso: una tipica soluzione è quella di passare ad una funzione produttrice una funzione callback che sarà invocata ogni volta che un nuovo dato sarà pronto (è il cosiddetto Principio di Hollywood: non chiamarci: ti chiameremo noi!). La soluzione è realizzabile in tutti i linguaggi che hanno almeno il concetto di puntatore a funzione, ma porta facilmente a complicare il flusso del codice. L'alternativa di pre-computare tutti i valori e fornirli su richiesta è spesso poco pratica (per esempio se dovessimo generare dati a seguito del parsing di un file enorme) o del tutto impossibile (se il generatore è in grado di fornire infiniti risultati).

Molti moderni linguaggi di programmazione forniscono un supporto esplicito per la costruzione di generatori. Riprendiamo ad esempio un problema fresco e costruiamo, in Python, un generatore della successione di Fibonacci:

def genera_fibonacci():
    """Genera i numeri della successione di Fibonacci."""
    a, b = 1, 1
    while True:
        yield a
        a, b = b, a + b

La funzione genera_fibonacci(), quando invocata, costruisce una nuova istanza di generatore, ma non genera alcun numero. Ogni volta che al generatore viene richiesto un nuovo elemento, il codice della funzione viene eseguito fino all'istruzione yield, che restituisce al chiamante un valore. Ma, a differenza di un convenzionale return, la funzione non viene conclusa: il valore delle variabili locali e la posizione corrente vengono mantenute nello stato del generatore. Quando viene richiesto il valore successivo, l'esecuzione riprende dall'istruzione successiva allo yield, con il contenuto delle variabili locali immutato. La generazione di numeri può terminare con un return oppure uscendo dallo scope della funzione. Nell'esempio della successione di Fibonacci ovviamente ciò non può succedere e il generatore continuerà a produrre tutti i numeri che il cliente gli chiederà.

I numeri possono essere generati ad esempio dando in pasto il generatore ad un consumatore for each:

>>> for n in genera_fibonacci():
...     if n > 100: break
...     print n,
...
1 1 2 3 5 8 13 21 34 55 89

I generatori sono apparsi in alcuni linguaggi di programmazione dagli anni '70 (tra questi nel CLU e nell'Icon) e sono attualmente disponibili in numerosi linguaggi di concezione moderna (per esempio in Ruby e in C#)

Un generatore infinito di numeri primi

L'idea alla base del generatore infinito di primi è quella di non conservare i "buchi" del crivello in una matrice, ma di associare ad ogni numero primo trovato un "piolo", e di memorizzare per ogni piolo la posizione corrispondente al prossimo suo multiplo da analizzare, facendoli avanzare man mano che serve.

Possiamo rappresentare lo stato del crivello con una lista ordinata di coppie (m, n), dove ogni coppia rappresenta un piolo: n è il numero primo che il piolo rappresenta, m è il prossimo numero, multiplo di n, che il piolo cancellerà.

Il primo elemento della lista, ovvero quello con m più piccolo, sarà il prossimo non-primo da non restituire. Per esempio, supponendo di aver generato i numeri primi 2, 3, 5, 7, lo stato del crivello sarà:

ci troviamo quindi ad esaminare l'8, che vediamo essere multiplo di 2. Di conseguenza non restituiremo tale numero e aggiorneremo le posizioni a:

Neanche il 9 è primo, quindi:

e non lo è neanche il 10:

Ora tocca all'11: poiché non compare nel primo elemento della lista ne deduciamo che si tratta un numero primo. Possiamo quindi emetterlo ed aggiungere un piolo per rappresentarlo. Già abbiamo osservato che la prima volta che ci sarà necessario usarlo sarà per valutare il suo quadrato, mentre ai suoi multipli precedenti ci penseranno gli altri pioli:

Come si vede, al passo successivo dovremo spostare 2 pioli per superare il buco in posizione 12, in quanto sia 2 che 3 sono suoi fattori primi.

In termini di prestazioni, l'aggiornamento del crivello può essere pesante: l'operazione più frequente, ovvero lo spostamento dei pioli, costa o(n) sia implementando il crivello con una lista linkata (perché dovremmo farne una scansione potenzialmente completa per trovare la nuova posizione), sia implementandolo con un array (essendo ordinato, potremmo trovare il punto di inserimento in o(log(n)) passi procedendo per bisezione, ma poi dovremmo spostare fino ad n elementi per posizionare il piolo). Si può fare di meglio?

La struttura di dati heap

Un heap (in particolare un min heap) è un albero in cui il valore dei figli di un nodo è sempre maggiore del valore del nodo (un max heap gode della proprietà inversa). In un heap binario inoltre tutti i livelli devono essere pieni, tranne al più l'ultimo che deve essere riempito da sinistra a destra. L'albero nell'immagine è un esempio di heap binario.

Figura 1.

L'aggiunta di un elemento all'heap si chiama push: per effettuarla, prima si mette il nuovo nodo (nell'esempio il "5") temporaneamente in coda all'albero,

Figura 2. Inserimento di un elemento nell'heap.

quindi si inverte la sua posizione con quella del padre finché la caratteristica dell'heap non tornerà valida, ovvero ci fermeremo quando il valore del padre sarà minore del valore del nodo.

Figura 3.

In maniera simile si può rimuovere l'elemento minimo dello heap (operazione di pop): prima si sostituisce la radice con l'ultimo nodo,

Figura 4. Rimozione di un elemento dall'heap.

quindi si fa percolare il nodo dalla radice fino ad un punto in cui i suoi figli avranno valore maggiore, tornando dunque ad un heap valido.

Figura 5.

In entrambi i casi, il numero delle inversioni da fare è proporzionale alla profondità dell'albero, quindi il loro costo è o(log(n)). Ovviamente la ricerca del valore minimo è o(1): il minimo è sempre nella radice dell'albero. Queste caratteristiche lo rendono ottimo per costruire code con priorità, o per lo scheduling di eventi (gli eventi futuri possono essere inseriti in un heap: il prossimo che accadrà è sempre quello nella radice).

Poiché le operazioni nello heap sono costruite sull'inversione di coppie di elementi, è facile ed efficiente memorizzare i nodi dello heap in un array: i figli dell'elemento in posizione k sono sempre nelle posizioni (2 k + 1) e (2 k + 2) (considerando k = 0 per la prima posizione). L'albero dell'ultimo esempio è rappresentato dall'array:

L'heap è una struttura dati ben studiata, ma spesso non ne sono ovvi gli utilizzi. Si possono trovare implementazioni dello heap in tutti i linguaggi di programmazione.

Avrete già capito che lo heap fa esattamente al caso nostro: come in una lista ordinata, il primo elemento è sempre il minore ed è quello da valutare per stabilire la primalità del numero in esame. Ma lo spostamento dei pioli può essere realizzato con soli o(log(n)) scambi di posizione.

Mettiamo tutto insieme...

Questa è una semplice implementazione di un generatore infinito di numeri primi in linguaggio Python.

from heapq import heappush, heapreplace

def genera_primi():
    """Generatore infinito di numeri primi."""

    # Generiamo il primo valore a parte, per non dover gestire il caso
    # particolare in cui l'heap è vuoto, che si presenta solo ora.
    yield 2

    # heap di coppie (numero da eliminare, numero primo eliminatore).
    todel = [ (4, 2) ]

    n = 3
    while True:
        if todel[0][0] != n:
            # Il numero attuale non è la testa dell'heap: il numero è primo!
            yield n
            heappush(todel, (n*n, n))   # nuovo piolo

        else:
            # Il numero non è primo: sposta tutti i pioli dei suoi fattori.
            while todel[0][0] == n:
                p = todel[0][1]
                heapreplace(todel, (n+p, p))
                # heapreplace è una versione ottimizzata di un pop
                # del minimo seguito da un push, per cui la dimensione
                # dell'heap resta invariata.

        n += 1

Questa implementazione è ancora passibile di numerose ottimizzazioni: per esempio si può far avanzare n a 2 a 2, saltando direttamente tutti i numeri pari e quindi evitando di muovere il più ballerino dei pioli, quello del 2. Si può generalizzare questa ottimizzazione prendendo in considerazione k primi minimi e setacciando i numeri facendo avanzare n secondo un pattern di salti pre-computato. I margini di miglioramento sono ampi ma, ancora una volta, in questo esempio si è privilegiata la chiarezza. Inutile dire anche che più raffinati metodi sono stati sviluppati per generare numeri primi (ad esempio il crivello di Atkin), che però non hanno ancora risolto quello che è considerato il Sacro Graal della Teoria dei Numeri, ovvero una formula per calcolare l'i-esimo numero primo.

Con questo strumento è ora facile generare i fattori primi di un numero: per restare in tema anche questo sarà un generatore, questa volta però destinato ad esaurire gli elementi generati. Implementa l'algoritmo visto all'inizio dell'articolo, per cui il numero da fattorizzare viene diviso per i suoi divisori man mano che essi vengono trovati.

def genera_fattori_primi(n):
    """Genera i fattori primi del numero n."""

    # Sappiamo che il massimo fattore primo non può essere maggiore della
    # radice quadrata del numero. Un caso particolare è quando n è primo:
    # lo gestiremo alla fine.
    top = ceil(sqrt(n))

    for m in genera_primi():

        # Un fattore può essere ripetuto (es. 12 = 2 x 2 x 3):
        # itera finché possibile.
        while not n % m:

            # m è un fattore primo di n: genera il numero e dividi m.
            yield m
            n //= m

            # Ora dobbiamo cercare i fattori primi di un numero più piccolo:
            # il massimo può essere abbassato.
            top = ceil(sqrt(n))

        if m > top:
            # Abbiamo superato il massimo prestabilito: alla fine
            # dell'esecuzione rimane sempre un numero primo:
            # sarà l'ultimo valore da restutuire.
            if n != 1: yield n

            return  # termine del generatore

Possiamo consumare interamente questo generatore per conoscere il maggiore dei fattori del numero proposto dal problema:

>>> for n in genera_fattori_primi(317584931803): pass
...
>>> print n
3919

L'esecuzione è rapidissima, merito del fatto che, a dispetto dell'enormità del numero iniziale, abbiamo dovuto generare i soli numeri primi fino a 3919, ovvero soli...

>>> print len(list(takewhile(lambda n: n <= 3919, genera_primi())))
543

probabilmente un numero minore di qualunque stima avremmo potuto fare. Questo è il tipo di guadagno ottenibile con la tecnica di valutazione pigra, concetto chiave in alcuni linguaggi marcatamente funzionali quali l'Haskell (grazie alla sua particolare strategia di valutazione delle espressioni, call by need), ma alla portata di molti altri linguaggi grazie all'uso dei generatori.

I numeri primi sono uno degli elementi fondamentali della Teoria dei Numeri, ciononostante molte loro proprietà sfuggono ancora alla comprensione. La loro distribuzione sembra casuale, eppure in essa si celano pattern esoterici e misteriosi, come si vede nell'immagine finale ottenuta disponendo i numeri primi lungo una spirale quadrata. Incontrerete anche quest'ultima più avanti, sul cammino del Progetto Eulero.

Figura 6. Disposizione a spirale dei numeri primi.