from gmpy2 import mpz
#
# Fast doubling Fibonacci algorithm (Python)
#
# Copyright (c) 2015 Project Nayuki
# All rights reserved. Contact Nayuki for licensing.
# https://www.nayuki.io/page/fast-fibonacci-algorithms
#
# (Public) Returns F(n).
def fib_nayuki(n):
# (Private) Returns the tuple (F(n), F(n+1)).
def _fib(n):
if n == 0:
return mpz(0), mpz(1)
a, b = _fib(n >> 1)
c = a * (b * 2 - a)
d = a * a + b * b
if (n & 1):
return d, c + d
return c, d
a, b = _fib(n >> 1)
if (n & 1):
return a * a + b * b
return a * (b * 2 - a)
def fib_primo1(n):
def fib_inner(n):
'''Returns F[n] and L[n]'''
if n == 0:
return mpz(0), mpz(2)
u, v = fib_inner(n >> 1)
q = (n & 2) - 1
u, v = u * v, v * v + 2*q
if (n & 1):
u1 = (u + v) >> 1
return u1, 2*u + u1
return u, v
u, v = fib_inner(n >> 1)
if (n & 1):
q = (n & 2) - 1
u1 = (u + v) >> 1
return v * u1 + q
return u * v
def fib_takahashi(n):
def fib_inner(n):
'''Returns F[n] and L[n]'''
if n == 0:
return mpz(0), mpz(2)
u, v = fib_inner(n >> 1)
q = (n & 2) - 1
t = u * u
# F[m + 1]
u1 = (u + v) >> 1
u = 2*(u1 * u1 + q) - 3*t
v = 5*t - 2*q
if (n & 1):
u1 = (u + v) >> 1
return u1, u1 + 2*u
return u, v
u, v = fib_inner(n >> 1)
if (n & 1):
q = (n & 2) - 1
u1 = (u + v) >> 1
return v * u1 + q
return u * v
def fib_primo2(n):
def fib_inner(n):
'''Returns F[n] and F[n+1]'''
if n == 0:
return mpz(0), mpz(1)
u, v = fib_inner(n >> 1)
q = (n & 2) - 1
u *= u
v *= v
if (n & 1):
return u + v, 3*v - 2*(u - q)
return 2*(v + q) - 3*u, u + v
u, v = fib_inner(n >> 1)
# L[m]
l = 2*v - u
if (n & 1):
q = (n & 2) - 1
return v * l + q
return u * l
from timeit import repeat
a, b = 0, 1
for n in range(1000):
assert(a == fib_nayuki(n))
assert(a == fib_primo1(n))
assert(a == fib_takahashi(n))
assert(a == fib_primo2(n))
a, b = b, a + b
runs = 5
for test in ['fib_primo2']:
print test, 'small values:',
print '%f'%min(repeat(
setup = 'from __main__ import %s'%test,
stmt = 'for n in range(25000): %s(n)'%test,
number = 1, repeat = runs)), 'secs per run'
print test, 'large values:',
print '%f'%(min(repeat(
setup = "from __main__ import %s\n"%test +
"from random import seed, randint\n" +
"seed('fibonacci')",
stmt = '%s(randint(50000000, 100000000))'%test,
number = runs, repeat = 3))/runs), 'secs per value'
print