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

From Rosetta Code
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:

Task

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.

Related tasks


C#[edit]

Matrix exponentiation[edit]

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(1), B.GetLength(0)];
for (int i = 0; i < A.GetLength(0); ++i)
{
for (int j = 0; j < B.GetLength(1); ++j)
{
for (int k = 0; k < B.GetLength(0); ++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[edit]

Matrix exponentiation[edit]

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.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 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[edit]

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[edit]

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

Haskell[edit]

Matrix exponentiation[edit]

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) -> String
fromChrono (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 -> String
formatIntegral sep = reverse.intercalate sep.bagOf 3.reverse.show
 
formatAns :: Integer -> String
formatAns 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 -> 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
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[edit]

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) -> String
fromChrono (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 -> String
formatIntegral sep = reverse.intercalate sep.bagOf 3.reverse.show
 
formatAns :: Integer -> String
formatAns p = start ++ (startEnd 20 (sizeFib p) num)
where
start = "Fibonacci("++ (formatIntegral "_" p) ++ ") = "
num = fib p
 
seeFib :: Integer -> String
seeFib 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 -> String
startEnd ndigit len num | len <= ndigit = show num
startEnd ndigit len num | len <= 2*ndigit = let (us,vs) = genericSplitAt ndigit (show num) in us ++ " ... " ++ vs
startEnd ndigit len num = start ++ " ... " ++ end
where
end = show.rem num $ 10 ^ ndigit
start = show.quot num $ 10 ^ (len - ndigit + 1)
 
phi :: Double
phi = (1 + sqrt 5)/2
log10phi = logBase 10 phi
halflog10five =(logBase 10 5)/2
 
sizeFib :: Integral a => a -> a
sizeFib p = ceiling $ (fromIntegral p)*log10phi - halflog10five
 
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
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[edit]

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) -> String
fromChrono (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 -> String
formatIntegral sep = reverse.intercalate sep.bagOf 3.reverse.show
 
formatAns :: Integer -> String
formatAns p = start ++ (startEnd 20 (sizeFib p) num)
where
start = "Fibonacci("++ (formatIntegral "_" p) ++ ") = "
num = fib p
 
seeFib :: Integer -> String
seeFib 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 -> String
startEnd ndigit len num | len <= ndigit = show num
startEnd ndigit len num | len <= 2*ndigit = let (us,vs) = genericSplitAt ndigit (show num) in us ++ " ... " ++ vs
startEnd ndigit len num = start ++ " ... " ++ end
where
end = show.rem num $ 10 ^ ndigit
start = show.quot num $ 10 ^ (len - ndigit + 1)
 
phi :: Double
phi = (1 + sqrt 5)/2
log10phi = logBase 10 phi
halflog10five =(logBase 10 5)/2
 
sizeFib :: Integral a => a -> a
sizeFib p = ceiling $ (fromIntegral p)*log10phi - halflog10five
 
fib :: Integer -> Integer
fib 0 = 0
fib 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 matrix
square :: 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 commute
mult :: 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 1
power m 1 = m
power 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[edit]

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[edit]

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)
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))
 
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[edit]

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[edit]

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

matrix exponentiation (2^16)[edit]

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 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
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[edit]

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

Sidef[edit]

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