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

Go

Library: GMP(Go wrapper)

This uses matrix exponentiation to calculate the (2^32)nd Fibonacci number which has more than 897 million digits! To improve performance, I've used a GMP wrapper rather than Go's native 'big.Int' type.

I have not attempted to calculate the (2^64)nd Fibonacci number which appears to be well out of reach. <lang go>package main

import (

   "fmt"
   big "github.com/ncw/gmp"
   "time"

)

type vector = []*big.Int type matrix []vector

var (

   zero = new(big.Int)
   one  = big.NewInt(1)

)

func (m1 matrix) mul(m2 matrix) matrix {

   rows1, cols1 := len(m1), len(m1[0])
   rows2, cols2 := len(m2), len(m2[0])
   if cols1 != rows2 {
       panic("Matrices cannot be multiplied.")
   }
   result := make(matrix, rows1)
   temp := new(big.Int)
   for i := 0; i < rows1; i++ {
       result[i] = make(vector, cols2)
       for j := 0; j < cols2; j++ {
           result[i][j] = new(big.Int)
           for k := 0; k < rows2; k++ {
               temp.Mul(m1[i][k], m2[k][j])
               result[i][j].Add(result[i][j], temp)
           }
       }
   }
   return result

}

func identityMatrix(n uint64) matrix {

   if n < 1 {
       panic("Size of identity matrix can't be less than 1")
   }
   ident := make(matrix, n)
   for i := uint64(0); i < n; i++ {
       ident[i] = make(vector, n)
       for j := uint64(0); j < n; j++ {
           ident[i][j] = new(big.Int)
           if i == j {
               ident[i][j].Set(one)
           }
       }
   }
   return ident

}

func (m matrix) pow(n *big.Int) matrix {

   le := len(m)
   if le != len(m[0]) {
       panic("Not a square matrix")
   }
   switch {
   case n.Cmp(zero) == -1:
       panic("Negative exponents not supported")
   case n.Cmp(zero) == 0:
       return identityMatrix(uint64(le))
   case n.Cmp(one) == 0:
       return m
   }
   pow := identityMatrix(uint64(le))
   base := m
   e := new(big.Int).Set(n)
   temp := new(big.Int)
   for e.Cmp(zero) > 0 {
       temp.And(e, one)
       if temp.Cmp(one) == 0 {
           pow = pow.mul(base)
       }
       e.Rsh(e, 1)
       base = base.mul(base)
   }
   return pow

}

func fibonacci(n *big.Int) *big.Int {

   if n.Cmp(zero) == 0 {
       return zero
   }
   m := matrix{{one, one}, {one, zero}}
   m = m.pow(n.Sub(n, one))
   return m[0][0]

}

func main() {

   start := time.Now()
   n := new(big.Int)
   n.Lsh(one, 32)
   s := fibonacci(n).String()
   fmt.Printf("The digits of the 2^32nd Fibonacci number (%d) are:\n", len(s))
   fmt.Printf("  First twenty : %s\n", s[0:20])
   fmt.Printf("  Final twenty : %s\n", s[len(s)-20:])
   fmt.Printf("\nTook %s\n", time.Since(start))

}</lang>

Output:

Timings are for an Intel Core i7 8565U machine, using Go 1.13.7 on Ubuntu 18.04:

The digits of the 2^32nd Fibonacci number (897595080) are:
  First twenty : 61999319689381859818
  Final twenty : 39623735538208076347

Took 11m5.991187993s

Although Go supports big.Float, the precision needed to calculate the (2^32nd) Fibonacci number makes the use of Binet's formula impractical. I have therefore used the same method as the Julia entry for my alternative approach which is more than twice as quick as the matrix exponentiation method.

Translation of: Julia

<lang go>package main

import (

   "fmt"
   big "github.com/ncw/gmp"
   "time"

)

var (

   zero  = new(big.Int)
   one   = big.NewInt(1)
   two   = big.NewInt(2)
   three = big.NewInt(3)

)

func lucas(n *big.Int) *big.Int {

   var inner func(n *big.Int) (*big.Int, *big.Int)
   inner = func(n *big.Int) (*big.Int, *big.Int) {
       if n.Cmp(zero) == 0 {
           return zero, one
       }
       t, q, r := new(big.Int), new(big.Int), new(big.Int)
       u, v := inner(t.Rsh(n, 1))
       t.And(n, two)
       q.Sub(t, one)
       u.Mul(u, u)
       v.Mul(v, v)
       t.And(n, one)
       if t.Cmp(one) == 0 {
           t.Sub(u, q)
           t.Mul(two, t)
           r.Mul(three, v)
           return u.Add(u, v), r.Sub(r, t)
       } else {
           t.Mul(three, u)
           r.Add(v, q)
           r.Mul(two, r)
           return r.Sub(r, t), u.Add(u, v)
       }
   }
   t, q, l := new(big.Int), new(big.Int), new(big.Int)
   u, v := inner(t.Rsh(n, 1))
   l.Mul(two, v)
   l.Sub(l, u) // Lucas function
   t.And(n, one)
   if t.Cmp(one) == 0 {
       q.And(n, two)
       q.Sub(q, one)
       t.Mul(v, l)
       return t.Add(t, q)
   }
   return u.Mul(u, l)

}

func main() {

   start := time.Now()
   n := new(big.Int)
   n.Lsh(one, 32)
   s := lucas(n).String()
   fmt.Printf("The digits of the 2^32nd Fibonacci number (%d) are:\n", len(s))
   fmt.Printf("  First twenty : %s\n", s[0:20])
   fmt.Printf("  Final twenty : %s\n", s[len(s)-20:])
   fmt.Printf("\nTook %s\n", time.Since(start))

}</lang>

Output:

The digits of the 2^32nd Fibonacci number (897595080) are:
  First twenty : 61999319689381859818
  Final twenty : 39623735538208076347

Took 5m2.722505927s

Julia

Because Julia uses the GMP library for its BigInt type, a BigInt cannot be larger than about 2^(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-nd 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]

binetfibonacci(n) = ((1+sqrt(big"5"))^n-(1-sqrt(big"5"))^n)/(sqrt(big"5")*big"2"^n) # not used here

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

function lucasfibonacci(n)

   function inner(n)
       if n == 0
           return big"0", big"1"
       end
       u, v = inner(n >> 1)
       q = (n & 2) - 1
       u *= u
       v *= v
       return isodd(n) ? BigInt(u + v), BigInt(3 * v - 2 * (u - q)) :
           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(64)), 45), rpad(firstlast(l2s(64)), 45))

</lang>

Output:

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