Przemierzając kiedyś ziemskie połacie w zachwycie i zadumie minąłem rzeźbę, której fotografię znajdziemy na Wikipedii. Posiadam całą kolekcję fotografii wszelkich posągów zasłużonych myślicieli, lecz fotografii tej konkretnej, niepozornej statui mi brakuje. Stałem tak chwilę w zadumie i nie mogłem skojarzyć... Leonardo da Pisa, posąg inny niż wszystkie. Cóż, dopiero później myśl wyparła bezmyśl, jednak w tym mieście pełnym uliczek i historii nie mogłem odnaleźć już miejsca, w którym ta rzeźba się znajduje. Szukałem długo, a lokalna społeczność usilnie kierowała mnie do innego posągu.
Analizowany w artykule kod dostępny jest w repozytorium git.
Problem 2: Even Fibonacci numbers
http://projecteuler.net/problem=2
Each new term in the Fibonacci sequence is generated by adding the previous two terms.
By starting with 1 and 2, the first 10 terms will be:
1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...
By considering the terms in the Fibonacci sequence whose values do not exceed four million,
find the sum of the even-valued terms.
Naszym zadaniem jest znalezienie sumy tych elementów ciągu, które są parzyste, a same elementy ciągu mają nie przekraczać zakresu 4 * 10 ** 6 (4 miliony). "Nie przekraczać" czyli: mniejsze lub równe.
Strategia
Problem można rozwiązać iteracyjnie lub rekurencyjnie. Zrobimy tak i tak. Ogólna zasada jednak jest taka: jeśli problem ma rozwiązanie iteracyjne, rozwiązujemy iteracyjnie. W przypadku rekurencji rozwijamy stos, zatem nie każdy problem da się rozwiązać rekurencyjnie z uwagi na ograniczenia maszyn (oraz naszego czasu).
Iteracja:
Dla przypomnienia: iteracja to przeskakiwanie po kolejnych elementach.
>>> L = [1, 2, 3]
>>> for i in L: print(i)
...
1
2
3
Rekurencja:
Rekurencja to samopowtarzalność/samodefinicja. Rekurencję łatwo zobaczyć. Bierzemy telefon, włączamy przedni aparat (ten do selfie) i kierujemy telefon ekranem w stronę lustra. W lustrze oglądamy co wyświetla. Pojawi się wtedy coś równie hipnotycznego jak Ummagumma.
Funkcja jest rekurencyjna kiedy wywołuje samą siebie.
>>> def f(): f()
...
>>> f()
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "<stdin>", line 1, in f
File "<stdin>", line 1, in f
File "<stdin>", line 1, in f
[Previous line repeated 996 more times]
RecursionError: maximum recursion depth exceeded
Jeśli ktoś woli...:
>>> import sys
>>> sys.setrecursionlimit(16)
>>> def f(level=0):
... print("#" * level)
... f(level + 1) # recursive!
...
>>> f()
#
##
###
####
#####
######
#######
########
#########
##########
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "<stdin>", line 3, in f
File "<stdin>", line 3, in f
File "<stdin>", line 3, in f
[Previous line repeated 8 more times]
File "<stdin>", line 2, in f
RecursionError: maximum recursion depth exceeded while calling a Python object
Leonardo myślał o królikach. Let's follow the white rabbit.
Python
Na Pythonie się najłatwiej uczy i ćwiczy.
Python: iteracyjnie
Iteracyjnie - czyli z pętlą. Użyjemy jednak generatora (słowo kluczowe: yield). Generator fibonacci_gen() produkuje kolejne elementy ciągu. Moglibyśmy umieścić w nim klauzulę warunkową if sprawdzającą parzystość, jednak możemy to również odfiltrować w późniejszym etapie. Dzięki temu można tego generatora użyć w innych programach bez zmiany kodu.
#!/usr/bin/env python
import sys
def fibonacci_gen(limit):
prev = 0
curr = 1
item = 0
while item <= limit:
item = prev + curr
yield item
prev = curr
curr = item
limit = int((sys.argv[1] if len(sys.argv) > 1 else (4 * 10 ** 6)))
print(sum([_ for _ in fibonacci_gen(limit) if _ % 2 == 0]))
W ciągu Leonarda mamy liczby:
[0, 1, 1, 2, 3, 5, 8, 13, ...]
W algorytmie śledzimy więc poprzednią i obecną liczbę, dodajemy je i sumę tę za każdym razem wydajemy (yield). Limit dla ciągu możemy podać przy wywołaniu tego skryptu (int(sys.argv[1])). Sumę elementów parzystych załatwiamy w ostatniej linii kodu. Sumujemy iterowane elementy, ale pod warunkiem, że są parzyste.
Underscore (_) tutaj to najzwyklejsza nazwa zmiennej. Programiści Python stosują tę dziwną nazwę w celu zaznaczenia, iż nie jest ona istotą zagadnienia, ani też nie będziemy jej później używać. Moglibyśmy ją nazwać i, x, n, item, czy element_of_a_squence. Nie ma to jednak żadnego znaczenia, takoż _.
Sprawdzamy program:
$ python3 iterative.py
4613732
Mamy więc już wynik. Zapamiętujemy sobie wartość, pozostałe rozwiązania mają dać taki sam wynik.
Z ciekawości możemy też sprawdzić sumy dla większego zakresu - wszakże odbieramy w kodzie argument z linii poleceń.
$ python3 iterative.py 1000000000
1485607536
$ python3 iterative.py 100000000000000000000000000000
171679151392093647435137529168
Możemy też policzyć sumę parzystych elementów ciągu Fibonacciego w zakresie szacunkowej liczby wszystkich cząstek we Wszechświecie.
$ python iterative.py $(expr '3.28 * 10 ^ 80' | bc -l | cut -d. -f 1 | tr -d '[:punct:][:cntrl:]')
442259386782137518167658571404249416738352889129833053963717692207734295884996994
bc umie liczyć, ale formatuje wyniki.
$ expr '3.28 * 10 ^ 80' | bc -l
32800000000000000000000000000000000000000000000000000000000000000000\
0000000000000.00
W tym przypadku akurat tego nie chcemy. Nie interesuje nas też część ułamkowa (dlatego cut). tr usuwa później backslash i znaki kontrolne (tutaj newline).
Teraz "na czysto":
$ expr '3.28 * 10 ^ 80' | bc -l | cut -d. -f 1 | tr -d '[:punct:][:cntrl:]'
328000000000000000000000000000000000000000000000000000000000000000000000000000000
Całość to argument naszego programu. W tym zakresie policzyliśmy powyżej sumę parzystych elementów ciągu. Impressive.
Python: rekurencyjnie
Tutaj Project Euler szykuje już dla nas niespodzianki i zaprasza do eksploracji głębokiego i pięknego świata wiedzy i odkryć. Omówimy zatem rozwiązanie, a okazję wykorzystamy by odkryć odrobinę teorii.
Rekurencyjna postać funkcji ciągu wygląda tak:
def fib_rec(n):
# F(n) = F(n-1) + F(n-2)
if n <= 1:
return n
return fib_rec(n - 1) + fib_rec(n - 2)
"Idziemy" od maximum, a każdy element ciągu to suma wyniku tej funkcji dla poprzednika i jego poprzednika. Kiedy dojdziemy do 1, to przestajemy.
Używalibyśmy tej funkcji w ten sposób:
>>> for i in range(10):
... print(fib_rec(i))
...
0
1
1
2
3
5
8
13
21
34
Liczbę 4 milionów przekroczymy przy n == 35.
>>> [fib_rec(i) for i in range(35)]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, ..., 3524578, 5702887]
Parzyste elementy wyfiltrujemy w ten sposób:
>>> list(filter(lambda n: n % 2 == 0, [fib_rec(i) for i in range(35)]))
[0, 2, 8, 34, 144, 610, 2584, 10946, 46368, 196418, 832040, 3524578]
Wynik otrzymamy poprawny:
>>> sum(list(filter(lambda n: n % 2 == 0, [fib_rec(i) for i in range(35)])))
4613732
Ale...:
To rozwiązanie jest bardzo wolne.
>>> t = time.monotonic(); sum(list(filter(lambda n: n % 2 == 0, [fib_rec(i) for i in range(35)]))); time.monotonic() - t
4613732
4.805196590954438
U mnie zajęło 5 sekund. To sporo, gdyż rozwiązanie iteracyjne z generatorem dla limitu równego liczbie wszystkich cząstek we Wszechświecie zajęło 200 mikrosekund:
>>> t = time.monotonic(); sum([_ for _ in fibonacci_gen(3.28 * 10 ** 80) if _ % 2 == 0]); time.monotonic() - t
442259386782137518167658571404249416738352889129833053963717692207734295884996994
0.00020728004164993763
Całkiem nieźle, nie? Kiedyś pisałem w C++ wielowątkowy filtr SMS z blacklistą/whitelistą i łańcuchem warunków - na jeden SMS przypadało średnio 3.5 mikrosekundy. Systemy real time są rozważane w milisekundach i mikrosekundach.
Wracamy do rekurencji. Poniżej nieco dodatkowej teorii.
Stos
Jeśli nasza funkcja rekurencyjna fib_rec() otrzyma dużą wartość liczbową, to tych wywołań wewnątrz będzie mnóstwo. Dla każdego n - 1 oraz każdego n - 2 będziemy wywoływać tę funkcję, a tym samym tworzyć kolejne ramki na stosie. Sprawdźmy ile tych wywołań jest:
def fib_rec(n):
fib_rec.counter += 1
if n <= 1:
return n
return fib_rec(n - 1) + fib_rec(n - 2)
fib_rec.counter = 0 # yes, it's just an object.
fib_rec(32)
print(fib_rec.counter)
Wykonujemy (dla liczby 32, jak wyżej w kodzie):
$ python count_calls.py
7049155
7+ milionów wywołań. Im większa liczba, tym wolniej. Poniżej przykładowa różnica 4 (czyli limit 32 vs 36) i przyrost czasu:
$ time python recursive.py 32
2178309
real 0m0.764s
user 0m0.760s
sys 0m0.004s
$ time python recursive.py 36
14930352
real 0m5.141s
user 0m5.122s
sys 0m0.008s
Ograniczenie stosu:
Jeśli będziemy wykonywać tę funkcję dla dużej wartości n, to python zgłosi wyjątek RecursionError. W przeciwnym wypadku system zabiłby sam proces Python.
Traceback (most recent call last):
File "/tmp/count_calls.py", line 10, in <module>
fib_rec(1024)
File "/tmp/count_calls.py", line 6, in fib_rec
return fib_rec(n - 1) + fib_rec(n - 2)
[Previous line repeated 995 more times]
File "/tmp/count_calls.py", line 4, in fib_rec
if n <= 1:
RecursionError: maximum recursion depth exceeded in comparison
Stos w Python
W przykładzie powyżej Python przerwał naszą funkcję. Co się stanie, jeśli zmienimy ograniczenie rekursji na wyższe?
Najpierw zobaczmy, że to C:
Traceback (most recent call last):
File "/tmp/count_calls.py", line 10, in <module>
sys.setrecursionlimit(102400000000000)
OverflowError: Python int too large to convert to C int
Tutaj inty to inne inty. Ale o tym było wczoraj. Teraz o rekurencji: zmieniamy limit rekursji, liczymy dla n=100k, a w systemie zmieniamy rozmiar stosu na 4096k.
Kod:
import sys
def fib_rec(n):
if n <= 1:
return n
return fib_rec(n - 1) + fib_rec(n - 2)
sys.setrecursionlimit(1234567890)
fib_rec(100000)
Wykonanie:
$ ulimit -s 4096
$ python3 count_calls.py
Segmentation fault (core dumped)
Ups... nadepnęliśmy na paszczę węża. Pozostaje nam pochylić się nad nim i zapytać co się stało:
$ gdb -q $(command -v python3) core.dev.lan.python3.1640658530
[...]
Core was generated by 'python3 count_calls.py'.
Program terminated with signal SIGSEGV, Segmentation fault.
#0 0x000000000050f61e in _PyEval_EvalFrameDefault ()
(gdb) bt
#0 0x000000000050f61e in _PyEval_EvalFrameDefault ()
Backtrace stopped: Cannot access memory at address 0x7ffc7fcd40f8
(gdb) frame view 0x000000000050f61e
#0 0x0000000000000000 in ?? (warning: Unable to restore previously selected frame.
)
Nie każdy problem da się rozwiązać rekurencyjnie. Jeśli problem ma rozwiązanie iteracyjne, to stosujemy rozwiązanie iteracyjne.
Tail call: TRE/TCO
Pisząc funkcje rekurencyjne możemy stosować TRE (Tail Recursion Elimination), a ogólnie to TCO - Tail Call Optimization. Chodzi o to, aby przypadek rekurencyjny (return func()) był na samym końcu.
def func(x):
if condition:
return ...
return func(x + 1) # here! recursive! tail.
To pozwala na eliminację nadmiarowych ramek stosu. W Python tego nie ma, co więcej - zmierzymy to teraz:
import time
def fib_rec(n):
if n <= 1:
return n
return fib_rec(n - 1) + fib_rec(n - 2) # tail call
def fib_rec_no_tail(n):
if n > 1:
return fib_rec(n - 1) + fib_rec(n - 2)
return n
def timeit(func, *args, **kwargs):
t0 = time.monotonic()
func(*args, **kwargs)
dt = time.monotonic() - t0
print(dt, func.__qualname__)
return dt
timeit(fib_rec, 42)
timeit(fib_rec_no_tail, 42)
Wykonuję:
$ python3 calls.py
8.308601832017303 fib_rec (37,) {}
8.29210988804698 fib_rec_no_tail (37,) {}
Jeśli ktoś jest zainteresowany szczegółami, to:
- ktoś ma fajnego bloga i tłumaczy how-tail-call-optimization-works,
- tutaj Guido tłumaczy dlaczego w Python tego nie ma.
Podsumowanie
$ time ./iterative.py
4613732
real 0m0.031s
user 0m0.024s
sys 0m0.008s
$ time ./recursive.py
4613732
real 0m5.093s
user 0m5.075s
sys 0m0.008s
Uuuuuuuu. Pierwsze jest szybkie jak kobra, a to drugie jest jak... Ouroboros. Taki rekurenycjny wąż, który zjada swój własny ogon.
C
O C i C++ mówiliśmy w ramach integralności danych. Tutaj nie będziemy się męczyć. Robimy "na piechotę", iteracyjnie, z ifem w środku. Można też doczytać o niecodziennych dodatkach do formatowania (PRIu64).
Kod:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#define __STDC_FORMAT_MACROS
#include <inttypes.h>
int main(int argc, char** argv)
{
uint64_t maximum = argc > 1 ? strtoull(argv[1], 0, 10) : 4000000;
if(errno) {
fprintf(stderr, strerror(errno));
return 1;
}
uint64_t last = 0;
uint64_t next = 1;
uint64_t result = 0;
uint64_t sum_even = 0;
uint64_t sum_all = 0;
/* printf("# Maximum fib number: %ld\n", maximum); */
while((result = last + next) <= maximum) {
last = next;
next = result;
if(result % 2 == 0) {
sum_even += result;
}
sum_all += result;
/* printf(". %20" PRIu64 "\n", result); */
}
printf("SUM of even: %" PRIu64 "\n", sum_even);
/* printf("TOTAL : %" PRIu64 "\n", sum_all); */
return 0;
}
Kompiluję, uruchamiam, patrzymy:
$ cc -o fib main.c -Wall -Wextra -Werror
$ time ./fib
SUM of even: 4613732
real 0m0.001s
user 0m0.001s
sys 0m0.000s
No... ultra szybka rzecz. Tak jak i w Python - tutaj również dorobiłem możliwość przekazania zakresu (argv[1]). Domyślnie jest 4 miliony, ale jeśli przekażemy argument, to możemy sprawdzić wyższe wartości:
$ time ./fib 12345678901234567890
SUM of even: 15970217317495049952
real 0m0.001s
user 0m0.001s
sys 0m0.000s
Nadal ultra szybki. Sprawdzamy więc dla liczby cząstek we Wszechświecie:
$ ./fib $(expr '3.28 * 10 ^ 80' | bc -l | cut -d. -f 1 | tr -d '[:punct:][:cntrl:]')
Numerical result out of range
Boom. Nie ma. Big Bang. Już z tym walczyliśmy. Tym razem wiemy już, że się da, ale wracamy do Project Euler. Euler to też Leonhard. "Leonard z Bazylei".
Lisp
(defun fib-gen (&optional (max 4000000) (prev 0) (next 1) (acc ()))
(if (> prev max)
(nreverse acc)
(fib-gen max next (+ prev next) (cons prev acc) )))
(prin1 (apply '+ (remove-if-not #'evenp (fib-gen))))
defun to oryginał Pythonowego def. &optional to oryginał Pythonowego *args. cons to construct (tworzy listę, tu wszystko jest listą). apply to to samo co dzisiaj w Pythonie map(). Tutaj funkcją jest +, a remove-if-not to to samo co filter() w Python. Funkcją dla tego filtra jest evenp czyli predykat "czy parzyste?". W środku funkcji fib-gen jest if/else. Tak, tutaj funkcje w nazwach mogą mieć kreski, slashe, itd. Jemu wszystko jedno, pierwsze to funkcja, reszta to argumenty. Troszkę trzeba pogłówkować, ale otwiera oczy i umysły.
Odpalamy:
$ time clisp -q solution.lisp
4613732
real 0m0.008s
user 0m0.004s
sys 0m0.004s
Umie liczyć i nie tylko. Może dzisiaj już troszkę stary (na polskie realia - emeryt), ale wciąż żywy i pewnie będzie z nami kiedy inne już zdążą zardzewieć.
C++
long double już męczyliśmy wcześniej, nie ma co pokazywać.
Można do rozwiązania podejść tak, ale przedmiotowy problem jest natury pętli for/while. To może - dla odmiany, wartości dydaktycznej oraz entuzjastów -funroll-loops - zrobimy "BEZ PĘTLI".
#include <iostream>
int main()
{
int limit = 4000000;
int prev = 0;
int curr = 1;
int temp = 0;
int total = 0;
fib:
temp = prev + curr;
if (temp > limit) goto end;
total += (temp % 2 == 0) ? temp : 0;
prev = curr;
curr = temp;
goto fib;
end:
std::cout << total << std::endl;
}
Kompiluję, uruchamiam, patrzymy:
$ c++ loopless.cxx
$ ./a.out
4613732
Pięknie. Zero myślenia. Zajrzyjmy mu pod maskę:
$ stat -c "%s" a.out
16592
$ echo "disass main" | gdb --quiet a.out > A.txt
$
Włączam optymalizacje. Kompiluję, uruchamiam, patrzymy:
$ c++ loopless.cxx -O3
$ ./a.out
4613732
$ stat -c "%s" a.out
16624
Oo! Większy!
$ echo "disass main" | gdb --quiet a.out > B.txt
$
Zróbmy benchmark na koniec roku 2021.
$ c++ loopless.cxx
$ time for x in `seq 2021`; do ./a.out; done > /dev/null
real 0m2.867s
user 0m2.220s
sys 0m0.794s
$ c++ loopless.cxx -O3
$ time for x in `seq 2021`; do ./a.out; done > /dev/null
real 0m2.878s
user 0m2.182s
sys 0m0.845s
Oo! Większy i wolniejszy. "Faced with a loop optimizer that has some brains, I decided to stop messing around" ...
bc
Tutaj raczej nie ma czego tłumaczyć.
define fib(limit) {
prev = 0;
curr = 1;
temp = 0;
total = 0;
while (temp <= limit) {
temp = prev + curr;
if (temp % 2 == 0) {
total += temp;
}
prev = curr;
curr = temp;
}
return total;
}
print "EULER: ", fib(4000000), "\n"
print "UNIVERSE:\n"
fib(3.28 * 10 ^ 80)
Wow, nawet markdown go lubi.
$ cat fib.bc | bc -q
EULER: 4613732
UNIVERSE:
44225938678213751816765857140424941673835288912983305396371769220773\
4295884996994
Teraz będzie benchmark, zatem bez Universe. Wyłączamy Wszechświat. W poniższej pętli trochę czasu zżera kotek, ale bez niego bc wchodzi w tryb interaktywny. Nie ma co kombinować, sprawdzamy:
$ time for x in `seq 2021`; do cat fib.bc | bc -q; done > /dev/null
real 0m3.849s
user 0m4.556s
sys 0m1.998s
Sam kotek zjada ~2 sekundy, a dochodzi jeszcze exec bc.
$ time for x in `seq 2021`; do cat fib.bc > /dev/null; done
real 0m1.959s
user 0m1.421s
sys 0m0.670s
Kiedy bc nic nie robi, to i tak samo uruchomienie zajmuje dużo czasu:
$ time for x in `seq 2021`; do echo | bc -q > /dev/null; done
real 0m3.892s
user 0m3.526s
sys 0m1.497s
Hmmm, impressive. Podoba mi się. Bierzemy go do następnych odcinków.
Podsumowanie
Zsumowaliśmy parzyste liczby ciągu Leonarda da Pisa. Mam nadzieję, że można było się czegoś ciekawego dowiedzieć. Ja się co nieco dowiedziałem.
W następnym odcinku mamy mnożenie i liczby pierwsze. Yay, can't wait!
Jeśli ktoś dysponuje własną fotografią całego posągu Leonardo da Pisa, będę wdzięczny za podesłanie jej na email wraz z licencją. Mam jednak nadzieję, że kiedyś sam tam wrócę i tym razem go odnajdę.
W tym artykule zostawię zatem jedynie foto dwóch heptadekagonów (tutaj jest oryginał [2.8M]).
Dzisiaj nie było słonika, ale wróci.