Fibonacci matrix-exponentiation

The Fibonacci sequence defined with matrix-exponentiation:

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

Task

Write a program using matrix exponentiation to generate Fibonacci numbers: Fibonacci(2¹⁶), Fibonacci(2³²) and Fibonacci(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 or another method.

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

Haskell

Matrix exponentiation

<lang Haskell>import System.CPUTime (getCPUTime)

seeFib :: Integer -> String seeFib n = start ++ xs ++ " ... " ++ (show . rem x $ (10^20))

   where
   start = "Fibonacci(2^" ++(show n) ++") = "
   x = fib (2^n)
   xs = take 20 $ show x

fib :: Integer -> Integer fib 0 = 0 -- this line is necessary because "something ^ 0" returns "fromInteger 1", which unfortunately -- in our case is not our multiplicative identity (the identity matrix) but just a 1x1 matrix of 1 fib n = (last . head . unMat) (Mat [[1, 1], [1, 0]] ^ n)

mult :: Num a => a -> a -> a mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs) . zipWith (flip (map . (*))) vss) uss

newtype Mat a = Mat

 { unMat :: a
 } deriving (Eq,Show)

instance Num a => Num (Mat a) where

 negate xm = Mat $ map (map negate) $ unMat xm
 xm + ym = Mat $ zipWith (zipWith (+)) (unMat xm) (unMat ym)
 xm * ym =  Mat $ mult (unMat xm) (unMat ym)
 fromInteger n = Mat fromInteger n
 abs = undefined
 signum = undefined

main :: IO () main =

   do
   startTime <- getCPUTime
   mapM_ (putStrLn.seeFib) [16,32]
   finishTime <- getCPUTime
   let diff = ceiling $ fromIntegral (finishTime - startTime) / 1000000000000
   let (minute,second) = quotRem diff 60
   putStrLn $ "Took: " ++ (show minute) ++ "min" ++ (show second) ++ "s"

</lang>

Output:

Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5:

Fibonacci(2^16) = 73199214460290552832 ... 97270190955307463227
Fibonacci(2^32) = 61999319689381859818 ... 39623735538208076347
Took: 5min31s

Matrix exponentiation - printing alternative

<lang Haskell>import System.CPUTime (getCPUTime)

phi :: Double phi = (1 + sqrt 5)/2 log10phi = logBase 10 phi halflog10five = (logBase 10 5)/2

seeFib :: Integer -> String seeFib n = start ++ xs

   where
   start = "Fibonacci(2^" ++(show n) ++") = "
   p = 2^n
   x = fib p
   fp = fromIntegral p
   len = floor $ fp*log10phi - halflog10five
   xs = if len <= 40 then (show x) else (show . quot x $ 10^(len-19)) ++ " ... " ++ (show . rem x $ 10^20)

fib :: Integer -> Integer fib 0 = 0 -- this line is necessary because "something ^ 0" returns "fromInteger 1", which unfortunately -- in our case is not our multiplicative identity (the identity matrix) but just a 1x1 matrix of 1 fib n = (last . head . unMat) (Mat [[1, 1], [1, 0]] ^ n)

mult :: Num a => a -> a -> a mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs) . zipWith (flip (map . (*))) vss) uss

newtype Mat a = Mat

 { unMat :: a
 } deriving (Eq,Show)

instance Num a => Num (Mat a) where

 negate xm = Mat $ map (map negate) $ unMat xm
 xm + ym = Mat $ zipWith (zipWith (+)) (unMat xm) (unMat ym)
 xm * ym =  Mat $ mult (unMat xm) (unMat ym)
 fromInteger n = Mat fromInteger n
 abs = undefined
 signum = undefined

main :: IO () main =

   do
   startTime <- getCPUTime
   mapM_ (putStrLn.seeFib) [16,32]
   finishTime <- getCPUTime
   let diff = ceiling $ fromIntegral (finishTime - startTime) / 1000000000000
   let (minute,second) = quotRem diff 60
   putStrLn $ "Took: " ++ (show minute) ++ "min" ++ (show second) ++ "s"

</lang>

Output:

Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5:

Fibonacci(2^16) = 73199214460290552832 ... 97270190955307463227
Fibonacci(2^32) = 61999319689381859818 ... 39623735538208076347
Took: 4min19s

Java

Performed the task to use Matrix multiplication to compute Fibonacci numbers.

Implemented fib and fibMod.

<lang java> import java.math.BigInteger; import java.util.Arrays;

public class FibonacciMatrixExponentiation {

   public static void main(String[] args) {
       BigInteger mod = BigInteger.TEN.pow(20);
       for ( int exp : Arrays.asList(32, 64) ) {
           System.out.printf("Last 20 digits of fib(2^%d) = %s%n", exp, fibMod(BigInteger.valueOf(2).pow(exp), mod));
       }
       
       for ( int i = 1 ; i <= 7 ; i++ ) {
           BigInteger n = BigInteger.TEN.pow(i);
           System.out.printf("fib(%,d) = %s%n", n, displayFib(fib(n)));
       }
   }
   
   private static String displayFib(BigInteger fib) {
       String s = fib.toString();
       if ( s.length() <= 40 ) {
           return s;
       }
       return s.substring(0, 20) + " ... " + s.subSequence(s.length()-20, s.length());
   }

   //  Use Matrix multiplication to compute Fibonacci numbers.
   private static BigInteger fib(BigInteger k) {
       BigInteger aRes = BigInteger.ZERO;
       BigInteger bRes = BigInteger.ONE;
       BigInteger cRes = BigInteger.ONE;
       BigInteger aBase = BigInteger.ZERO;
       BigInteger bBase = BigInteger.ONE;
       BigInteger cBase = BigInteger.ONE;
       while ( k.compareTo(BigInteger.ZERO) > 0 ) {
           if ( k.mod(BigInteger.valueOf(2)).compareTo(BigInteger.ONE) == 0 ) {
               BigInteger temp1 = aRes.multiply(aBase).add(bRes.multiply(bBase));
               BigInteger temp2 = aBase.multiply(bRes).add(bBase.multiply(cRes));
               BigInteger temp3 = bBase.multiply(bRes).add(cBase.multiply(cRes));
               aRes = temp1;
               bRes = temp2;
               cRes = temp3;
           }
           k = k.shiftRight(1);
           BigInteger temp1 = aBase.multiply(aBase).add(bBase.multiply(bBase));
           BigInteger temp2 = aBase.multiply(bBase).add(bBase.multiply(cBase));
           BigInteger temp3 = bBase.multiply(bBase).add(cBase.multiply(cBase));
           aBase = temp1;
           bBase = temp2;
           cBase = temp3;
       }
       return aRes;
   }

   //  Use Matrix multiplication to compute Fibonacci numbers.
   private static BigInteger fibMod(BigInteger k, BigInteger mod) {
       BigInteger aRes = BigInteger.ZERO;
       BigInteger bRes = BigInteger.ONE;
       BigInteger cRes = BigInteger.ONE;
       BigInteger aBase = BigInteger.ZERO;
       BigInteger bBase = BigInteger.ONE;
       BigInteger cBase = BigInteger.ONE;
       while ( k.compareTo(BigInteger.ZERO) > 0 ) {
           if ( k.mod(BigInteger.valueOf(2)).compareTo(BigInteger.ONE) == 0 ) {
               BigInteger temp1 = aRes.multiply(aBase).add(bRes.multiply(bBase)).mod(mod);
               BigInteger temp2 = aBase.multiply(bRes).add(bBase.multiply(cRes)).mod(mod);
               BigInteger temp3 = bBase.multiply(bRes).add(cBase.multiply(cRes)).mod(mod);
               aRes = temp1;
               bRes = temp2;
               cRes = temp3;
           }
           k = k.shiftRight(1);
           BigInteger temp1 = aBase.multiply(aBase).add(bBase.multiply(bBase)).mod(mod);
           BigInteger temp2 = aBase.multiply(bBase).add(bBase.multiply(cBase)).mod(mod);
           BigInteger temp3 = bBase.multiply(bBase).add(cBase.multiply(cBase)).mod(mod);
           aBase = temp1;
           bBase = temp2;
           cBase = temp3;
       }
       return aRes.mod(mod);
   }

} </lang>

Output:

Last 20 digits of fib(2^32) = 39623735538208076347
Last 20 digits of fib(2^64) = 17529800348089840187

fib(10) = 55
fib(100) = 354224848179261915075
fib(1,000) = 43466557686937456435 ... 76137795166849228875
fib(10,000) = 33644764876431783266 ... 66073310059947366875
fib(100,000) = 25974069347221724166 ... 49895374653428746875
fib(1,000,000) = 19532821287077577316 ... 68996526838242546875
fib(10,000,000) = 11298343782253997603 ... 86998673686380546875

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^16 ", rpad(firstlast(m2s(16)), 45), rpad(firstlast(l2s(16)), 45),

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

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^16  73199214460290552832...97270190955307463227  73199214460290552832...97270190955307463227  73199214460290552832...97270190955307463227
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

Perl 6

Following the general approach of Sidef, and relying on Perl for fibmod function. Borrowed/adapted routines from Ramanujan's_constant task to allow FatRat calculations throughout. Does not quite meet task spec, as n^64 results not working yet. <lang perl6>use Math::Matrix; use Inline::Perl5; my $p5 = Inline::Perl5.new(); $p5.use( 'Math::AnyNum' );

constant D = 53; # set the size of FatRat calcluations

matrix exponentiation

sub fibonacci ($n) {

   my $M = Math::Matrix.new( [[1,1],[1,0]] );
   ($M ** $n)[0][1]

}

calculation of 𝑒

sub postfix:<!> (Int $n) { (constant f = 1, |[\×] 1..*)[$n] } sub 𝑒 (--> FatRat) { sum map { FatRat.new(1,.!) }, ^D }

calculation of π

sub π (--> FatRat) {

   my ($a, $n, $g, $z, $pi) = 1, 1, sqrt(1/2.FatRat), 0.25;

   for ^5 {
       given [ ($a + $g)/2, sqrt $a × $g ] {
           $z -= (.[0] - $a)**2 × $n;
           $n += $n;
           ($a, $g) = @$_;
           $pi = ($a ** 2 / $z).substr: 0, 2 + D
       }
   }
   $pi.FatRat

}

square-root: accepts/return FatRat

multi sqrt(FatRat $r --> FatRat) {

   FatRat.new: sqrt($r.nude[0] × 10**(D×2) div $r.nude[1]), 10**D

}

square-root: accepts/return Int

multi sqrt(Int $n --> Int) {

   my $guess = 10**($n.chars div 2);
   my $iterator = { ( $^x   +   $n div ($^x) ) div 2 };
   my $endpoint = { $^x == $^y|$^z };
   min ($guess, $iterator … $endpoint)[*-1, *-2]

}

arithmetic-geometric mean: accepts/returns FatRat

sub AG-mean(FatRat $a is copy, FatRat $g is copy --> FatRat) {

   ($a, $g) = ($a + $g)/2, sqrt $a × $g until $a - $g < 10**-D;
   $a

}

override built-in definitions with 'FatRat' versions

constant 𝑒 = &𝑒(); constant π = &π();

approximation of natural log, accepts any numeric, returns FatRat

sub log-approx ($x --> FatRat) {

   constant ln2 = 2 * [+] (((1/3).FatRat**(2*$_+1))/(2*$_+1) for 0..D); # 1/3 = (2-1)/(2+1)
   π / (2 × AG-mean(1.FatRat, 2.FatRat**(2-D)/$x)) - D × ln2

}

power function, with exponent less than zero: accepts/returns FatRat

multi infix:<**> (FatRat $base, FatRat $exp is copy where * < 0 --> FatRat) {

   constant ε = 10**-D;
   my ($low, $high) = 0.FatRat, 1.FatRat;
   my $mid  = $high / 2;
   my $acc  = my $sqr = sqrt($base);
   $exp = -$exp;

   while (abs($mid - $exp) > ε) {
       $sqr = sqrt($sqr);
       if ($mid <= $exp) { $low  = $mid; $acc ×=   $sqr }
       else              { $high = $mid; $acc ×= 1/$sqr }
       $mid = ($low + $high) / 2
   }

   (1/$acc).substr(0, D).FatRat

}

sub binet_approx (Int $n --> FatRat) {

   constant PHI =   sqrt(1.25.FatRat) + 0.5 ;
   constant IHP = -(sqrt(1.25.FatRat) - 0.5);
   $n × log-approx(PHI) - log-approx(PHI - IHP)

}

sub nth_fib_first_k_digits ($n,$k) {

   my $f     = binet_approx($n);
   my $log10 = log-approx(10);
   floor( 𝑒**($f - $log10×(floor($f/$log10 + 1))) × 10**$k)

}

my &nth_fib_last_k_digits =

   $p5.run('sub { my($n,$k) = @_; Math::AnyNum::fibmod($n, 10**$k) }');

matrix exponentiation is very inefficient, n^64 not feasible

for (16, 32) -> $n {

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

}

this way is much faster, but not yet able to handle 2^64 case

for 16, 32 -> $n {

   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

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

matrix exponentiation (2^16)

Somewhat closer to the original specification, but certainly not recommended for 2^32, let alone 2^64...
Contains copies of routines from Matrix-exponentiation_operator#Phix, but modified to use gmp. <lang Phix>include mpfr.e

constant ZERO = mpz_init(0),

        ONE = mpz_init(1)

function identity(integer n)

   sequence res = repeat(repeat(ZERO,n),n)
   for i=1 to n do
       res[i][i] = ONE
   end for
   return res

end function

function matrix_mul(sequence a, sequence b)

   if length(a[1])!=length(b) then
       crash("matrices cannot be multiplied")
   end if
   sequence c = repeat(repeat(ZERO,length(b[1])),length(a))
   for i=1 to length(a) do
       for j=1 to length(b[1]) do
           for k=1 to length(a[1]) do
               mpz cij = mpz_init()
               mpz_mul(cij,a[i][k],b[k][j])
               mpz_add(cij,cij,c[i][j])
               c[i][j] = cij
           end for
       end for
   end for
   return c

end function

function matrix_exponent(sequence m, integer n)

   integer l = length(m)
   if l!=length(m[1]) then crash("not a square matrix") end if
   if n<0 then crash("negative exponents not supported") end if
   sequence res = identity(l)
   while n do
       if and_bits(n,1) then
           res = matrix_mul(res,m)
       end if
       n = floor(n/2)
       if n=0 then exit end if -- (avoid unnecessary matrix_mul)
       m = matrix_mul(m,m)
   end while
   return res

end function

function fibonacci(integer n)

   sequence m = {{ONE, ONE},
                 {ONE, ZERO}}
   m = matrix_exponent(m,n-1)
   mpz res = m[1][1]
   return res

end function

atom t0 = time() mpz fn = fibonacci(power(2,16)) string s = mpz_get_str(fn) string e = elapsed(time()-t0) printf(1,"fibonnaci(2^16) = %s [%s]\n",{shorten(s),e})

for i=1 to 7 do

   t0 = time()
   integer n = power(10,i)
   fn = fibonacci(n)
   s = mpz_get_str(fn)
   t0 = time()-t0
   e = iff(t0>0.1?" ["&elapsed(t0)&"]":"")
   printf(1,"fibonnaci(%,d) = %s%s\n",{n,shorten(s),e})

end for</lang>

Output:

fibonnaci(2^16) = 73199214460290552832...97270190955307463227 (13,696 digits) [0s]
fibonnaci(10) = 55
fibonnaci(100) = 354224848179261915075
fibonnaci(1,000) = 43466557686937456435...76137795166849228875 (209 digits)
fibonnaci(10,000) = 33644764876431783266...66073310059947366875 (2,090 digits)
fibonnaci(100,000) = 25974069347221724166...49895374653428746875 (20,899 digits)
fibonnaci(1,000,000) = 19532821287077577316...68996526838242546875 (208,988 digits) [0.2s]
fibonnaci(10,000,000) = 11298343782253997603...86998673686380546875 (2,089,877 digits) [9.0s]

Change that loop to 8 and a 9 year old 3.3GHz i3 also eventually gets:

fibonnaci(100,000,000) = 47371034734563369625...06082642167760546875 (20,898,764 digits) [14 minutes and 24s]

Clearly 2^32 (897 million digits, apparently) is a tad out of bounds, let alone 2^64.

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