1 #!/usr/bin/python3
    2 #
    3 # Written by Marek Behun
    4 # Oct 23, 2015
    5 # Public domain
    6 
    7 # φ = (1+√5)/2
    8 # ψ = (1-√5)/2
    9 # fib(n) = (φⁿ-ψⁿ) / √5
   10 
   11 def fib1(n):
   12         # Computation directly from definition
   13         # This is inaccurate for large n, because we are using floating point
   14         # arithmetic, which has limited precision
   15         from math import sqrt
   16         φ = (1 + sqrt(5)) / 2
   17         ψ = (1 - sqrt(5)) / 2
   18         return int((φ**n - ψ**n) / sqrt(5))
   19 
   20 # Let's try this by representing numbers of form a+b√5 as pairs (a, b)
   21 # and define multiplication, substraction and integer exponentiation
   22 # of such numbers
   23 
   24 def mul(x, y):
   25         # x = a+b√5
   26         # y = c+d√5
   27         # x·y = (a+b√5)·(c+d√5) = ac+ad√5+bc√5+bd√5√5 = (ac+5bd)+(ad+bc)√5
   28         (a, b) = x # x = a+b√5
   29         (c, d) = y # y = c+d√5
   30         return (a*c+5*b*d, a*d+b*c)
   31 
   32 def pow(x, n):
   33         # xⁿ = (a+b√5)ⁿ
   34         # We use exponentiation by squaring
   35         if n == 0:
   36                 # x⁰ = 1 = 1+0√5
   37                 return (1, 0)
   38         elif n == 1:
   39                 # x¹ = x
   40                 return x
   41         else:
   42                 if n % 2 == 0:
   43                         # for even n, xⁿ = (x²)ᵐ, m = n/2
   44                         return pow(mul(x, x), n // 2)
   45                 else:
   46                         # for odd n, xⁿ = x·xⁿ¯¹ = x·(x²)ᵐ, m = (n-1)/2
   47                         # for odd n, (n-1) // 2 = n // 2
   48                         return mul(x, pow(mul(x, x), n // 2))
   49 
   50 # But we can do exponentiation by squaring without recursion!
   51 def pow(x, n):
   52         if n == 0:
   53                 return (1, 0)
   54         y = (1, 0)
   55         while n > 1:
   56                 if n % 2 == 0:
   57                         x = mul(x, x)
   58                         n //= 2
   59                 else:
   60                         y = mul(x, y)
   61                         x = mul(x, x)
   62                         n //= 2
   63         return mul(x, y)
   64 
   65 def fib2(n):
   66         # Theorem: If (a+b√5)ⁿ = c+d√5, then (a-b√5)ⁿ = c-d√5.
   67         # Proof (by induction):
   68         #    For n = 0 and n = 1 the statement is true.
   69         #    For n ← n-1 we have:
   70         #        We know that (a+b√5)ⁿ¯¹ = c+d√5, then (a-b√5)ⁿ¯¹ = c-d√5.
   71         #        Now,
   72         #            (a+b√5)ⁿ = (c+d√5)·(a+b√5) = (ac+5bd)+(ad+bc)√5
   73         #        and
   74         #            (a-b√5)ⁿ = (c-d√5)·(a-b√5) = (ac+5bd)-(ad+bc)√5.
   75         #    Q.E.D.
   76 
   77         # fib(n) = (φⁿ - ψⁿ) / √5 = ((2φ)ⁿ - (2ψ)ⁿ) / (2ⁿ√5)
   78         # But using the theorem above for
   79         #    2φ = 1+√5
   80         #    2ψ = 1-√5
   81         # we know that if (2φ)ⁿ = a+b√5, then (2ψ)ⁿ = a-b√5, and so
   82         # fib(n) = ((2φ)ⁿ - (2ψ)ⁿ) / (2ⁿ√5) = 2b√5 / (2ⁿ√5) = b/2ⁿ¯¹
   83 
   84         # Therefore it is sufficient to compute only (2φ)ⁿ and take
   85         # the coefficient of √5:
   86         φ2 = (1, 1)
   87         (a, b) = pow(φ2, n)
   88 
   89         return 2*b // 2**n
   90 
   91 def fib3(n):
   92         # Inlining the exponentiation and removing code repetition we get
   93         if n == 0:
   94                 return 0
   95         x = (1, 1)
   96         y = (1, 0)
   97         _n = n
   98         while _n > 1:
   99                 if _n % 2 == 1:
  100                         y = mul(x, y)
  101                 x = mul(x, x)
  102                 _n //= 2
  103         a, b = mul(x, y)
  104         return b // 2**(n-1)
  105 
  106 def fib4(n):
  107         # We see that at the last step in fib3 we divide by 2ⁿ¯¹.
  108         # That means that there are n-1 redundant bits in the result
  109         # before we remove them. Those bits were added during the
  110         # execution of the while cycle, 2ⁱ bits in the x variable in
  111         # the i-th iteration starting from 1 and 1 bit in the y
  112         # variable everytime y is changed. Why?
  113         #
  114         # Both x₊ and y₊ at the end are some integer powers of x.
  115         # If x₀ and y₀ are x and y at the beginning and x₊ and y₊ at
  116         # the end, both x₊ and y₊ are some integer powers of x₀:
  117         #    x₊ = x₀ᵏ   and   y₊ = x₀ᵐ
  118         #
  119         #
  120         # It would save space in the variables if we get rid of these
  121         # bits during the cycle. By the theorems below, we can do this by
  122         # halving x at the end of each iteration and by halving y at the
  123         # end of each iteration in which y is changed but the first
  124         # (after the first iteration y has odd coefficients), and
  125         # then halve y after the cycle if it was changed at least once.
  126         #
  127         # Theorem: For every integer n ≥ 1 there are integers a and b
  128         #    with the same parity such that (1+√5)ⁿ = 2ⁿ¯¹·(a+b√5).
  129         # Proof (by induction):
  130         #    For n = 1 obviously holds (a = 1, b = 1, both odd).
  131         #    For n ← n-1 we have:
  132         #        (1+√5)ⁿ¯¹ = 2ⁿ¯²·(a+b√5)
  133         #        (1+√5)ⁿ = 2ⁿ¯²·(a+b√5)·(1+√5) = 2ⁿ¯²·((a+5b)+(a+b)√5)
  134         #        We know that a and b have the same parity.
  135         #        If they are both odd, there are some integers c and d
  136         #        such that a = 2c-1 and b = 2d-1. Then
  137         #            a+5b = 2c-1+5·(2d-1) = 2·(c+5d-3)
  138         #            a+b  = 2c-1+2d-1 = 2·(c+d-1)
  139         #        and
  140         #            (1+√5)ⁿ = 2ⁿ¯¹·((c+5d-3)+(c+d-1)√5).
  141         #        Here c+5d-3 and c+d-1 have the same parity, because
  142         #        their difference, 4d-2, is even.
  143         #        If they are both even, there are some integers c and d
  144         #        such that a = 2c and b = 2d. Then
  145         #            a+5b = 2c+10d = 2·(c+5d)
  146         #            a+b  = 2c+2d = 2·(c+d)
  147         #        and
  148         #            (1+√5)ⁿ = 2ⁿ¯¹·((c+5d)+(c+d)√5).
  149         #        Here c+5d and c+d have the same parity as well, because
  150         #        their difference, 4d, is even.
  151         #    Q.E.D.
  152         #
  153         # Theorem: If a and b are odd integers, then there are odd
  154         #    integers c and d such that (a+b√5)²/2 = c+d√5.
  155         # Proof:
  156         #    We have (a+b√5)² = (a²+5b²)+(2ab)√5, thus
  157         #        c = (a²+5b²)/2   and   d = ab.
  158         #    Oddity of c:
  159         #        Let a = 2e-1 and b = 2f-1, then
  160         #            c = (a²+5b²)/2 = ((4c²-4c+1)+5·(4d²-4d+1))/2 =
  161         #              = (4c²-4c+20d²-20d+6)/2 = 2c²-2c+20d²-20d+3 =
  162         #              = 2·(c²-c+10d²-10d)+3
  163         #        thus c is odd.
  164         #    Oddity of d:
  165         #        d is odd because it is a product of odd numbers.
  166         #    Q.E.D.
  167         #
  168         # Theorem: If A, B, C and D are odd integers, then there are
  169         #    integers E and F such that (A+B√5)·(C+D√5)/2 = E+F√5.
  170         # Proof:
  171         #    Let A = 2a-1, B = 2b-1, C = 2c-1 and D = 2d-1 for some
  172         #    integers a, b, c, d. Then we have
  173         #        (A+B√5)·(C+D√5) = ((2a-1) + (2b-1)√5) · ((2c-1) + (2d-1)√5) =
  174         #        = ((2a-1)·(2c-1) + 5·(2b-1)·(2d-1)) + ((2a-1)·(2d-1) + (2b-1)·(2c-1))√5 =
  175         #        = ((4ac-2a-2c+1) + 5·(4bd-2b-2d+1)) + ((4ad-2a-2d+1) + (4bc-2b-2c+1))√5 =
  176         #        =     (4ac-2a-2c+20bd-10b-10d+6)    +     (4ad-2a-2d+4bc-2b-2c+2)√5 =
  177         #        =      2·(2ac-a-c+10bd-5b-5d+3)     +      2·(2ad-a-d+2bc-b-c+1)√5
  178         #        = 2·(E+F)√5
  179         #    Q.E.D.
  180 
  181         if n == 0:
  182                 return 0
  183         x = (1, 1)
  184         y = (1, 0)
  185         j = 0
  186         while n > 1:
  187                 if n % 2 == 1:
  188                         a, b = mul(x, y)
  189                         if j == 1:
  190                                 a, b = a // 2, b // 2
  191                         j = 1
  192                         y = a, b
  193                 a, b = mul(x, x)
  194                 x = a // 2, b // 2
  195                 n //= 2
  196         a, b = mul(x, y)
  197         return b // 2**j
  198 
  199 def fib5(n):
  200         # After inlining multiplication
  201         if n == 0:
  202                 return 0
  203         a, b, c, d, j = 1, 1, 1, 0, 0
  204         while n > 1:
  205                 if n % 2 == 1:
  206                         # x·y = (a+b√5)·(c+d√5) = (ac+5bd)+(ad+bc)√5
  207                         c, d = a*c+5*b*d, a*d+b*c
  208                         if j == 1:
  209                                 c, d = c // 2, d // 2
  210                         j = 1
  211                 # x·x/2 = (a+b√5)·(a+b√5)/2 = (a²+5b²)/2 + ab√5
  212                 a, b = (a*a+5*b*b) // 2, a*b
  213                 n //= 2
  214         # x·y = (a+b√5)·(c+d√5) = (ac+5bd)+(ad+bc)√5
  215         return (a*d+b*c) // 2**j
  216 
  217 # Finally use bitwise operators instead of arithmetic
  218 def fib(n):
  219         "Returns the n-th Fibonacci number"
  220         if n == 0:
  221                 return 0
  222         a, b, c, d, j = 1, 1, 1, 0, 0
  223         while n > 1:
  224                 if n & 1 == 1:
  225                         c, d, j = a*c+5*b*d >> j, a*d+b*c >> j, 1
  226                 a, b, n = a*a+5*b*b >> 1, a*b, n >> 1
  227         return a*d+b*c >> j
  228 
  229 if __name__ == '__main__':
  230         for i in range(0,10):
  231                 print(fib(i))