I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# 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:

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

Write a program using matrix exponentiation to generate Fibonacci(n) for n equal to: 10, 100, 1_000, 10_000, 100_000, 1_000_000 and 10_000_000.

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

Extra

Generate Fibonacci(216 ), Fibonacci(232) and Fibonacci(264) using the same method or another one.

## C#

### Matrix exponentiation

using System;using System.IO;using System.Numerics;using System.Threading;using System.Diagnostics;using System.Globalization; namespace Fibonacci {    class Program    {        private static readonly BigInteger[,] F = { { BigInteger.One, BigInteger.One }, { BigInteger.One, BigInteger.Zero } };        private static NumberFormatInfo nfi  = new NumberFormatInfo { NumberGroupSeparator = "_" };        private static BigInteger[,] Multiply(in BigInteger[,] A, in BigInteger[,] B)        {            if (A.GetLength(1) != B.GetLength(0))            {                throw new ArgumentException("Illegal matrix dimensions for multiplication.");            }            var C = new BigInteger[A.GetLength(0), B.GetLength(1)];            for (int i = 0; i < A.GetLength(0); ++i)            {                for (int j = 0; j < B.GetLength(1); ++j)                {                    for (int k = 0; k < A.GetLength(1); ++k)                    {                        C[i, j] +=  A[i, k] * B[k, j];                    }                }            }            return C;        }        private static BigInteger[,] Power(in BigInteger[,] A, ulong n)        {            if (A.GetLength(1) != A.GetLength(0))            {                throw new ArgumentException("Not a square matrix.");            }            var C = new BigInteger[A.GetLength(0), A.GetLength(1)];            for (int i = 0; i < A.GetLength(0); ++i)            {                C[i, i] = BigInteger.One;            }            if (0 == n) return C;            var S = new BigInteger[A.GetLength(0), A.GetLength(1)];            for (int i = 0; i < A.GetLength(0); ++i)            {                for (int j = 0; j < A.GetLength(1); ++j)                {                    S[i, j] = A[i, j];                }            }            while (0 < n)            {                if (1 == n % 2) C = Multiply(C, S);                S = Multiply(S,S);                n /= 2;            }            return C;        }        public static BigInteger Fib(in ulong n)        {            var C = Power(F, n);            return C[0, 1];        }        public static void Task(in ulong p)        {            var ans = Fib(p).ToString();            var sp = p.ToString("N0", nfi);            if (ans.Length <= 40)            {                Console.WriteLine("Fibonacci({0}) = {1}", sp, ans);            }            else            {                Console.WriteLine("Fibonacci({0}) = {1} ... {2}", sp, ans[0..19], ans[^20..]);            }        }        public static void Main()        {            Stopwatch stopWatch = new Stopwatch();            stopWatch.Start();            for (ulong p = 10; p <= 10_000_000; p *= 10) {                Task(p);            }            stopWatch.Stop();            TimeSpan ts = stopWatch.Elapsed;            string elapsedTime = String.Format("{0:00}:{1:00}:{2:00}.{3:00}",                ts.Hours, ts.Minutes, ts.Seconds,                ts.Milliseconds / 10);            Console.WriteLine("Took " + elapsedTime);        }    }}
Output:
Fibonacci(10) = 55
Fibonacci(100) = 354224848179261915075
Fibonacci(1_000) = 4346655768693745643 ... 76137795166849228875
Fibonacci(10_000) = 3364476487643178326 ... 66073310059947366875
Fibonacci(100_000) = 2597406934722172416 ... 49895374653428746875
Fibonacci(1_000_000) = 1953282128707757731 ... 68996526838242546875
Fibonacci(10_000_000) = 1129834378225399760 ... 86998673686380546875
Took 00:04:00.92


## Go

### Matrix exponentiation

This uses matrix exponentiation to calculate the (2^16)th and (2^32)nd Fibonacci numbers the last of 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.

package main import (    "fmt"    big "github.com/ncw/gmp"    "time") type vector = []*big.Inttype 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 commatize(n uint64) string {    s := fmt.Sprintf("%d", n)    le := len(s)    for i := le - 3; i >= 1; i -= 3 {        s = s[0:i] + "," + s[i:]    }    return s} func main() {    start := time.Now()    n := new(big.Int)    for i := uint64(10); i <= 1e7; i *= 10 {        n.SetUint64(i)        s := fibonacci(n).String()        fmt.Printf("The digits of the %sth Fibonacci number (%s) are:\n",            commatize(i), commatize(uint64(len(s))))        if len(s) > 20 {            fmt.Printf("  First 20 : %s\n", s[0:20])            if len(s) < 40 {                fmt.Printf("  Final %-2d : %s\n", len(s)-20, s[20:])            } else {                fmt.Printf("  Final 20 : %s\n", s[len(s)-20:])            }        } else {            fmt.Printf("  All %-2d   : %s\n", len(s), s)        }        fmt.Println()    }     sfxs := []string{"nd", "th"}    for i, e := range []uint{16, 32} {        n.Lsh(one, e)        s := fibonacci(n).String()        fmt.Printf("The digits of the 2^%d%s Fibonacci number (%s) are:\n", e, sfxs[i],            commatize(uint64(len(s))))        fmt.Printf("  First 20 : %s\n", s[0:20])        fmt.Printf("  Final 20 : %s\n", s[len(s)-20:])        fmt.Println()    }     fmt.Printf("Took %s\n\n", time.Since(start))}
Output:

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

The digits of the 10th Fibonacci number (2) are:
All 2    : 55

The digits of the 100th Fibonacci number (21) are:
First 20 : 35422484817926191507
Final 1  : 5

The digits of the 1,000th Fibonacci number (209) are:
First 20 : 43466557686937456435
Final 20 : 76137795166849228875

The digits of the 10,000th Fibonacci number (2,090) are:
First 20 : 33644764876431783266
Final 20 : 66073310059947366875

The digits of the 100,000th Fibonacci number (20,899) are:
First 20 : 25974069347221724166
Final 20 : 49895374653428746875

The digits of the 1,000,000th Fibonacci number (208,988) are:
First 20 : 19532821287077577316
Final 20 : 68996526838242546875

The digits of the 10,000,000th Fibonacci number (2,089,877) are:
First 20 : 11298343782253997603
Final 20 : 86998673686380546875

The digits of the 2^16nd Fibonacci number (13,696) are:
First 20 : 73199214460290552832
Final 20 : 97270190955307463227

The digits of the 2^32th Fibonacci number (897,595,080) are:
First 20 : 61999319689381859818
Final 20 : 39623735538208076347

Took 11m0.255916206s


### 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
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 new(big.Int), big.NewInt(1)        }        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 commatize(n uint64) string {    s := fmt.Sprintf("%d", n)    le := len(s)    for i := le - 3; i >= 1; i -= 3 {        s = s[0:i] + "," + s[i:]    }    return s} func main() {    start := time.Now()    n := new(big.Int)    for i := uint64(10); i <= 1e7; i *= 10 {        n.SetUint64(i)        s := lucas(n).String()        fmt.Printf("The digits of the %sth Fibonacci number (%s) are:\n",            commatize(i), commatize(uint64(len(s))))        if len(s) > 20 {            fmt.Printf("  First 20 : %s\n", s[0:20])            if len(s) < 40 {                fmt.Printf("  Final %-2d : %s\n", len(s)-20, s[20:])            } else {                fmt.Printf("  Final 20 : %s\n", s[len(s)-20:])            }        } else {            fmt.Printf("  All %-2d   : %s\n", len(s), s)        }        fmt.Println()    }     sfxs := []string{"nd", "th"}    for i, e := range []uint{16, 32} {        n.Lsh(one, e)        s := lucas(n).String()        fmt.Printf("The digits of the 2^%d%s Fibonacci number (%s) are:\n", e, sfxs[i],            commatize(uint64(len(s))))        fmt.Printf("  First 20 : %s\n", s[0:20])        fmt.Printf("  Final 20 : %s\n", s[len(s)-20:])        fmt.Println()    }     fmt.Printf("Took %s\n\n", time.Since(start))}
Output:
As first version except time is now 5m4.13427997s.


### 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
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 commatize(n uint64) string {    s := fmt.Sprintf("%d", n)    le := len(s)    for i := le - 3; i >= 1; i -= 3 {        s = s[0:i] + "," + s[i:]    }    return s} func main() {    start := time.Now()    n := new(big.Int)       for i := uint64(10); i <= 1e7; i *= 10 {        n.SetUint64(i)        nn := new(gmp.Int)           nn.SetUint64(i)        fmt.Printf("\nThe digits of the %sth Fibonacci number are:\n", commatize(i))        nd2, nd3 := nd, nd        // These need to be preset for i == 10 & i == 100        // as there is no way of deriving the total length of the string using this method.        if i == 10 {            nd2 = 2        } else if i == 100 {            nd3 = 1        }        s1 := firstFibDigits(n, nd2)               if len(s1) < 20 {            fmt.Printf("  All %-2d   : %s\n", len(s1), s1)        } else {            fmt.Printf("  First 20 : %s\n", s1)            s2 := lastFibDigits(nn, nd3)            if len(s2) < 20 {                fmt.Printf("  Final %-2d : %s\n", len(s2), s2)            } else {                fmt.Printf("  Final 20 : %s\n", s2)            }        }     }     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))}
Output:
The digits of the 10th Fibonacci number are:
All 2    : 55

The digits of the 100th Fibonacci number are:
First 20 : 35422484817926191507
Final 1  : 5

The digits of the 1,000th Fibonacci number are:
First 20 : 43466557686937456435
Final 20 : 76137795166849228875

The digits of the 10,000th Fibonacci number are:
First 20 : 33644764876431783266
Final 20 : 66073310059947366875

The digits of the 100,000th Fibonacci number are:
First 20 : 25974069347221724166
Final 20 : 49895374653428746875

The digits of the 1,000,000th Fibonacci number are:
First 20 : 19532821287077577316
Final 20 : 68996526838242546875

The digits of the 10,000,000th Fibonacci number are:
First 20 : 11298343782253997603
Final 20 : 86998673686380546875

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 6.391894ms


### Matrix exponentiation

import System.CPUTime (getCPUTime)import Data.List main = do    startTime <- getCPUTime    mapM_ (putStrLn.formatAns).take 7.iterate (*10) $10 mapM_ (putStrLn.seeFib) [16,32] finishTime <- getCPUTime putStrLn$ "Took " ++ (took startTime finishTime) took t = fromChrono.chrono t fromChrono :: (Integer,Integer,Integer) -> StringfromChrono (m,s,ms) = show m ++ "m" ++ show s ++ "." ++ show ms ++ "s" chrono :: Integer -> Integer -> (Integer,Integer,Integer)chrono start end = (m,s,ms)    where    tera = 1000000000000    fdt = fromIntegral (end - start) / tera    dt = floor fdt    (m,s) = quotRem dt 60     ms = round $fromIntegral (round (fdt - (fromIntegral dt))*1000) / 1000 bagOf :: Int -> [a] -> [[a]]bagOf _ [] = []bagOf n xs = let (us,vs) = splitAt n xs in us : bagOf n vs formatIntegral :: Show a => String -> a -> StringformatIntegral sep = reverse.intercalate sep.bagOf 3.reverse.show formatAns :: Integer -> StringformatAns p = start ++ go x where start = "Fibonacci("++ (formatIntegral "_" p) ++ ") = " x = fib p tenPow20 = 10^20 tenPow40 = tenPow20^2 go u | u <= tenPow20 = show u go u | u <= tenPow40 = let (us,vs) = splitAt 20$ show u in us ++ " ... " ++ vs    go u = (take 20 $show u) ++ " ... " ++ (show . rem u$ 10^20) seeFib :: Integer -> StringseeFib 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 -> Integerfib 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 1fib 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
Output:

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

Fibonacci(10) = 55
Fibonacci(100) = 35422484817926191507 ... 5
Fibonacci(1_000) = 43466557686937456435 ... 76137795166849228875
Fibonacci(10_000) = 33644764876431783266 ... 66073310059947366875
Fibonacci(100_000) = 25974069347221724166 ... 49895374653428746875
Fibonacci(1_000_000) = 19532821287077577316 ... 68996526838242546875
Fibonacci(10_000_000) = 11298343782253997603 ... 86998673686380546875
Fibonacci(2^16) = 73199214460290552832 ... 97270190955307463227
Fibonacci(2^32) = 61999319689381859818 ... 39623735538208076347
Took 5m20.1s


### Matrix exponentiation - printing alternative

import System.CPUTime (getCPUTime)import Data.List main = do    startTime <- getCPUTime    mapM_ (putStrLn.formatAns).take 7.iterate (*10) $10 mapM_ (putStrLn.seeFib) [16,32] finishTime <- getCPUTime putStrLn$ "Took " ++ (took startTime finishTime) took t = fromChrono.chrono t fromChrono :: (Integer,Integer,Integer) -> StringfromChrono (m,s,ms) = show m ++ "m" ++ show s ++ "." ++ show ms ++ "s" chrono :: Integer -> Integer -> (Integer,Integer,Integer)chrono start end = (m,s,ms)    where    tera = 1000000000000    fdt = fromIntegral (end - start) / tera    dt = floor fdt    (m,s) = quotRem dt 60     ms = round $fromIntegral (round (fdt - (fromIntegral dt))*1000) / 1000 bagOf :: Int -> [a] -> [[a]]bagOf _ [] = []bagOf n xs = let (us,vs) = splitAt n xs in us : bagOf n vs formatIntegral :: Show a => String -> a -> StringformatIntegral sep = reverse.intercalate sep.bagOf 3.reverse.show formatAns :: Integer -> StringformatAns p = start ++ (startEnd 20 (sizeFib p) num) where start = "Fibonacci("++ (formatIntegral "_" p) ++ ") = " num = fib p seeFib :: Integer -> StringseeFib n = start ++ (startEnd 20 (sizeFib p) num) where start = "Fibonacci(2^" ++ (show n) ++") = " p = 2^n num = fib p startEnd :: (Integral a, Show a) => a -> a -> a -> StringstartEnd ndigit len num | len <= ndigit = show numstartEnd ndigit len num | len <= 2*ndigit = let (us,vs) = genericSplitAt ndigit (show num) in us ++ " ... " ++ vsstartEnd ndigit len num = start ++ " ... " ++ end where end = show.rem num$ 10 ^ ndigit    start = show.quot num $10 ^ (len - ndigit + 1) phi :: Doublephi = (1 + sqrt 5)/2log10phi = logBase 10 phihalflog10five =(logBase 10 5)/2 sizeFib :: Integral a => a -> asizeFib p = ceiling$ (fromIntegral p)*log10phi - halflog10five fib :: Integer -> Integerfib 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 1fib 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
Output:

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

Fibonacci(10) = 55
Fibonacci(100) = 35422484817926191507 ... 5
Fibonacci(1_000) = 4346655768693745643 ... 76137795166849228875
Fibonacci(10_000) = 3364476487643178326 ... 66073310059947366875
Fibonacci(100_000) = 2597406934722172416 ... 49895374653428746875
Fibonacci(1_000_000) = 1953282128707757731 ... 68996526838242546875
Fibonacci(10_000_000) = 1129834378225399760 ... 86998673686380546875
Fibonacci(2^16) = 7319921446029055283 ... 97270190955307463227
Fibonacci(2^32) = 6199931968938185981 ... 39623735538208076347
Took 3m58.1s


### Matrix exponentiation for a symmetric matrix

We will use a property of symmetric matrices which commute.

If X,Y are two symmetric matrices of same size and if they commute then X*Y is a symmetric matrix.

It means to compute Z = X*Y, only terms on and below the diagonal need to be computed (above = below).

At each step of the exponentiation of a symmetric matric, we multiply 2 symmetric matrices which commute.

-- https://yutsumura.com/symmetric-matrices-and-the-product-of-two-matrices/import System.CPUTime (getCPUTime)import Data.List main = do    startTime <- getCPUTime    mapM_ (putStrLn.formatAns).take 7.iterate (*10) $10 mapM_ (putStrLn.seeFib) [16,32] finishTime <- getCPUTime putStrLn$ "Took " ++ (took startTime finishTime) took t = fromChrono.chrono t fromChrono :: (Integer,Integer,Integer) -> StringfromChrono (m,s,ms) = show m ++ "m" ++ show s ++ "." ++ show ms ++ "s" chrono :: Integer -> Integer -> (Integer,Integer,Integer)chrono start end = (m,s,ms)    where    tera = 1000000000000    fdt = fromIntegral (end - start) / tera    dt = floor fdt    (m,s) = quotRem dt 60     ms = round $fromIntegral (round (fdt - (fromIntegral dt))*1000) / 1000 bagOf :: Int -> [a] -> [[a]]bagOf _ [] = []bagOf n xs = let (us,vs) = splitAt n xs in us : bagOf n vs formatIntegral :: Show a => String -> a -> StringformatIntegral sep = reverse.intercalate sep.bagOf 3.reverse.show formatAns :: Integer -> StringformatAns p = start ++ (startEnd 20 (sizeFib p) num) where start = "Fibonacci("++ (formatIntegral "_" p) ++ ") = " num = fib p seeFib :: Integer -> StringseeFib n = start ++ (startEnd 20 (sizeFib p) num) where start = "Fibonacci(2^" ++ (show n) ++") = " p = 2^n num = fib p startEnd :: (Integral a, Show a) => a -> a -> a -> StringstartEnd ndigit len num | len <= ndigit = show numstartEnd ndigit len num | len <= 2*ndigit = let (us,vs) = genericSplitAt ndigit (show num) in us ++ " ... " ++ vsstartEnd ndigit len num = start ++ " ... " ++ end where end = show.rem num$ 10 ^ ndigit    start = show.quot num $10 ^ (len - ndigit + 1) phi :: Doublephi = (1 + sqrt 5)/2log10phi = logBase 10 phihalflog10five =(logBase 10 5)/2 sizeFib :: Integral a => a -> asizeFib p = ceiling$ (fromIntegral p)*log10phi - halflog10five fib :: Integer -> Integerfib 0 = 0fib n = zeroOne (power (Mat 1 1 1 0)  n) data Mat a = Mat {zeroZero :: a, zeroOne :: a, oneZero :: a, oneOne :: a} deriving (Eq,Show) -- for a symmetric matrixsquare :: Num a => (Mat a) -> (Mat a)square (Mat x00 x01 x10 x11) = Mat y00 y10 y10 y11    where    y00 = y10 + y11 -- F_{n+1} = F_{n} + F_{n-1}    y10 = x10*(x00+x11)    y11 = x11*x11+x10*x10 -- for 2 symmetric matrices which commutemult :: Num a => (Mat a) -> (Mat a) -> (Mat a)mult (Mat x00 x01 x10 x11) (Mat y00 y01 y10 y11) = Mat xy00 xy01 xy01 xy11    where    xy00 =  xy01 + xy11 -- F_{n+1} = F_{n} + F_{n-1}    xy01 = x10*y00 + x11*y10    xy11 = x10*y01 + x11*y11 power :: Num a => (Mat a) -> Integer -> (Mat a)power _ n | n < 0 = error "Exception: Negative exponent"power _ 0 = Mat 1 0 0 1power m 1 = mpower m n = if even n then w else mult w m   where w = square.power m.quot n $2 Output: Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5: Fibonacci(10) = 55 Fibonacci(100) = 35422484817926191507 ... 5 Fibonacci(1_000) = 4346655768693745643 ... 76137795166849228875 Fibonacci(10_000) = 3364476487643178326 ... 66073310059947366875 Fibonacci(100_000) = 2597406934722172416 ... 49895374653428746875 Fibonacci(1_000_000) = 1953282128707757731 ... 68996526838242546875 Fibonacci(10_000_000) = 1129834378225399760 ... 86998673686380546875 Fibonacci(2^16) = 7319921446029055283 ... 97270190955307463227 Fibonacci(2^32) = 6199931968938185981 ... 39623735538208076347 Took 2m6.1s ## Java Performed the task to use Matrix multiplication to compute Fibonacci numbers. Implemented fib and fibMod.  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); } }  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. # 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)endlastfibmod(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 * lend 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))  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 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;}
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: Phix/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.)

include mpfr.empfr_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):=0setd({1,-1},mpz_init(1),cache) --     fib(1):=1setd({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 ifend 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
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.

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 resend 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 cend 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 resend function function fibonacci(integer n)    sequence m = {{ONE, ONE},                  {ONE, ZERO}}    m = matrix_exponent(m,n-1)    mpz res = m[1][1]    return resend 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
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.

## Raku

(formerly 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.

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 exponentiationsub 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 FatRatmulti sqrt(FatRat$r --> FatRat) {    FatRat.new: sqrt($r.nude[0] × 10**(D×2) div$r.nude[1]), 10**D} # square-root: accepts/return Intmulti 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 FatRatsub 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' versionsconstant 𝑒 = &𝑒();constant π = &π(); # approximation of natural log, accepts any numeric, returns FatRatsub 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 FatRatmulti 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 feasiblefor (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 casefor 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}
Output:
F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347
F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347

## Sidef

Computing the n-th Fibonacci number, using matrix-exponentiation (this function is also built-in):

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]

First and last 20 digits of the n-th Fibonacci number:

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}"    }}
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:

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}"}
Output:
F(2^16) = 73199214460290552832 ... 97270190955307463227
F(2^32) = 61999319689381859818 ... 39623735538208076347
F(2^64) = 11175807536929528424 ... 17529800348089840187