Progetto Eulero: Problema 2

Secondo problema di Eulero.La scorsa settimana abbiamo introdotto il primo di quella che sarà una lunga serie di articoli riguardanti Project Euler. Ringrazio coloro che hanno bloggato circa l’iniziativa o hanno commentato direttamente nel sito, mostrandoci soluzioni persino in GW-Basic e Lua.

Fatti i dovuti ringraziamenti, passiamo subito al problema numero 2, così enunciato:

Ogni nuovo termine nella sequenza di Fibonacci è generato dalla somma dei due termini precedenti. Partendo con 1 e 2, i primi dieci termini saranno:
1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …
Trova la somma di tutti i termini pari nella sequenza, che non eccedono un milione.

Osserviamo per un attimo la definizione matematica dei numeri di Fibonacci:

I numeri di Fibonacci

Un primo approccio naïve è quello di creare una funzione ricorsiva tratta da tale definizione, che restituisca il numero di Fibonacci per un dato indice passato come argomento. A quel punto basterà scrivere un ciclo nel quale verifichiamo che il termine della sequenza è pari aggiungendolo, in tal caso, alla somma. Il ciclo si interrompe nel momento in cui un temine trovato è maggiore di 106. In Ruby abbiamo:

def fib(n)
if n == 0 || n == 1
n
else
fib(n-1) + fib(n-2)
end
end

N = 1_000_000

sum, i = 0, 0
while (val = fib(i)) <= N
sum += val if val % 2 == 0
i += 1
end

puts "Somma: #{sum}"

L’esecuzione dello script richiede di media 7.173 secondi con Ruby 1.8.6 e 1.804 secondi con Ruby 1.9 (sul mio MacBook Pro 2.2 GHz). La versione corrente di Ruby (1.9) fornisce il risultato in circa un quarto del tempo, ma in entrambi i casi è facile intuire che per valori di N molto grandi, l’algoritmo adottato non è in grado di fornire un risultato nei 60 secondi previsti dalle regole del Progetto Eulero. Ad esempio anche solo per N = 108, Ruby 1.9 impiega già più di due minuti. E se volessimo calcolare la somma per N che vale 1010000? La soluzione vista finora sarebbe del tutto inadeguata.

Può essere provato che la funzione fib così definita, opera in tempo O(ϕ^n), dove ϕ (phi) è la sezione aurea:

La sezione aurea

In altre parole, la complessità computazionale cresce esponenzialmente al crescere dell’argomento. La crescita esponenziale del call stack generato è un collo di bottiglia indipendentemente dal linguaggio adottato. Il primo passo è dunque quello di ottimizzare la funzione utilizzata sino ad ora. Questo può essere ottenuto salvando all’interno di variabili gli ultimi due valori calcolati, per poi sommarli e generare un nuovo valore, ripetendo il ciclo finché necessario (a seconda di N). La parte relativa alla verifica del numero e la sua somma, rimane inalterata. Sempre in Ruby, questa diventa:

N = 1_000_000

a, b, sum = 1, 2, 0
while (val = a + b) <= N
a, b = b, val
sum += val if val % 2 == 0
end

puts "Somma: #{sum}"

Ruby 1.9 lo esegue in qualche dozzina di milionesimi di secondo, rispetto ai quasi due secondi precedenti. Provando ad imporre un limite di 1010000, che prima era inconcepibile, otteniamo un risultato corretto in soli 1.602 secondi, e stiamo sempre parlando di un linguaggio notoriamente lento. Un algoritmo così efficiente ci permette di sperimentare anche con le nuove caratteristiche del linguaggio, ottenendo comunque risultati in tempi accettabili. Ad esempio, lo script qui sopra può essere riscritto utilizzando la classe Fiber introdotta con Ruby 1.9, ed esegue in circa 1.7 secondi per N = 1010000. Questo il codice che risolve il problema iniziale per N = 106:

fib = Fiber.new do
f1 = f2 = 1
loop do
Fiber.yield f1
f1, f2 = f2, f1 + f2
end
end

N = 1_000_000

sum = 0
while (val = fib.resume) <= N
sum += val if val % 2 == 0
end

puts "Somma: #{sum}"

Soluzioni analoghe possono essere scritte in vari linguaggi con facilità. Ad esempio, in Python possiamo usare il seguente snippet che è praticamente identico a quello in Ruby.

a, b, sum = 1, 2, 0
while b <= 1000000:
if b % 2 == 0:
sum += b
a, b = b, a + b
print sum

In Lisp abbiamo:

(defun p ()
(let* ((a 0)
(b 1)
(fib (lambda ()
(prog1 a (psetf a b b (+ a b))))))
(loop for i = (funcall fib)
until (> i 1000000)
when (evenp i)
sum i)))

E la sua versione ottimizzata (che richiede due milionesimi di secondo):

(defun p2-optimized ()
(declare (optimize (speed 3) (debug 0) (safety 0)))
(let* ((a 0)
(b 1)
(fib (lambda ()
(the fixnum (prog1 a (psetf a b b (the fixnum (+ a b))))))))
(declare (fixnum a b))
(loop for i = (the fixnum (funcall fib))
until (> i 1000000)
when (evenp i)
sum i
of-type fixnum)))

Nel caso dei termini della successione che non eccedono 1010000, Python e Lisp (la prima versione) richiedono rispettivamente 1.326 e 0.264 secondi.

Per i più interessati, ecco una soluzione per il linguaggio funzionale Erlang:

-module(p2).
-compile(export_all).

fib(N) ->
{Fib, _} = fib(N, {1, 1}, {0, 1}), Fib.

fib(0, _, Pair) -> Pair;
fib(N, {Fib1, Fib2}, Pair) when N rem 2 == 0 ->
SquareFib1 = Fib1*Fib1,
fib(N div 2, {2*Fib1*Fib2 - SquareFib1, SquareFib1 + Fib2*Fib2}, Pair);
fib(N, {FibA1, FibA2}=Pair, {FibB1, FibB2}) ->
fib(N-1, Pair, {FibA1*FibB2 + FibB1*(FibA2 - FibA1), FibA1*FibB1 +
FibA2*FibB2}).

solve(N) ->
solve(0, 0, N).
solve(Fib, Num, Lim) when Fib < Lim, Fib rem 2 == 0 ->
Cur = fib(Num),
Fib + solve(Cur, Num + 1, Lim);
solve(Fib, Num, Lim) when Fib < Lim ->
Cur = fib(Num),
solve(Cur, Num + 1, Lim);
solve(Fib, Num, Lim) ->
0.

Per Haskell:

p2 = sum [ x | x <- takeWhile (<= 1000000) fibs, x `mod` 2 == 0]
where fibs = 1 : 1 : zipWith (+) fibs (tail fibs)

main = print p2

In C:

#include 

int main() {
 int sum  = 0; 
 int curr = 1; 
 int succ = 2;
 int max  = 1000000;
 int temp;
 while (curr <= max) {
   if((curr % 2) == 0) sum += curr;
   temp = curr;
   curr = succ;
   succ += temp;
 }
 printf("Somma: %d\n", sum);
 return 0;
}

In C#:

using System;

class Problema2
{
static void Main()
{
int sum = 0;
int curr = 1;
int succ = 2;
int max = 1000000;
int temp;
while (curr <= max)
{
if (curr % 2 == 0) sum += curr;
temp = curr;
curr = succ;
succ += temp;
}
Console.WriteLine("Somma: {0}", sum);
}
}

E infine in Standard ML:

fun sumEvenFib (a, b, c) =
if b > 1000000 then c
else if b mod 2 = 0 then sumEvenFib (b, a+b, c+b) else
sumEvenFib(b, a+b, c)

val () = print (Int.toString (sumEvenFib (1, 2, 0)))

Abbiamo dimostrato nuovamente che un algoritmo efficiente è la chiave di volta per il programmatore, soprattutto quando questi è alle prese con problemi le cui soluzioni più immediate e naïve operano con una complessità computazionale di tipo esponenziale, come è spesso il caso per problemi di tipo matematico (che in Project Euler abbondano). Ciò che ci ha permesso di ottimizzare così tanto la funzione iniziale è stato, nel nostro caso, il passaggio da un algoritmo ricorsivo ad uno iterativo. Potremmo fermarci qui, in fondo abbiamo trovato un algoritmo soddisfacente. Qualcuno nella mailing list del team ha però scritto: “Non ho trovato teoremi particolari per questo problema, esiste un teorema per ottenere un numero in una data posizione nella sequenza ma non la loro somma e men che meno la somma di solo quelli pari.”

In matematica quando il teorema non c’è, basta inventarlo. La formula a cui si riferiva il redattore che ha posto la domanda è la cosiddetta formula di Binet:

Dove phi è sempre la sezione aurea citata più sopra. Da questa formula, come possiamo trovare la somma dei soli termini pari? Osserviamo per un attimo la successione dei primi 20 numeri di Fibonacci:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, …

Notate nessun pattern nella distribuzione dei numeri pari? Facendo un po’ di attenzione si può notare che i termini pari sono presenti in posizioni multiple di 3, F(0) = 0, F(3) = 2, F(6) = 8 e così via. Possiamo dunque modificare la formula di Binet per identificare in maniera univoca solamente i membri pari della sequenza, scrivendola in funzione di x, dove n = 3x per qualsiasi valore di x in N (0 incluso):

Binet per numeri pari

A questo punto possiamo calcolare la sommatoria per i che va da 0 a (un generico) n ed ottenere la seguente formula:

Formula sommatoria dei termini pari di Fibonacci

La formula trovata ci permette di inserire la posizione di un termine nella sequenza di Fibonacci ed ottenere la somma di tutti i numeri pari presenti fino a quel punto; in maniera aritmetica, per cui risolvibile con carta e penna o comunque con una calcolatrice non programmabile.

Per risolvere aritmicamente il problema, ci rimane solo da determinare la posizione del termine che non eccede il milione. Utilizzando la formula di Binet modificata per i numeri pari è possibile stimare (è un equazione trascendente) che per x = 11 si ha il primo termine maggiore di 106. Useremo dunque nella formula trovata qui sopra, 10 come valore della variabile n, per ottenere lo stesso risultato calcolato inizialmente nei vari linguaggi proposti.

Il teorema (o meglio la formula) non c’era, ma l’abbiamo scovata. Rimane però una considerazione da fare. Anche senza calcolare la serie e ricercare matematicamente il valore di n che non fa eccedere alla formula di Binet il milione del problema iniziale, è comunque possibile usare un semplice script che impieghi la formula di Binet che ho modificato più sopra, dopo avere osservato il pattern dei termini positivi. In questo modo sarà il programma ad eseguire numericamente il calcolo della sommatoria per noi e a decidere “empiricamente” quando fermarsi. L’unico svantaggio di questo approccio ibrido è quello di avere a che fare con numeri in virgola mobile piuttosto che interi e calcolare una formula del genere è piuttosto “costoso”, per cui il vantaggio di una soluzione semi-aritmetica viene sminuito da questi fattori. Ovviamente la portata di questi svantaggi dipende sopratutto dal linguaggio scelto.

In ogni caso, l’algoritmo può essere implementato in Ruby abbastanza semplicemente, nel seguente modo:

include Math

def f(x)
(((1 + sqrt(5))/2)**(3*x) - ((1 - sqrt(5))/2)**(3*x))/sqrt(5)
end

N = 1_000_000

sum = i = 0
while (val = f(i)) <= N
sum += val
i += 1
end

puts "Somma: #{sum}"

Lascio ai lettori, come esercizio, l’implementazione di un algoritmo che utilizzi il primo approccio iterativo mostrato inizialmente, sfruttando però l’osservazione circa la posizione, a multipli di tre, dei termini pari della successione di Fibonacci.

Alla prossima sfida.

About Antonio Cangiano

Antonio lavora come Software Engineer & Technical Evangelist presso la IBM in Canada. È inoltre il direttore di Stacktrace.it, un internet marketing strategist, imprenditore del web, serial blogger, e autore di un paio di libri in inglese (recentemente Technical Blogging.) Puoi dare un'occhiata ai suoi progetti sulla sua homepage e seguendolo su Twitter.

Comments

  1. è una figata avervi nel feed reader 😉
    siete troppo hardcore xD

  2. Cavolo, mi impressionano molto quei 2 milionesimi di secondo in Lisp…

  3. Questo era davvero un problema semplice ma dai risvolti parecchio interessanti. Personalmente non conoscevo nemmeno la formula di Binet.

    Siete troppo forti!

  4. Gabriele says:

    Non c’è che dire, continuate così..Vi distinguete e appassionate!

  5. Sapevo della formula di Binet (che però credo sia un poco più “costosa” scritta in codice), ma sinceramente non avevo pensato di poter calcolare direttamente i numeri pari della successione in quel modo: notevole.
    Comunque la soluzione in Ruby 1.9 usando la nuova classe Fiber è molto simile all’utilizzo delle coroutine in Lua, anche perché all’atto pratico sono la stessa cosa (si tratta di semi-coroutine in entrambi i casi ed espongono un’interfaccia simile), comunque ecco come ho risolto il problema proprio in Lua:

    function fibonacci_sum(maxf)
        local fibonacci = coroutine.wrap(function ()
            local x, y = 0, 1
            while true do 
                coroutine.yield(x)
                x, y = y, x + y
            end
        end)
    
        local sum = 0
        for n in fibonacci do 
            if n < = maxf then
                if n % 2 == 0 then sum = sum + n end
            else 
                break
            end 
        end
        return sum
    end
    
    print("Somma: " .. fibonacci_sum(1000000))
    

    Rimanendo in tema Ruby (sempre 1.9, ma per pochi dettagli) voglio riportare la versione del codice *volutamente* meno autoesplicativa a cui sono giunto solo ed esclusivamente per il gusto di farmi del male

    def fibonacci_sum(max)
        sum, x, y = 0, 0, 1
        sum += x if x.even? while (proc ||= -> { x, y = y, x + y; x }).() < = max
        sum
    end
    
    puts "Somma: #{fibonacci_sum 10 ** 6}"
    

    Magari può suscitare curiosità quel proc ||= {...}, ma nella astrusità di quella singola riga c'è un suo motivo (se avete voglia di trovarlo, anche se OT... :))

    Ora volta però aggiungo anche F#, recentemente mi sono messo a giocare un po' con questo linguaggio e lo trovo decisamente simpatico (inoltre verrà ufficialmente supportato in VS nel prossimo futuro) per cui non potevo sfruttare occasione migliore per continuare, anche perché sono abbastanza soddisfatto del risultato che alla fine ho ottenuto (il tutto in forma "lazy"):

    #light
    let fibonacci_sum max = 
        (0, 1) 
            |> Seq.unfold (fun (x, y) -> 
                    if x < = max then Some (x, (y, x + y)) else None
                )
            |> Seq.filter (fun x -> x % 2 = 0)
            |> Seq.fold (+) 0
    
    printfn "%d" (fibonacci_sum 1000000)
    

    Per finire (e in PHP per far contento Fullo :p) la somma con la formula di Binet originale ma usando un passo incrementato a multipli di 3:

    function fibonacci_nth($x) {
        return (pow((1 + sqrt(5)) / 2, $x) - pow((1 - sqrt(5)) / 2, $x)) / sqrt(5);
    }
    
    function fibonacci_sum($max) {
        for ($n = $i = 0; $n < 1000000; $i += 3, $n = fibonacci_nth($i))
            $sum += $n;
        return $sum;
    }
    
    printf("Somma: %d", fibonacci_sum(1000000));
    

    PS: l'ho già detto vero che sono notevoli questi post? 🙂

  6. non ero a conoscenza del teorema di binet. In ogni caso avevo trovato la seguente soluzione (in python) iterativa che sfrutta anche la considerazione sulla posizione:

    tot = 0
    x = 0
    y = 1
    limit = 100000
    while True:
    if x%2 == 0 :
    tot+=x
    #se x è pari posso saltare un ciclo
    x,y=y,x+y
    #se y supera il limite esco
    elif y > limit:
    break
    elif y%2 == 0:
    tot+=y
    #se y è pari posso saltare un ciclo
    x,y=y,x+y
    x+=y
    #se x supera il limite esco subito
    if x > limit:
    break
    y+=x

    Con questa soluzione ottengo un tempo di esecuzione di circa 0.37, contro i circa 0.42 che ottengo facendo girare la soluzione python semplice.

  7. Una domanda sulla versione in erlang… ma perchè è stata implementata in quel modo? Non andava bene questo?

    -module(e2).
    -export([solve/0, solve/1, fib/1, fib/3]).

    % versione “semplificata” (non usata)
    fib(N) ->
    fib(N, 0, 1).

    fib(0, A, _) -> A;
    fib(1, A, B) -> A + B;
    fib(N, A, B) -> fib(N-1, A+B, A).

    % versione usata
    solve() ->
    solve(1000000).

    solve(N) ->
    solve(N, 0, 0, 1).

    solve(Max, Sum, A, B) ->
    F = A + B,
    if F < Max ->
    case F rem 2 of
    0 -> solve(Max, Sum+F, B, F); % pari -> sommo
    1 -> solve(Max, Sum, B, F) % dispari -> non sommo
    end;
    true ->
    Sum
    end.

    (spero nell’indentazione…)

  8. (…appunto…)

    -module(e2). 
    -export([solve/0, solve/1, fib/1, fib/3]).
    % versione "semplificata" (non usata) 
    fib(N) -> fib(N, 0, 1).
    fib(0, A, _) -> A; 
    fib(1, A, B) -> A + B; 
    fib(N, A, B) -> fib(N-1, A+B, A).
    % versione usata 
    solve() -> solve(1000000).
    solve(N) -> solve(N, 0, 0, 1).
    solve(Max, Sum, A, B) -> F = A + B, if F < Max -> case F rem 2 of 
      0 -> solve(Max, Sum+F, B, F); % pari -> sommo 
      1 -> solve(Max, Sum, B, F) % dispari -> non sommo 
    end; 
    true -> Sum 
    end.
    

    PS: se non funziona neanche questa volta, prego gli admin di cancellare questo messaggio

  9. @Zed: ciao sono io l’autore del codice erlang postato nell’articolo.

    Premetto che non sono nativo di erlang :), il mio linguaggio primario è ruby e mi diletto soltanto in erlang.

    Il codice relativo alla sequenza di fibonacci è complesso perché è cercavo la massima efficienza, ed ho utilizzato l’algoritmo più efficiente per determinare i numeri di fibonacci. Più informazioni a riguardo su literateprograms.

    Per quanto riguarda solve non mi piace usare case (che sicuramente è il modo più leggibile di procedere), quindi ho fatto un’ulteriore funzione sfruttando il pattern matching.

  10. mmm… magari non so calcolare i tempi di una funzione, ma…

    e2.erl

    -module(e2).
    -export([timeit/1]).
    
    timeit(N) ->
        printTime(giovanni, fun solve/1, N),
        printTime('   zed', fun solveZ/1, N).
    
    printTime(Atom, Fun, Par) ->
        Start = erlang:now(),
        OutPut = Fun(Par),
        Stop = erlang:now(),
        io:format("~w(~w) = ~w\t[time = ~w]\n", [Atom, Par, OutPut, timer:now_diff(Stop, Start)]).
    
    % ZeD
    solveZ(N) ->
        solveZ(N, 0, 0, 1).
    
    solveZ(Max, Sum, A, B) ->
        F = A + B,
        if F < Max ->
            case F rem 2 of
                0 -> solveZ(Max, Sum+F, B, F); % pari -> sommo
                1 -> solveZ(Max, Sum, B, F)    % dispari -> non sommo
            end;
        true ->
            Sum
        end.
    
    
    % Giovanni
    fib(N) ->
       {Fib, _} = fib(N, {1, 1}, {0, 1}), Fib.
    
    fib(0, _, Pair) -> Pair;
    fib(N, {Fib1, Fib2}, Pair) when N rem 2 == 0 ->
       SquareFib1 = Fib1*Fib1,
       fib(N div 2, {2*Fib1*Fib2 - SquareFib1, SquareFib1 + Fib2*Fib2}, Pair);
    fib(N, {FibA1, FibA2}=Pair, {FibB1, FibB2}) ->
       fib(N-1, Pair, {FibA1*FibB2 + FibB1*(FibA2 - FibA1), FibA1*FibB1 +
    FibA2*FibB2}).
    
    solve(N) ->
           solve(0, 0, N).
    solve(Fib, Num, Lim) when Fib < Lim, Fib rem 2 == 0 ->
           Cur = fib(Num),
           Fib + solve(Cur, Num + 1, Lim);
    solve(Fib, Num, Lim) when Fib < Lim ->
           Cur = fib(Num),
           solve(Cur, Num + 1, Lim);
    solve(Fib, Num, Lim) ->
           0.
    

    dalla console

    67> c('e2').
    ./e2.erl:50: Warning: variable 'Fib' is unused
    ./e2.erl:50: Warning: variable 'Lim' is unused
    ./e2.erl:50: Warning: variable 'Num' is unused
    {ok,e2}
    68> e2:timeit(9999999999999999999999999999999999999999999999999999999).
    giovanni(9999999999999999999999999999999999999999999999999999999) = 8713593760208533336540511604820729774562803052910629256     [time = 5832]
    '   zed'(9999999999999999999999999999999999999999999999999999999) = 8713593760208533336540511604820729774562803052910629256     [time = 330]
    ok
    69> e2:timeit(9999999999999999999999999999999999999999999999999999999).
    giovanni(9999999999999999999999999999999999999999999999999999999) = 8713593760208533336540511604820729774562803052910629256     [time = 5106]
    '   zed'(9999999999999999999999999999999999999999999999999999999) = 8713593760208533336540511604820729774562803052910629256     [time = 329]
    ok
    
  11. (niente da fare, sono tonto e non so come copincollare pezzi di codice… mi sa che la prossima volta uso pastebin o servizi simili)

  12. Antonio Cangiano says:

    Ciao ZeD,

    ho corretto la formattazione. Purtroppo appare in celestino, perché il nostro CMS crede sia stato inviato dall’autore dell’articolo (cioè me), quando in realtà l’ho solo modificato. Lo sistemiamo quanto prima. 🙂

  13. @ZeD: prova il timing esclusivamente delle funzioni per i numeri di fibonacci e cerca un numero di ordinalità alta. In questo modo saranno evidenti le differenze nell’implementazione delle funzioni di fibonacci.

  14. @giovanni
    Scusa, probabilmente sono stato troppo orbo, e solo ora vedo “decentemente” il tuo codice…

    …ma la tua solve ad ogni ciclo ri-calcola da zero l’i-esimo numero di fibonacci, senza salvare da qualche parte l’i-1-esimo e l’i-2-esimo? ora si’, capisco perchè la tua versione è molto più lenta della mia… il fatto è che la versione “più veloce” che hai preso da literateprograms è ottima se devi calcolare solo un generico numero, poichè dimezza i calcoli necessari (ad esempio, per calcolare F(50) non è necessario conoscere F(49)…F(27), bastano F(25) e F(26))

    …ma questo si trasforma in un peso nella tua solve, in quanto la maggior parte dei valori dovrà essere ricalcolata!

  15. Si Zed, è molto più lenta la mia solve, ma non mi sono preoccupato più di tanto perché avevo quel generatore di fibonacci “blazingly fast”. Il bello di Project Euler è che non devi affrontare eccessive ottimizzazioni su tutti i problemi.

    Progredento nei problemi ci si accorge che alcuni necessitano di soluzioni clever, ma questo non era uno di quelli 🙂

  16. Complimenti per il lavoro, bravi. Avrei una domanda, e un piccolissimo appunto. Prima la domanda: non ci sono problemi di approssimazioni utilizzando la formula di Binet? L’appunto, invece, è questo: il termine corretto per l’equazione sarebbe “trascendente”, non “trascendentale”.

    Complimenti ancora!

  17. Grazie per i complimenti, professore.

    Il problema dell’approssimazione nel calcolo in virgola mobile c’e’ e ne limita l’applicabilita’ per numeri sufficientemente grandi. Rimane pero’ possibile l’utilizzo della formula sfruttando librerie a precisione aritmetica arbitraria, ad esempio BigDecimal in Ruby e Java. Le operazioni con queste sono ovviamente piu’ lente di un double/float, per cui bisogna vedere se il gioco vale la candela a seconda del problema che si intende risolvere.

    Grazie anche per l’appunto circa il “trascendentale”. L’ho corretto. Vivendo all’estero da piu’ di quattro anni, l’inglesismo mi scappa ogni tanto. 🙂

Policy per i commenti: Apprezzo moltissimo i vostri commenti, critiche incluse. Per evitare spam e troll, e far rimanere il discorso civile, i commenti sono moderati e prontamente approvati poco dopo il loro invio.