Fibonacci matrix-exponentiation

Fibonacci matrix-exponentiation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The Fibonacci sequence defined with matrix-exponentiation:

      F0 = 0 
      F1 = 1 

${\begin{pmatrix}F_{n+1}\\F_{n}\end{pmatrix}}={\begin{pmatrix}1&1\\1&0\end{pmatrix}}^{n}{\begin{pmatrix}F_{1}\\F_{0}\end{pmatrix}}$

Task

Write a function using matrix-exponentiation to generate Fibonacci number for n = 2³² and n = 2⁶⁴.

Display only the 20 first digits and 20 last digits of each Fibonacci number.

Extra

Generate the same Fibonacci numbers using Binet's formula below or another method.

$F_{n}={\cfrac {\varphi ^{n}-(-\varphi )^{-n}}{\sqrt {5}}}$

$\varphi ={\cfrac {1+{\sqrt {5}}}{2}}$

Related tasks

•   Fibonacci sequence
•   Matrix-exponentiation operator

Julia

Because Julia uses the GMP library for its BigInt type, a BigInt cannot be larger in bits than about 2^37. This prevents generation of the 2^64-th fibonacchi number, due to BigInt overflow. The Binet method actually overflows even with the 2^32-rd fibonacchi number, so the Lucas method is used as the alternative method.

<lang julia>const b = [big"1" 1; 1 0] matrixfibonacci(n) = n == 0 ? 0 : n == 1 ? 1 : (b^(n+1))[2,2]

1. This is derived from the Schönhage-Strassen algorithm using Lucas numbers

function lucasfibonacci(n)

   function inner(n)
       if n == 0
           return zero(n), one(n)
       end
       u, v = inner(n >> 1)
       q = (n & 2) - 1
       u *= u
       v *= v
       if isodd(n)
           return u + v, 3 * v - 2 * (u - q)
       end
       return BigInt(2 * (v + q) - 3 * u), BigInt(u + v)
   end
   u, v = inner(n >> 1)
   l = 2*v - u # the lucas function
   if isodd(n)
       q = (n & 2) - 1
       return v * l + q
   end
   return u * l

end

m2s(bits) = string(matrixfibonacci(big"2"^bits)) l2s(bits) = string(lucasfibonacci(big"2"^bits)) firstlast(s) = (length(s) < 44 ? s : s[1:20] * "..." * s[end-20+1:end])

println("N", " "^23, "Matrix", " "^40, "Lucas\n", "-"^97) println("2^32 ", rpad(firstlast(m2s(32)), 45), rpad(firstlast(l2s(32)), 45))

1. println("2^64 ", rpad(firstlast(m2s(44)), 45), rpad(firstlast(l2s(44)), 45))

</lang>

Output:

N                       Matrix                                        Lucas
-------------------------------------------------------------------------------------------------
2^32  61999319689381859818...39623735538208076347  61999319689381859818...39623735538208076347