Numeri primi in Haskell

Rieccoci con la serie di post sulla programmazione funzionale, finalmente (o purtroppo, dipende dal vostro gradimento). Oggi cercheremo di realizzare un semplice test per i numeri primi in Haskell. Lo scopo è di essere didattici, non di essere efficienti, quindi i matematici e i super esperti portino pazienza, per favore!

Vediamo prima di fare un po di sommario.

Puntate precedenti:

Ripasso di Haskell

Sì, perché l’ultima volta è stata un sacco di tempo fa, e allora facciamo un po’ di ripassino.

I linguaggi funzionali prediligono le funzioni. Tutto è funzione, l’utente definisce funzioni, che possono prendere in ingresso altre funzioni, e ne può fare delle applicazioni parziali per ottenere delle funzioni più semplici. In particolare in Haskell per definire una funzione basta scrivere un’equazione. Ad esempio, ecco la funzione che data una lista ci restituisce un’altra lista contenente solo i numeri pari:

pari :: [Integer] -> [Integer]
pari l = [ x | x <- l, even x ]

Vi ricordo che la prima linea è la type signature, cioè la dichiarazione del tipo di dato che passiamo come argomento in input e del tipo di dato che diamo come output; la seconda linea è invece la definizione vera e propria della funzione. In questo caso abbiamo usato la list comprehension, ricordate? Leggiamolo come se fosse notazione matematica: “tutti gli x tali che x appartiene a l, e x è pari”.

Numeri primi

Per continuare il ripasso, proviamo timidamente a scrivere una funzione che ci dica se un dato intero è primo o no. Per cominciare, ci serve sapere se un numero divide un altro. Si dice che un numero divide esattamente un altro, se il resto della divisione è zero. Quindi, dato n il numero da dividere e x il divisore, deve essere vero che:

n `mod` x == 0

Dove mod è la funzione che calcola il resto. A onor del vero bisognerebbe scrivere

(mod n x) == 0

perché Haskell usa la notazione prefissa. Però chi ha inventato Haskell sa bene che a volte la notazione infissa (quella degli operatori) è più comoda. Per cui ha dato la possibilità di usare entrambe. Se avete una funzione di due argomenti, potete usare la notazione infissa se mettete la funzione tra “backtick”. Quindi, scrivere mod n x è equivalente a scrivere n `mod` x. Chiaro? (se la risposta è no, non vi preoccupate, ci saranno altri esempi nel seguito).

Per fare le cose un po’ più leggibili, definiamoci una funzione chiamata divides:

divides :: Integer -> Integer -> Bool
x `divides` n = n `mod` x == 0

Ok, adesso vogliamo sapere se un numero è primo. La definizione dice che un numero n è primo se è divisibile solo per se stesso e per l’unità. Quindi, dobbiamo provare a dividere n per tutti i numeri da 2 a n-1, e vedere se viene fuori che qualcuno è divisibile.

La lista di numeri da 2 a n-1 si scrive semplicemente [2..n-1]. La funzione divides ci restituisce dei valori booleani, ovvero True oppure False. Se io applico la funzione divides a tutti i membri della lista [2..n-1] ottengo una lista di True e False. Se tutti i valori della lista risultante sono False, allora il numero è primo, mentre se almeno uno è True, allora il numero è composto. Ci siamo fino a qui? Bene, guardiamo il codice che esprime questo mio ragionamento:

isprime :: Integer -> Bool
isprime n = not $ or [ x `divides` n | x <- [2..n-1] ] 

Spiegazione, cominciando dal fondo: qui vedete una bella list comprehension per costruire la lista di booleani (tutti i risultati della funzione x `divides` n tali che x appartiene alla lista [2..n-1]). La funzione or, come ci si aspetterebbe, prende una lista di booleani e ne fa l’or logico. La funzione not restituisce la negazione logica del suo input e ci da il risultato. Direi che ci siamo… no aspetta, e quel dollaro? Il dollaro serve per evitare di mettere troppe parentesi. Infatti, Haskell legge da sinistra verso destra, e se non mettete le parentesi potrebbe incasinarsi: il dollaro serve a dirgli “aspetta che arrivi il risultato di quello che sta a destra del dollaro prima di valutare la funzione che sta a sinistra del dollaro”. Tecnicamente: cambia l’ordine di precedenza delle funzioni, o meglio associa a destra invece che a sinistra.

Ed ecco il programmino risultante:

module Main where
import System.Environment

main :: IO()
main = do args <- getArgs
          putStrLn ( show $ isprime (read $ args!!0) )

divides :: Integer -> Integer -> Bool
x `divides` n = n `mod` x == 0

isprime :: Integer -> Bool
isprime n = not $ or [ x `divides` n | x <- [2..n-1] ] 

Per provare il codice: salvate il testo qui sopra in un file chiamato ‘primo.hs’ (oppure come meglio credete). Quindi, compilate il file da linea di comando con :

ghc primo.hs

e otterrete il file eseguibile primo, che potete eseguire così:

./primo 1237

e dovrebbe restituirvi True.

(Esercizio: che cosa succede se togliete il dollaro senza mettere alcuna parentesi? provare per capire…)

Let e where

Un matematico storcerà il naso a vedere il mio algoritmo. Inefficiente! Infatti, per provare che un numero è primo, è sufficiente testare i numeri da 2 fino alla radice quadrata di n. E allora facciamolo, che aspettiamo?

Per fare la radice quadrata di n si usa la funzione sqtr che però prende un double e restituisce un double (insomma, quasi). Mentre il nostro n è intero: e quindi se scrivo sqrt n il compilatore protesta: come vi ho già detto, Haskell è fortemente, maledettamente tipato, e se non gli dai i tipi giusti non ti compila. E quindi, bisogna trasformare n in un double, e poi il risultato di sqrt di nuovo in intero. Ecco come si fa:

floor $sqrt $fromIntegral n 

Al solito, i dollari evitano le parentesi.

Potremmo mettere l’espressione con la radice quadrata direttamente in linea al posto di n-1, ma sarebbe orribile a leggersi. Potremmo definire una nuova funzione:

sqn n = floor $sqrt $fromIntegral n
isprime n = not $ or [ x `divides` n | x <- [2..sqn n] ] 

ma definire una nuova funzione sembra eccessivo, specialmente se questa funzione viene usata solo dalla funzione isprime. E allora usiamo il let:

isprime :: Integer -> Bool
isprime n = let sqn = floor $sqrt $fromIntegral n  
            in not $ or [ x `divides` n | x <- [2..sqn] ] 

Let serve a definire delle funzioni “locali” alla funzione isprime. Si usa come nella notazione matematica, e quindi si legge così: “sia (let) sqn la radice quadrata di n, nell’espressione not …”. Quindi, dopo il let si mettono le definizioni locali (anche più di una, su righe diverse), quindi la parola chiave in, e infine si mette l’espressione che usa le definizioni precedenti. Altro esempio:

mytan x = let s = sin x
              c = cos x
          in s / c

definisce la tangente di un angolo. Invece di let si può usare un costrutto alternativo che si chiama where:

mytan x = s / c
          where s = sin x
                c = cos x

che si legge: “la funzione mytan x è definita come s / c, dove (where) s è definito come il seno di x, e c è definito come il coseno di x”. Si legge bene, vero?

C’è una differenza importante fra usare let e where. Let può essere messo dovunque ci sia un’espressione, e quindi anche dentro un’espressione più grande. Per esempio, riscriviamo la nostra sqn

sqn n = floor (let i = fromIntegral n in sqrt i)

Come potete vedere, abbiamo inserito una let/in in un’espressione più grande. La where invece può solo essere legata a una funzione, e quindi la seguente cosa dà un errore di sintassi:

sqn n = floor (sqrt i where i = fromIntegral n)

Se proprio volete usare il where, dovete scrivere:

sqn n = floor s 
        where s = sqrt (fromIntegral n)

Lo so, non è chiarissimo, ma nel seguito spero di inserire qualche altro esempio.

Il secondo programmino:

module Main where
import System.Environment

main :: IO()
main = do args <- getArgs
          putStrLn ( show $ isprime (read $ args!!0) )

x `divides` n = n `mod` x == 0

isprime :: Integer -> Bool
isprime n = let sqn = floor $sqrt $fromIntegral n
            in not $ or [ x `divides` n | x <- [2..sqn] ] 

Al solito, fermatevi e prendetevi 30 secondi di tempo per compilare, eseguire, etc. E magari, provate a usare il where invece del let.

Ricorsione!

Adesso siamo un po’ più efficienti, ma ancora lontani da una buona soluzione. Per esempio, possiamo controllare la divisibilità solo per i primi contenuti tra 2 e la radice quadrata di n. Infatti, non ha senso controllare la divibilità per 4 una volta che il numero non è più divisibile per 2. Ed è inutile controllare la divisibilità per 9 se il numero non è divisibile per 3, ecc.

Ma la lista di primi fino alla radice quadrata di n non ce l’abbiamo! Idea: non potremmo utilizzare la stessa definizione di isprime per calcolarci questa lista?

Allora, adesso proviamo a ragionare top/down. Abbiamo già una funzione isprime decente (anche se non ottimizzata), proviamo ad utilizzarla per calcolare la lista (infinita) di tutti i primi. Come vi ho spiegato qualche tempo fa Haskell permette di definire liste infinite, ovvero liste che potenzialmente contengono un numero infinito di elementi. Haskell è lazy, e quindi non è che effettua un calcolo infinito: semplicemente, calcola gli elementi man mano che servono. Ecco la lista infinita di primi a partire dal 2:

primes :: [Integer]
primes = (filter isprime [2..])

Abbiamo usato la funzione filter che prende come argomenti una funzione (in questo caso isprime), e una lista (in questo caso la lista infinita degli interi maggiori o uguali a 2), e applica la funzione a tutti gli elementi della lista. Il risultato è una lista che contiene tutti e solo gli elementi per i quali isprime ha restituito True, ovvero la lista dei primi.

Ci siamo quasi. Adesso, per testare la primalità, ci basta controllare la divisibilità per tutti i numeri primi fino alla readice quadrata di n. Ecco il codice:

isprime n = let sqn = floor $sqrt $fromIntegral n
            in not $ or [ x `divides` n | x <- takeWhile (<=sqn) primes ]

Spiegazione: la funzione takeWhile prende due paramteri, una funzione e una lista. La prima funzione è una sezione dell’operatore <= e ci restituisce True quando un numero è minore o uguale di sqn. Il secondo argomento è una lista, nel nostro caso la lista infinita primes. La funzione takeWhile applica la funzione a tutti gli elementi della lista finché la funzione restituisce True, poi si ferma. Il risultato è quindi tutta la prima parte della lista fino a che troviamo un numero pari a sqn.

Vi si sta attorcigliando il cervello? No problem, accade sempre con casi di ricorsione doppia incrociata con scappellamento a destra! Qui abbiamo che la funzione isprime ha bisogno della lista primes per essere calcolata, e quest’ultima si calcola a partire dalla funzione isprime. Proviamo a vedere se funziona:

module Main where
import System.Environment

main :: IO()
main = do args <- getArgs
          putStrLn ( show $ isprime (read $ args!!0) )

x `divides` n = n `mod` x == 0

isprime :: Integer -> Bool
isprime n = let sqn = floor $sqrt $fromIntegral n
            in not $ or [ x `divides` n | x <- takeWhile (<=sqn) primes ] 

primes :: [Integer]
primes = (filter isprime [2..])

Salvate in un file chiamiato primi2.hs, compilate e eseguite. Ed ecco quello salta fuori:

lipari@lipari$ ./primi 11
primi: <<Loop>>

Ahi ahi. Che è successo? E’ successo che Haskell si è incasinato e non riesce a venire fuori da questa ricorsione doppia con lista infinita. In particolare, quando ha cercato di calcolare il primo elemento della lista primes, è andato in loop. Non avremmo esagerato? Beh, sì un po. Vediamo perché e magari salterà fuori che abbiamo capito qualcosa in più di come funziona il tutto.

  1. Per prima cosa, Haskell cerca di calcolare la fuzione isprime sul numero che gli abbiamo passato sulla linea di comando, che è 11. Per fare questo cerca di risolvere la list comprehension, e in particolare cerca di estrarre il primo elemento della lista primes per darlo in pasto alla takeWhile.
  2. Bisogna quindi calcolare il primo elemento (solo il primo!) della lista primes, che corrisponde al primo elemento della lista risultante dalla funzione filter
  3. Per fare questo, deve prima calcolare isprime 2.
  4. Torniamo su: per calcolare isprime 2 ha bisogno del primo elemento della primes…
  5. Haskell a questo punto si accorge di essere finito in un ciclo infinito, e alza le braccia, stampando l’errore sul terminale e uscendo dal programma.

L’errore qui è che lui non riesce a partire, non ha un punto di inizio. In matematica, si direbbe che abbiamo un’induzione senza una base; in informatica, abbiamo una ricorsione senza una condizione di terminazione. Notate che Haskell non vi pianta il calcolatore fino al fatidico out of memory. Si accorge del loop infinito, e esce subito.

La soluzione è più semplice di quanto possiate immaginare. Se Haskell non riesce a calcolarsi il primo elemento della lista primes, possiamo suggerirlo noi! La definizione di primes diventa:

primes = 2:(filter isprime [3..])

Il primo della lista lo mettiamo noi, gli altri glieli lasciamo calcolare. Adesso abbiamo la base induttiva.

Proviamo a rifare il ragionamento di prima. Dopo il punto 2, il punto 3 ci dice che il primo elemento di primes è 2, che è minore di sqn. Quindi, x prende il valore 2. 11 non è divisibile per 2, quindi il primo elemento della lista da dare in pasto alla funzione or è un bel False. A questo punto dobbiamo calcolare il secondo elemento di primes. Per far questo dobbiamo calcolare isprime 3. La radice di 3 è sqn=1.73…, quindi l’espressione takeWhile (<sqn) primes è una lista vuota. L’or di una lista vuota da come risultato False, e quindi isprime 3 è True. Quindi il secondo elemento di primes è 3, e possiamo tornare a calcolare isprimes di 11. 3 non divide 11, quindi adesso la lista di booleani contiene due False, e ci tocca andare a calcolare il terzo elemento di primes… che penso ci possiamo fermare qui.

A questo punto, dopo tutto questo sforzo mentale, provate a compilare e lanciare il programma!

Ultima piccola modifica. Nel calcolare la lista dei primi, noi sappiamo già che i pari non sono primi, e quindi potremmo applicare una piccola ottimizzazione:

primes = 2:(filter isprime [3,5..])

La lista [3,5,..] è un modo elegante e molto haskelliano di rappresentare i numeri dispari da 3 in su.

Conclusioni

Naturalmente, se proprio vogliamo calcolare la lista dei numeri primi, ci sono metodi parecchio più efficienti, come il crivello di Eratostene, ma questo magari lo vediamo la prossima volta che oggi penso di avervi incriccato la cervice abbastanza. E chissà, magari, forse, prima o poi diamo uno sguardo anche al crivello di Atkin (ah, gli studenti di oggi!)

Note

L’immagine è presa da questo blog.

Annunci
Post a comment or leave a trackback: Trackback URL.

Commenti

  • zar  On 15 febbraio 2012 at 22:38

    Uh, ricorsività incrociata!

  • juhan  On 16 febbraio 2012 at 10:44

    Alla fine il programma è breve ma certo che la sintassi è molto diversa a quella dei linguaggi normali.
    Ma sembra anche leggibile, se e quando ti abitui. Bello il $ anche se le parentesi non hanno mai spaventato i lispisti. A proposito: inizialmete McCarthy aveva proposto anche la notazione infissa ma poi è stata abbandonata (anche per via della comodità di trattare la lista in funzione del primo elemento (car)) (ecco vedi due paretesi!). –OK, OT 😉

    • glipari  On 17 febbraio 2012 at 14:45

      Ma in realtà non è tanto OT!
      Uno dei problemi in effetti è superare lo scoglio iniziale della sintassi. Un’altro problema grosso è capire come funziona il type system di Haskell, che non è tanto banale. Finora l’ho nascosto sotto il tappeto, ma si avvicina il momento in cui dovrò parlare della monade, e di come si scrive codice “sequenziale” in haskell… ma intanto godiamoci la ricorsione incrociata.

  • peppe, blablabla... (@9peppe)  On 14 maggio 2014 at 09:42

    (c’è baco. secondo quel programma 9 è primo, hai dimenticato il <=)

    • glipari  On 14 maggio 2014 at 10:41

      Acc, l’avevo messo nel primo listato e mi è saltato sul secondo! grazie mille per il debug!

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione / Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione / Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione / Modifica )

Google+ photo

Stai commentando usando il tuo account Google+. Chiudi sessione / Modifica )

Connessione a %s...

%d blogger hanno fatto clic su Mi Piace per questo: