Fibonacci matrix-exponentiation

The Fibonacci sequence defined with matrix-exponentiation:

      F₀ = 0 
      F₁ = 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

Go

Matrix exponentiation

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)th Fibonacci number which appears to be well out of reach using this approach.

Library: GMP(Go wrapper)

<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

Lucas method

Although Go supports big.Float, the precision needed to calculate the (2^32)nd 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

Fibmod method

This uses the Sidef entry's 'Fibmod' approach to enable the (2^64)th Fibonacci number to be processed. As Go lacks such a function, I have translated the Julia version. I have also had to pull in a third party library to provide functions (such as Log) which Go's big.Float implementation lacks.

The speed-up compared to the other approaches is astonishing!

Library: bigfloat

Translation of: Sidef

Translation of: Julia

<lang go>package main

import (

   "fmt"
   "github.com/ALTree/bigfloat"
   "github.com/ncw/gmp"
   "math/big"
   "time"

)

const (

   nd = 20  // number of digits to be displayed at each end
   pr = 128 // precision to be used

)

var (

   one  = gmp.NewInt(1)
   two  = gmp.NewInt(2)
   ten  = gmp.NewInt(10)
   onef = big.NewFloat(1).SetPrec(pr)
   tenf = big.NewFloat(10).SetPrec(pr)
   ln10 = bigfloat.Log(tenf)

)

func fibmod(n, nmod *gmp.Int) *gmp.Int {

   if n.Cmp(two) < 0 {
       return n
   }
   fibmods := make(map[string]*gmp.Int)
   var f func(n *gmp.Int) *gmp.Int
   f = func(n *gmp.Int) *gmp.Int {
       if n.Cmp(two) < 0 {
           return one
       }
       ns := n.String()
       if v, ok := fibmods[ns]; ok {
           return v
       }
       k, t, u, v := new(gmp.Int), new(gmp.Int), new(gmp.Int), new(gmp.Int)
       k.Quo(n, two)
       t.And(n, one)
       if t.Cmp(one) != 0 {
           t.Set(f(k))
           t.Mul(t, t)
           v.Sub(k, one)
           u.Set(f(v))
           u.Mul(u, u)
       } else {
           t.Set(f(k))
           v.Add(k, one)
           v.Set(f(v))
           u.Sub(k, one)
           u.Set(f(u))
           u.Mul(u, t)
           t.Mul(t, v)
       }
       t.Add(t, u)
       fibmods[ns] = t.Rem(t, nmod)
       return fibmods[ns]
   }
   w := new(gmp.Int)
   w.Sub(n, one)
   return f(w)

}

func binetApprox(n *big.Int) *big.Float {

   phi, ihp := big.NewFloat(0.5).SetPrec(pr), big.NewFloat(0.5).SetPrec(pr)
   root := big.NewFloat(1.25).SetPrec(pr)
   root.Sqrt(root)
   phi.Add(root, phi)
   ihp.Sub(root, ihp)
   ihp.Neg(ihp)
   ihp.Sub(phi, ihp)
   ihp = bigfloat.Log(ihp)
   phi = bigfloat.Log(phi)
   nn := new(big.Float).SetPrec(pr).SetInt(n)
   phi.Mul(phi, nn)
   return phi.Sub(phi, ihp)

}

func firstFibDigits(n *big.Int, k int) string {

   f := binetApprox(n)
   g := new(big.Float).SetPrec(pr)
   g.Quo(f, ln10)
   g.Add(g, onef)
   i, _ := g.Int(nil)
   g.SetInt(i)
   g.Mul(ln10, g)
   f.Sub(f, g)
   f = bigfloat.Exp(f)
   p := big.NewInt(int64(k))
   p.Exp(big.NewInt(10), p, nil)
   g.SetInt(p)
   f.Mul(f, g)
   i, _ = f.Int(nil)
   return i.String()[0:k]

}

func lastFibDigits(n *gmp.Int, k int) string {

   p := gmp.NewInt(int64(k))
   p.Exp(ten, p, nil)
   return fibmod(n, p).String()[0:k]

}

func main() {

   start := time.Now()
   n := new(big.Int)
   o := big.NewInt(1)
   ord := []string{"th", "nd", "th"}
   for i, p := range []uint{16, 32, 64} {
       n.Lsh(o, p)
       nn := new(gmp.Int)
       nn.Lsh(one, p)
       fmt.Printf("\nThe digits of the 2^%d%s Fibonacci number are:\n", p, ord[i])
       fmt.Printf("  First %d : %s\n", nd, firstFibDigits(n, nd))
       fmt.Printf("  Final %d : %s\n", nd, lastFibDigits(nn, nd))
   }
   fmt.Printf("\nTook %s\n", time.Since(start))

}</lang>

Output:

The digits of the 2^16th Fibonacci number are:
  First 20 : 73199214460290552832
  Final 20 : 97270190955307463227

The digits of the 2^32nd Fibonacci number are:
  First 20 : 61999319689381859818
  Final 20 : 39623735538208076347

The digits of the 2^64th Fibonacci number are:
  First 20 : 11175807536929528424
  Final 20 : 17529800348089840187

Took 3.244828ms

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># Here is the matrix Fibonacci formula as specified to be used for the solution. const b = [big"1" 1; 1 0] matrixfibonacci(n) = n == 0 ? 0 : n == 1 ? 1 : (b^(n+1))[2,2]

This exact Binet Fibonacci formula is not used due to BigFloat exponent size limitations.

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

Use the exponent size limiting variant of the Binet formula seen in the Sidef example.

function firstbinet(bits, ndig=20)

   logφ =  big"2"^bits * log(10, (1 + sqrt(BigFloat(5.0))) / 2)
   mantissa = logφ - trunc(logφ) + ndig + 1
   return string(BigInt(round((10^mantissa - 10^(-mantissa)) / sqrt(BigFloat(5.0)))))[1:ndig]

end

The fibmod function has no builtin in Julia, so here is one.

function fibmod(n::BigInt, nmod::BigInt)

   n < 2 && return n
   fibmods = Dict{BigInt, BigInt}()
   function f(n::BigInt)
       n < 2 && return 1
       haskey(fibmods, n) && return fibmods[n]
       k = div(n, 2)
       fibmods[n] = iseven(n) ?
           (f(k) * f(k) + f(k - 1) * f(k - 1)) % nmod :
           (f(k) * f(k + 1) + f(k - 1) * f(k)) % nmod
   end
   f(n - 1)

end lastfibmod(bits, ndig=21) = string(fibmod(big"2"^bits, big"10"^(ndig+1)))

See Wikipedia on Lucas function for the algorithm below.
inner -> F(n/2), F(n/2 - 1), L(n) = F(n) + 2F(n-1), and L(n/2) * F(n/2) = F(n)

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) < 40 ? s : s[1:20] * "..." * s[end-20+1:end])

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

   rpad(firstlast(firstbinet(32) * lastfibmod(32)), 45))

println("2^64 ", " "^90, rpad(firstlast(firstbinet(64) * lastfibmod(64)), 45))

</lang>

Output:

N                       Matrix                                        Lucas                                        Mod
-------------------------------------------------------------------------------------------------------------------------------------------------
2^32  61999319689381859818...39623735538208076347  61999319689381859818...39623735538208076347  61999319689381859818...39623735538208076347
2^64                                                                                            11175807536929528424...17529800348089840187

Perl

Translation of: Sidef

<lang perl>use strict; use warnings;

use Math::AnyNum qw(:overload fibmod floor); use Math::MatrixLUP;

sub fibonacci {

   my $M = Math::MatrixLUP->new([ [1, 1], [1,0] ]);
   (@{$M->pow(shift)})[0][1]

}

for my $n (16, 32) {

   my $f = fibonacci(2**$n);
   printf "F(2^$n) = %s ... %s\n",  substr($f,0,20), $f % 10**20;

}

sub binet_approx {

   my($n) = @_;
   use constant PHI =>   sqrt(1.25) + 0.5 ;
   use constant IHP => -(sqrt(1.25) - 0.5);
   (log(PHI)*$n - log(PHI-IHP))

}

sub nth_fib_first_k_digits {

   my($n,$k) = @_;
   my $f = binet_approx($n);
   floor(exp($f - log(10)*(floor($f / log(10) + 1))) * 10**$k)

}

sub nth_fib_last_k_digits {

   my($n,$k) = @_;
   fibmod($n, 10**$k);

}

print "\n"; for my $n (16, 32, 64) {

   my $first20 = nth_fib_first_k_digits(2**$n, 20);
   my $last20  = nth_fib_last_k_digits(2**$n, 20);
   printf "F(2^$n) = %s ... %s\n", $first20, $last20;

}</lang>

Output:

F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347

F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347
F(2^64) = 11175807536929528424 ... 17529800348089840187

Phix

Library: mpfr

Translation of: Sidef

Since I don't have a builtin fibmod, I had to roll my own, and used {n,m} instead of(/to mean) 2^n+m, thus avoiding some 2^53 native atom limits on 32-bit.

(mpz and mpfr variables are effectively pointers, and therefore simply won't work as expected/needed should you try and use them as keys to a cache.) <lang Phix>include mpfr.e mpfr_set_default_prec(-40)

constant PHI = mpfr_init(1.25),

        IHP = mpfr_init(),
        HLF = mpfr_init(0.5),
        L10 = mpfr_init(10)

mpfr_sqrt(PHI,PHI) mpfr_sub(IHP,PHI,HLF) mpfr_add(PHI,PHI,HLF) mpfr_neg(IHP,IHP) mpfr_log(L10,L10)

procedure binet_approx(mpfr r, integer n)

   mpfr m = mpfr_init()
   mpfr_ui_pow_ui(m,2,n) -- (n as in 2^n here)
   mpfr_log(r,PHI)
   mpfr_mul(r,r,m)
   mpfr_sub(m,PHI,IHP)
   mpfr_log(m,m)
   mpfr_sub(r,r,m)
   m = mpfr_free(m)

end procedure

function nth_fib_first_k_digits(integer n, k)

   mpfr {f,g,h} = mpfr_inits(3)
   binet_approx(f,n)
   mpfr_div(g,f,L10)
   mpfr_add_si(g,g,1)
   mpfr_floor(g,g)
   mpfr_mul(g,L10,g)
   mpfr_sub(f,f,g)
   mpfr_exp(f,f)
   mpfr_ui_pow_ui(h,10,k)
   mpfr_mul(f,f,h)
   mpfr_floor(f,f)
   string fmt = sprintf("%%%d.0Rf",k)
   return mpfr_sprintf(fmt,f)

end function

integer cache = new_dict() -- key of {n,m} means 2^n+m (where n and m are integers, n>=0, m always 0 or -1) setd({0,0},mpz_init(0),cache) -- aka fib(0):=0 setd({1,-1},mpz_init(1),cache) -- fib(1):=1 setd({1,0},mpz_init(1),cache) -- fib(2):=1

procedure mpz_fibmod(mpz f, p, integer n, m=0)

   sequence key = {n,m}
   if getd_index(key,cache)!=NULL then
       mpz_set(f,getd(key,cache))
   else
       mpz {f1,f2} = mpz_inits(2)
       n -= 1
       mpz_fibmod(f1,p,n)
       mpz_fibmod(f2,p,n,-1)
       if m=-1 then
           -- fib(2n-1) = fib(n)^2 + fib(n-1)^2
           mpz_mul(f1,f1,f1)
           mpz_mul(f2,f2,f2)
           mpz_add(f1,f1,f2)
       else
           -- fib(2n) = (2*fib(n-1)+fib(n))*fib(n)
           mpz_mul_si(f2,f2,2)
           mpz_add(f2,f2,f1)
           mpz_mul(f1,f2,f1)
       end if
       mpz_mod(f1,f1,p) 
       setd(key,f1,cache)
       mpz_set(f,f1)   
   end if

end procedure

function nth_fib_last_k_digits(integer n, k)

   mpz {f,p} = mpz_inits(2)
   mpz_ui_pow_ui(p,10,k)
   mpz_fibmod(f,p,n)
   return mpz_get_str(f)

end function

constant tests = {16,32,64} for i=1 to length(tests) do

   integer n = tests[i]
   string first20 = nth_fib_first_k_digits(n, 20),
          last20 = nth_fib_last_k_digits(n, 20)
   printf(1,"F(2^%d) = %s ... %s\n",{n,first20,last20})

end for</lang>

Output:

F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347
F(2^64) = 11175807536929528424 ... 17529800348089840187

Sidef

Computing the n-th Fibonacci number, using matrix-exponentiation (this function is also built-in): <lang ruby>func fibonacci(n) {

   ([[1,1],[1,0]]**n)[0][1]

}

say 15.of(fibonacci) #=> [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377]</lang>

First and last 20 digits of the n-th Fibonacci number: <lang ruby>for n in (16, 32) {

   var f = fibonacci(2**n)

   with (20) {|k|
       var a = (f // 10**(f.ilog10 - k + 1))
       var b = (f % 10**k)
       say "F(2^#{n}) = #{a} ... #{b}"
   }

}</lang>

Output:

F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347

More efficient approach, using Binet's formula for computing the first k digits, combined with the built-in method fibmod(n,m) for computing the last k digits: <lang ruby>func binet_approx(n) {

   const PHI =  (1.25.sqrt + 0.5)
   const IHP = -(1.25.sqrt - 0.5)
   (log(PHI)*n - log(PHI-IHP))

}

func nth_fib_first_k_digits(n,k) {

   var f = binet_approx(n)
   floor(exp(f - log(10)*(floor(f / log(10) + 1))) * 10**k)

}

func nth_fib_last_k_digits(n,k) {

   fibmod(n, 10**k)

}

for n in (16, 32, 64) {

   var first20 = nth_fib_first_k_digits(2**n, 20)
   var last20  = nth_fib_last_k_digits(2**n, 20)
   say "F(2^#{n}) = #{first20} ... #{last20}"

}</lang>

Output:

F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347
F(2^64) = 11175807536929528424 ... 17529800348089840187