Liczby Fibonacciego - implementacja#
Liczby Fibonacciego to ciąg liczb całkowitych określony następującą formułą rekurencyjną:
\(F_0=0\) , \(F_1 =1 \) and \(F_{n+2} = F_{n+1}+F_{n}\) for \(n\geq 0\)
Ten ciąg znany jest od czasów średniowiecza i pojawia się w wielu nieoczekiwanych miejscach (np. w biologii, w problemach związanych z liczbą złotego podziału itd.). Liczba (stosunek) złotego podziału to \(\frac{1+\sqrt{5}}{2}\) i może być przybliżona jako iloraz dwóch kolejnych liczb Fibonacciego.
W tym pliku korzystając z wiedzy zdobytej na temat ciągu i jego własności (formuły jawne, postać macierzowa i tożsamości), spróbujemy napisać implementację (kod w języku Python), który będzie wydajnie generował liczby Fibonacciego dla dużych wykładników.
%timeit fibonacci(10000) #szybkie sprawdzenie wydajności wbudowanej w SageMath procedury obliczania liczb Fibonacciego
6.7 µs ± 272 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
def Fib1(n): #implementacja iteracyjna
a=0
b=1
c=0
for i in [1..n]:
c=a+b
a=b
b=c
return a
%timeit Fib1(10000) #wydajność naszej funkcji jest raczej kiepska w porównaniu z wbudowaną!
3.57 ms ± 154 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Formuła Bineta#
Liczby Fibonacciego mogą być wyrażone w terminach formuły Bineta, która podaje postać jawną kolejnych elementów ciągu Fibonacciego. Zaletą tej formuły jest niezależność postaci \(n\)-tej liczby Fibonacciego od jej poprzedników.
\(F_{n}=\frac{1}{\sqrt{5}}\left( \left(\frac{1+\sqrt{5}}{2}\right)^{n}-\left(\frac{1-\sqrt{5}}{2}\right)^{n}\right)\)
R.<x>=QQ[] #inicjalizujemy pierścień wielomianów o współczynnikach wymiernych (QQ) i zmiennej x
K.<phi>=NumberField(x^2-x-1) #inicjalizujemy ciało liczbowe Q(phi) określone przez wielomian x^2-x-1
R
Univariate Polynomial Ring in x over Rational Field
K
Number Field in phi with defining polynomial x^2 - x - 1
phi^2-phi-1
0
print((x^2-x-1).roots(K)) #SageMath może wygenerować listę pierwiastków wielomianu w zadanym ciele
print((x^2-5).roots(K))
print((x^2-7).roots(K)) # używaj print(...) zamiast print ... -> jesteśmy w Python 3.x!
[(phi, 1), (-phi + 1, 1)]
[(2*phi - 1, 1), (-2*phi + 1, 1)]
[]
alpha=phi
beta=1-phi
def BinetFib(n): #implementacja formuły Bineta
return 1/(2*phi-1)*(alpha^n-beta^n)
[BinetFib(n)-Fib1(n) for n in [1..10]] #szybkie sprawdzenie zgodności obu funkcji
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
%timeit BinetFib(10000) #czasowa wydajność nowej procedury jest znacznie lepsza niż poprzedniej formuły iteracyjnej!
52.8 µs ± 2.98 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Numeryczna formuła Bineta#
W formule Bineta wyrażenia symboliczne możemy zamienić na ich numeryczne odpowiedniki. Sprawdzimy jak to wpłynie na wydajność generowania liczb Fibonacciego.
embed1=K.embeddings(RealField(10000)) #lista zanurzeń ciała K w liczby rzeczywiste (są dwa)
embed1[0]((phi*2-1)^2)
5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
1/phi^1000
-43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875*phi + 70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501
golden=embed1[1](phi) #generujemy wartości liczbowe naszych abstrakcyjnych wersji wyrażeń symbolicznych
goldenhat=embed1[0](phi)
sq5=sqrt(RealField(10000)(5))
def BinetNum(n): #numeryczna formuła Bineta
return round(1/sq5*(golden^n-goldenhat^n))
[BinetNum(i)-Fib1(i) for i in [1..20]]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
%timeit BinetNum(10000) #to jest wolniejsze niż poprzednia implementacja!
215 µs ± 10.3 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit BinetFib(100)
12.6 µs ± 449 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit BinetNum(100)
116 µs ± 8.54 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
print goldenhat
print abs(goldenhat/sq5)
print abs(goldenhat^2/sq5)
-0.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137484754088075386891752126633862223536931793180060766726354433389086595939582905638322661319928290267880675208766892501711696207032221043216269548626296313614438149758701220340805887954454749246185695364864449241044320771344947049565846788509874339442212544877066478091588460749988712400765217057517978834166256249407589069704000281210427621771117778053153171410117046665991466979873176135600670874807101317952368942752194843530567830022878569978297783478458782289110976250030269615617002504643382437764861028383126833037242926752631165339247316711121158818638513316203840052221657912866752946549068113171599343235973494985090409476213222981017261070596116456299098162905552085247903524060201727997471753427775927786256194320827505131218156285512224809394712341451702237358057727861600868838295230459264787801788992199027077690389532196819861514378031499741106926088674296226757560523172777520353613936210767389376455606060592165894667595519004005559089502295309423124823552122124154440064703405657347976639723949499465845788730396230903750339938562102423690251386804145779956981224457471780341731264532204163972321340444494873023154176768937521030687378803441700939544096279558986787232095124268935573097045095956844017555198819218020640529055189349475926007348522821010881946445442223188913192946896220023014437702699230078030852611807545192887705021096842493627135925187607778846658361502389134933331223105339232136243192637289106705033992822652635562090297986424727597725655086154875435748264718141451270006023890162077732244994353088999095016803281121943204819643876758633147985719113978153978074761507722117508269458639320456520989698555678141069683728840587461033781054443909436835835813811311689938555769754841491445341509129540700501947754861630754226417293946803673198058618339183285991303960720144559504497792120761247856459161608370594987860069701894098864007644361709334172709191433650137157660114803814306262380514321173481510055901345610118007905063814215270930858809287570345050780814545881990633612982798141174533927312080928972792221329806429468782427487401745055406778757083237310975915117762978443284747908176518097787268416117632503861211291436834376702350371116330725869883258710336322238109809012110198991768414917512331340152733843837234500934786049792945991582201258104598230925528721241370436149102054718554961180876426576511060545881475604431784798584539731286301625448761148520217064404111660766950597757832570395110878230827106478939021115691039276838453863333215658296597731034360323225457436372041244064088826737584339536795931232213437320995749889469956564736007295999839128810319742631251797141432012311279551894778172691415891177991956481255800184550656329528598591000908621802977563789259991649946428193022293552346674759326951654214021091363018194722707890122087287361707348649998156255472811373479871656952748900814438405327483781378246691744422963491470815700735254570708977
0.276393202250021030359082633126872376455938164038847427572910275458947907436219510058558559162121772503049182384922621649574673227555292613641363987846654729113322182680812083418872335467736014341946423864958246621499657660758593555791356746090274740737277112370048259755931838822409109050150762860927027110151791135845731010590086830642298025132111557491024586704381682307850002257519846956588496404233166748750118482186059199943757914475645776444389369365717976590666801706604025364772879865825038579736409526211449561031293886433995424286004340443304308243542177804749993946076876599499071323512447027794323374633392551414649473766932150536657775768236272297336759231989555668417426649410690186377365680131352805301002981918104757355403796547785880776708740180367418889582950419295187959654400505649314444814442748761135834498973756368742897555038121057531709659552528388454427679826232340953908147042439642201560194584461922093560636027697124393700051778614782265140754648487895365444495929277212757846522124708878787881566821066480896199198888182099540938115375035289575575169111987059318868530404672055210100106830842253920753819249932012287579515261949722639170844008603755108505643931653747093559167205535731911101025395369164646212495793862524239311659812091180744088202642553580975146212885380590980808631196488960236156395871894188962130104814798530295435797823610710911555362217361410620755995397112459460153984393829477638490961422458995780631501274572814962478444230668327699522173013333755378932153572751361472542178658993201435469472887581940402715054480454868982769024912850347056371709745998795221967584453551001129382200180996639343775611359036071224648273370402856177204369204385047698455576498346108272135908695802060288864371786063254231882507793243789111218112632832837237737662012288846049031701710931698174091859899610449027673849154716541210639265360388276332163342801739207855971088099100441575847750428708167678325881002427986059621180227198471127658133165458161713269972568467977039237138747523897135765303697988819730877976398418987237156945813828238142485930989843837090823601873277403440371765093214537583814205441555734038714106243514502519650988918644248583352537804816976447404311343050418364696380442546316776473499227757741712633124659529925776733854826023348257932735552378038197577960201646317016497533731969453231232553099813042790041410801683559748379080353814894255751725912770179589056289007763824714684697787890823704879113643040283092053742739674910247770295956587119177667846609880448433485920977824353834578704212195776861792144632309227333356868340680453793127935354908512725591751187182234652483132092640813753557312535800850022106008687052798540800032174237936051473749640571713597537744089621044365461716821764401608703748839963089868734094280281799818275639404487242148001670010714361395541289530665048134609669157195781727396361055458421975582542527658530270000368748905437725304025668609450219837112318934503243724350661651115407301705836859852949085858205
0.170820393249936908922752100619382870632185507883457717281269173623156277691341469824324322513634682490852452845232135051275980317334122159075908036460035812660033451957563749743382993596791956974160728405125260135501027017724219332625929761729175777788168662889855220732204483532772672849547711417218918669544626592462806968229739508073105924603665327526926239886854953076449993227440459130234510787300499753749644553441822400168726256573062670666831891902846070227999594880187923905681360402524884260790771421365651316906118340698013727141986978670087075269373466585750018161769370201502786029462658916617029876099822345756051578699203548390026672695291183107989722304031332994747720051767929440867902959605941584096991054245685727933788610356642357669873779458897743331251148742114436121036798483052056665556671753716592496503078730893771307334885636827404871021342414834636716960521302977138275558872681073395319416246614233719318091916908626818899844664155653204577736054536313903666512212168361726460433625873363636355299536800557311402403335453701377185653874894131273274492664038822043394408785983834369699679507473238237738542250203963137261454214150832082487467974188734674483068205038758719322498383392804266696923813892506061362512618412427282065020563726457767735392072339257074561361343858227057574106410533119291530812384317433113609685555604409113692606529167867265333913347915768137732013808662621619538046818511567084527115732623012658105496176281555112564667307995016901433480959998733863203539281745915582373464023020395693591581337254178791854836558635393051692925261448958830884870762003614334097246639346996611853399457010081968673165922891786326055179888791431468386892386844856904633270504961675183592273912593819133406884641810237304352476620268632666345662101501488286787013963133461852904894867204905477724420301168652916978452535850376368082203918835171003509971594782376432086735702698675272456748713875496965022356992716041821136459318404586617025600503625514860190082294596068882288583757428308592704088906033540807366070804743038288529162558515285572542207030468488727529194380167789678884704720356387248557383675332797883857681269456492441047033244067254249942386585549070657787065970848744905910858672361049670579502316726774862100626021410222669798435521929955226201793342865885407266119395061048950507398804091640306302340700560871629875767594949320754862758938555317232744822261689461232831132976708525855945906636327528885362659070879150723838771780975269256689112130238642466996460170358654699542237066526938496263887363412669414623566103072317999929394977958638620616193935274461823224746438453296042550603722077558739328062392597449933681973938841604377599903477286191845578751078284859207386767731136866903614849534706795173888753480110730393797717159154600545173081786538273555994989967856915813376131408004855596170992528412654817810916833624734073252372417024409189998893753283686824087922994171649340488663043196490268826948015046653778094882489420441152742425386
def BinetNum2(n): # możemy w drugim podejściu wykorzystać fakt, że drugi człon formuły z phihat jest co do moduły mniejszy od 1
if (n <=1):
return n
else:
if (n % 2 == 0):
return floor(1/sq5*(golden^n))
else:
return ceil(1/sq5*(golden^n))
[i for i in [1..7200] if not (BinetNum2(i)==BinetFib(i))]
[7195, 7197, 7199]
%timeit BinetNum2(10000) #to jest szybsze niż BinetNum1, ale nadal wolniejsze niż wersja symboliczna
120 µs ± 7.34 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
def BinetNum3(n,prec): #implementacja pozwalająca zarządzać precyzją liczb
sq5=sqrt(RealField(prec)(5))
golden=(1+sq5)/2
goldenhat=(1-sq5)/2
if (n <=1):
return n
else:
if (n % 2 == 0):
return floor(1/sq5*(golden^n))
else:
return ceil(1/sq5*(golden^n))
%timeit BinetNum3(100,100)
The slowest run took 7.45 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 10.3 µs per loop
%timeit fibonacci(100)
The slowest run took 21.45 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 1.63 µs per loop
BinetNum3(100,100)-fibonacci(100)
0
%timeit BinetFib(100)
The slowest run took 7.45 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 25.2 µs per loop
%timeit BinetNum3(10000,7700)#z odpowiednio dobraną precyzją możemy uzyskać implementację szybszą niż symboliczna
#ale z drugiej strony potrzebujemy dostatecznej precyzji, aby wyniki były poprawne
1000 loops, best of 3: 224 µs per loop
%timeit fibonacci(10000)
The slowest run took 4.95 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 11.5 µs per loop
BinetNum3(10000,7700)-fibonacci(10000)
0
Czy istnieje szybsza implementacja niż BinetNum?#
m=matrix(2,2,[1,1,1,0]) #inicjalizujemy macierz o wymiarze dwa i współczynnikach 1,1,1,0
print m
[1 1]
[1 0]
[m^n-matrix(2,2,[fibonacci(n+1), fibonacci(n),fibonacci(n),fibonacci(n-1)]) for n in [0..10]] #ciekawa tożsamość
[
[0 0] [0 0] [0 0] [0 0] [0 0] [0 0] [0 0] [0 0] [0 0] [0 0]
[0 0], [0 0], [0 0], [0 0], [0 0], [0 0], [0 0], [0 0], [0 0], [0 0],
[0 0]
[0 0]
]
Ćwiczenie
Udowodnij równość
\(\left(\begin{array}{cc} 1 & 1 \\ 1 & 0 \end{array}\right)^n = \left(\begin{array}{cc} F_{n+1} & F_{n} \\ F_{n} & F_{n-1} \end{array}\right)\)
Szybkie potęgowanie#
def pow_sq(x, n): #szybkie potęgowanie
if n < 0:
x = x^(-1)
n = -n
if n == 0:
return x*(x^(-1))
y = matrix(2,2,[1,0,0,1])
while n > 1:
if n % 2 == 0:
x = x * x
n = n / 2;
else:
y = x * y;
x = x * x;
n = (n - 1) / 2
return x * y
m=matrix(2,2,[1,1,1,0]) #liczby Fibonacciego obliczane metodą macierzową
def FibMat(n):
return pow_sq(m,n)[0][1]
%timeit FibMat(10000)
128 µs ± 7.7 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit BinetNum3(10000,10000)
1000 loops, best of 3: 358 µs per loop
%timeit fibonacci(10000)
7.66 µs ± 760 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit BinetFib(10000)
10000 loops, best of 3: 107 µs per loop
%timeit FibMat(50000)
1000 loops, best of 3: 751 µs per loop
%timeit fibonacci(50000)
10000 loops, best of 3: 109 µs per loop
%timeit BinetFib(50000) #dla dostatecznie dużych indeksów nasza implementacja macierzowa działa lepiej niż formuła Bineta
1000 loops, best of 3: 577 µs per loop
Aby dalej przyspieszyć działanie naszej procedury możemy wykorzystać dwie znane nam tożsamości:
\(F(2k) = F(k)(2F(k+1)-F(k))\)
\(F(2k+1)=F(k+1)^2+F(k)^2\)
matgen=matrix(2,2,[1,1,1,0]) #poprawiona wersja implementacji macierzowej
def FibMatv2(n):
if n % 2 == 0:
m= n/2
S=pow_sq(matgen,m)
Fk=S[0][1]
FSk=S[0][0]
return Fk*(2*FSk-Fk)
else:
m=(n-1)/2
S=pow_sq(matgen,m)
Fk=S[0][1]
FSk=S[0][0]
return Fk^2+FSk^2
%timeit FibMat(100000)
1.13 ms ± 77.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit FibMatv2(100000)
573 µs ± 61.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit fibonacci(100000)
159 µs ± 6.73 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit BinetFib(100000)
744 µs ± 52.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
import timeit as TI
#def FibMatTime(n,iter1):
# return TI.timeit('FibMat(%d)'%(n),"from __main__ import FibMat",number=iter1)/iter1
def FunTime(n,iter1,funname):
return TI.timeit('%s(%d)'%(funname,n),"from __main__ import %s"%(funname),number=iter1)/iter1
li1=list_plot([(i,FunTime(1000*i,100,"FibMat")) for i in [1..100]],color='red');
li2=list_plot([(i,FunTime(1000*i,100,"fibonacci")) for i in [1..100]],color='green');
li3=list_plot([(i,FunTime(1000*i,100,"FibMatv2")) for i in [1..100]],color='blue');
li4=list_plot([(i,FunTime(1000*i,100,"BinetFib")) for i in [1..100]],color='purple');
li1+li2+li3+li4
