#!/usr/bin/python3 # # Written by Marek Behun # Oct 23, 2015 # Public domain # φ = (1+√5)/2 # ψ = (1-√5)/2 # fib(n) = (φⁿ-ψⁿ) / √5 def fib1(n): # Computation directly from definition # This is inaccurate for large n, because we are using floating point # arithmetic, which has limited precision from math import sqrt φ = (1 + sqrt(5)) / 2 ψ = (1 - sqrt(5)) / 2 return int((φ**n - ψ**n) / sqrt(5)) # Let's try this by representing numbers of form a+b√5 as pairs (a, b) # and define multiplication, substraction and integer exponentiation # of such numbers def mul(x, y): # x = a+b√5 # y = c+d√5 # x·y = (a+b√5)·(c+d√5) = ac+ad√5+bc√5+bd√5√5 = (ac+5bd)+(ad+bc)√5 (a, b) = x # x = a+b√5 (c, d) = y # y = c+d√5 return (a*c+5*b*d, a*d+b*c) def pow(x, n): # xⁿ = (a+b√5)ⁿ # We use exponentiation by squaring if n == 0: # x⁰ = 1 = 1+0√5 return (1, 0) elif n == 1: # x¹ = x return x else: if n % 2 == 0: # for even n, xⁿ = (x²)ᵐ, m = n/2 return pow(mul(x, x), n // 2) else: # for odd n, xⁿ = x·xⁿ¯¹ = x·(x²)ᵐ, m = (n-1)/2 # for odd n, (n-1) // 2 = n // 2 return mul(x, pow(mul(x, x), n // 2)) # But we can do exponentiation by squaring without recursion! def pow(x, n): if n == 0: return (1, 0) y = (1, 0) while n > 1: if n % 2 == 0: x = mul(x, x) n //= 2 else: y = mul(x, y) x = mul(x, x) n //= 2 return mul(x, y) def fib2(n): # Theorem: If (a+b√5)ⁿ = c+d√5, then (a-b√5)ⁿ = c-d√5. # Proof (by induction): # For n = 0 and n = 1 the statement is true. # For n ← n-1 we have: # We know that (a+b√5)ⁿ¯¹ = c+d√5, then (a-b√5)ⁿ¯¹ = c-d√5. # Now, # (a+b√5)ⁿ = (c+d√5)·(a+b√5) = (ac+5bd)+(ad+bc)√5 # and # (a-b√5)ⁿ = (c-d√5)·(a-b√5) = (ac+5bd)-(ad+bc)√5. # Q.E.D. # fib(n) = (φⁿ - ψⁿ) / √5 = ((2φ)ⁿ - (2ψ)ⁿ) / (2ⁿ√5) # But using the theorem above for # 2φ = 1+√5 # 2ψ = 1-√5 # we know that if (2φ)ⁿ = a+b√5, then (2ψ)ⁿ = a-b√5, and so # fib(n) = ((2φ)ⁿ - (2ψ)ⁿ) / (2ⁿ√5) = 2b√5 / (2ⁿ√5) = b/2ⁿ¯¹ # Therefore it is sufficient to compute only (2φ)ⁿ and take # the coefficient of √5: φ2 = (1, 1) (a, b) = pow(φ2, n) return 2*b // 2**n def fib3(n): # Inlining the exponentiation and removing code repetition we get if n == 0: return 0 x = (1, 1) y = (1, 0) _n = n while _n > 1: if _n % 2 == 1: y = mul(x, y) x = mul(x, x) _n //= 2 a, b = mul(x, y) return b // 2**(n-1) def fib4(n): # We see that at the last step in fib3 we divide by 2ⁿ¯¹. # That means that there are n-1 redundant bits in the result # before we remove them. Those bits were added during the # execution of the while cycle, 2ⁱ bits in the x variable in # the i-th iteration starting from 1 and 1 bit in the y # variable everytime y is changed. Why? # # Both x₊ and y₊ at the end are some integer powers of x. # If x₀ and y₀ are x and y at the beginning and x₊ and y₊ at # the end, both x₊ and y₊ are some integer powers of x₀: # x₊ = x₀ᵏ and y₊ = x₀ᵐ # # # It would save space in the variables if we get rid of these # bits during the cycle. By the theorems below, we can do this by # halving x at the end of each iteration and by halving y at the # end of each iteration in which y is changed but the first # (after the first iteration y has odd coefficients), and # then halve y after the cycle if it was changed at least once. # # Theorem: For every integer n ≥ 1 there are integers a and b # with the same parity such that (1+√5)ⁿ = 2ⁿ¯¹·(a+b√5). # Proof (by induction): # For n = 1 obviously holds (a = 1, b = 1, both odd). # For n ← n-1 we have: # (1+√5)ⁿ¯¹ = 2ⁿ¯²·(a+b√5) # (1+√5)ⁿ = 2ⁿ¯²·(a+b√5)·(1+√5) = 2ⁿ¯²·((a+5b)+(a+b)√5) # We know that a and b have the same parity. # If they are both odd, there are some integers c and d # such that a = 2c-1 and b = 2d-1. Then # a+5b = 2c-1+5·(2d-1) = 2·(c+5d-3) # a+b = 2c-1+2d-1 = 2·(c+d-1) # and # (1+√5)ⁿ = 2ⁿ¯¹·((c+5d-3)+(c+d-1)√5). # Here c+5d-3 and c+d-1 have the same parity, because # their difference, 4d-2, is even. # If they are both even, there are some integers c and d # such that a = 2c and b = 2d. Then # a+5b = 2c+10d = 2·(c+5d) # a+b = 2c+2d = 2·(c+d) # and # (1+√5)ⁿ = 2ⁿ¯¹·((c+5d)+(c+d)√5). # Here c+5d and c+d have the same parity as well, because # their difference, 4d, is even. # Q.E.D. # # Theorem: If a and b are odd integers, then there are odd # integers c and d such that (a+b√5)²/2 = c+d√5. # Proof: # We have (a+b√5)² = (a²+5b²)+(2ab)√5, thus # c = (a²+5b²)/2 and d = ab. # Oddity of c: # Let a = 2e-1 and b = 2f-1, then # c = (a²+5b²)/2 = ((4c²-4c+1)+5·(4d²-4d+1))/2 = # = (4c²-4c+20d²-20d+6)/2 = 2c²-2c+20d²-20d+3 = # = 2·(c²-c+10d²-10d)+3 # thus c is odd. # Oddity of d: # d is odd because it is a product of odd numbers. # Q.E.D. # # Theorem: If A, B, C and D are odd integers, then there are # integers E and F such that (A+B√5)·(C+D√5)/2 = E+F√5. # Proof: # Let A = 2a-1, B = 2b-1, C = 2c-1 and D = 2d-1 for some # integers a, b, c, d. Then we have # (A+B√5)·(C+D√5) = ((2a-1) + (2b-1)√5) · ((2c-1) + (2d-1)√5) = # = ((2a-1)·(2c-1) + 5·(2b-1)·(2d-1)) + ((2a-1)·(2d-1) + (2b-1)·(2c-1))√5 = # = ((4ac-2a-2c+1) + 5·(4bd-2b-2d+1)) + ((4ad-2a-2d+1) + (4bc-2b-2c+1))√5 = # = (4ac-2a-2c+20bd-10b-10d+6) + (4ad-2a-2d+4bc-2b-2c+2)√5 = # = 2·(2ac-a-c+10bd-5b-5d+3) + 2·(2ad-a-d+2bc-b-c+1)√5 # = 2·(E+F)√5 # Q.E.D. if n == 0: return 0 x = (1, 1) y = (1, 0) j = 0 while n > 1: if n % 2 == 1: a, b = mul(x, y) if j == 1: a, b = a // 2, b // 2 j = 1 y = a, b a, b = mul(x, x) x = a // 2, b // 2 n //= 2 a, b = mul(x, y) return b // 2**j def fib5(n): # After inlining multiplication if n == 0: return 0 a, b, c, d, j = 1, 1, 1, 0, 0 while n > 1: if n % 2 == 1: # x·y = (a+b√5)·(c+d√5) = (ac+5bd)+(ad+bc)√5 c, d = a*c+5*b*d, a*d+b*c if j == 1: c, d = c // 2, d // 2 j = 1 # x·x/2 = (a+b√5)·(a+b√5)/2 = (a²+5b²)/2 + ab√5 a, b = (a*a+5*b*b) // 2, a*b n //= 2 # x·y = (a+b√5)·(c+d√5) = (ac+5bd)+(ad+bc)√5 return (a*d+b*c) // 2**j # Finally use bitwise operators instead of arithmetic def fib(n): "Returns the n-th Fibonacci number" if n == 0: return 0 a, b, c, d, j = 1, 1, 1, 0, 0 while n > 1: if n & 1 == 1: c, d, j = a*c+5*b*d >> j, a*d+b*c >> j, 1 a, b, n = a*a+5*b*b >> 1, a*b, n >> 1 return a*d+b*c >> j if __name__ == '__main__': for i in range(0,10): print(fib(i))