Fibonacci matrix-exponentiation: Difference between revisions
m (→{{header|Go}}: Corrected typo) |
m (→{{header|Wren}}: Minor tidy) |
||
(76 intermediate revisions by 13 users not shown) | |||
Line 1: | Line 1: | ||
{{draft task|Matrices}} |
{{draft task|Matrices}} |
||
[[Category:Mathematics]] |
|||
The '''Fibonacci sequence''' defined with [[wp:Fibonacci_number#Matrix_form|matrix-exponentiation]]: |
The '''Fibonacci sequence''' defined with [[wp:Fibonacci_number#Matrix_form|matrix-exponentiation]]: |
||
::::: <math>\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}.</math> |
|||
<big><big> F<sub>0</sub> = 0 </big></big> |
|||
<big><big> F<sub>1</sub> = 1 </big></big> |
|||
:<math> |
|||
\begin{pmatrix} |
|||
F_{n+1} \\ |
|||
F_{n} |
|||
\end{pmatrix} |
|||
= |
|||
\begin{pmatrix} |
|||
1 & 1 \\ |
|||
1 & 0 |
|||
\end{pmatrix}^{n} |
|||
\begin{pmatrix} |
|||
F_{1} \\ |
|||
F_{0} |
|||
\end{pmatrix} |
|||
</math> |
|||
;Task: |
;Task: |
||
Write a |
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 |
|||
:::* 10,000,000 |
|||
Display only the 20 first digits and 20 last digits of each Fibonacci number. |
|||
Only display the first '''20''' decimal digits and the last '''20''' decimal digits of each Fibonacci number. |
|||
;Extra: |
;Extra: |
||
Generate |
Generate Fibonacci(2<sup>16 </sup>), Fibonacci(2<sup>32</sup>) and Fibonacci(2<sup>64</sup>) using the same method or another one. |
||
:<math>F_{n} = \cfrac{\varphi^n-(-\varphi)^{-n}}{\sqrt{5}}</math> |
|||
:<math>\varphi = \cfrac{1+\sqrt{5}}{2}</math> |
|||
;Related tasks: |
;Related tasks: |
||
* [[Fibonacci sequence]] |
* [[Fibonacci sequence]] |
||
* [[Matrix-exponentiation operator]] |
* [[Matrix-exponentiation operator]] |
||
<br><br> |
|||
=={{header|C sharp|C#}}== |
|||
===Matrix exponentiation=== |
|||
<syntaxhighlight lang="csharp">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); |
|||
} |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
===Matrix exponentiation=== |
|||
{{libheader|GMP(Go wrapper)}} |
|||
This uses matrix exponentiation to calculate the (2^32)nd Fibonacci |
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) |
I have not attempted to calculate the (2^64)th Fibonacci number which appears to be well out of reach using this approach. |
||
{{libheader|GMP(Go wrapper)}} |
|||
<lang go>package main |
|||
<syntaxhighlight lang="go">package main |
|||
import ( |
import ( |
||
Line 131: | Line 232: | ||
m = m.pow(n.Sub(n, one)) |
m = m.pow(n.Sub(n, one)) |
||
return m[0][0] |
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 |
|||
} |
} |
||
Line 136: | Line 246: | ||
start := time.Now() |
start := time.Now() |
||
n := new(big.Int) |
n := new(big.Int) |
||
for i := uint64(10); i <= 1e7; i *= 10 { |
|||
n.Lsh(one, 32) |
|||
n.SetUint64(i) |
|||
s := fibonacci(n).String() |
|||
fmt.Printf("The digits of the 2^32nd Fibonacci number (%d) are:\n", len(s)) |
|||
fmt.Printf(" |
fmt.Printf("The digits of the %sth Fibonacci number (%s) are:\n", |
||
commatize(i), commatize(uint64(len(s)))) |
|||
if len(s) > 20 { |
|||
fmt.Printf("\nTook %s\n", time.Since(start)) |
|||
fmt.Printf(" First 20 : %s\n", s[0:20]) |
|||
}</lang> |
|||
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)) |
|||
}</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
Timings are for an Intel Core i7 8565U machine, using Go 1. |
Timings are for an Intel Core i7 8565U machine, using Go 1.14 on Ubuntu 18.04: |
||
<pre> |
<pre> |
||
The digits of the |
The digits of the 10th Fibonacci number (2) are: |
||
All 2 : 55 |
|||
First twenty : 61999319689381859818 |
|||
Final twenty : 39623735538208076347 |
|||
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 |
Took 11m0.255916206s |
||
</pre> |
</pre> |
||
<br> |
<br> |
||
===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. |
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. |
||
<br> |
<br> |
||
{{trans|Julia}} |
{{trans|Julia}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 176: | Line 342: | ||
inner = func(n *big.Int) (*big.Int, *big.Int) { |
inner = func(n *big.Int) (*big.Int, *big.Int) { |
||
if n.Cmp(zero) == 0 { |
if n.Cmp(zero) == 0 { |
||
return |
return new(big.Int), big.NewInt(1) |
||
} |
} |
||
t, q, r := new(big.Int), new(big.Int), new(big.Int) |
t, q, r := new(big.Int), new(big.Int), new(big.Int) |
||
Line 209: | Line 375: | ||
} |
} |
||
return u.Mul(u, l) |
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 |
|||
} |
} |
||
Line 214: | Line 389: | ||
start := time.Now() |
start := time.Now() |
||
n := new(big.Int) |
n := new(big.Int) |
||
for i := uint64(10); i <= 1e7; i *= 10 { |
|||
n.Lsh(one, 32) |
|||
n.SetUint64(i) |
|||
s := lucas(n).String() |
|||
fmt.Printf("The digits of the 2^32nd Fibonacci number (%d) are:\n", len(s)) |
|||
fmt.Printf(" |
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)) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
As first version except time is now 5m4.13427997s. |
|||
</pre> |
|||
<br> |
|||
===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! |
|||
{{libheader|bigfloat}} |
|||
{{trans|Sidef}} |
|||
{{trans|Julia}} |
|||
<syntaxhighlight 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 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)) |
fmt.Printf("\nTook %s\n", time.Since(start)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Haskell}}== |
|||
===Matrix exponentiation=== |
|||
<syntaxhighlight lang="haskell">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</syntaxhighlight> |
|||
{{out}} |
|||
Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5: |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
===Matrix exponentiation - printing alternative === |
|||
<syntaxhighlight lang="haskell">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</syntaxhighlight> |
|||
{{out}} |
|||
Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5: |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
===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. |
|||
<syntaxhighlight lang="haskell">-- 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</syntaxhighlight> |
|||
{{out}} |
|||
Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5: |
|||
<pre>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</pre> |
|||
=={{header|Java}}== |
|||
Performed the task to use Matrix multiplication to compute Fibonacci numbers. |
|||
Implemented fib and fibMod. |
|||
<syntaxhighlight 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); |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren]]''' |
|||
<br> |
|||
'''Works with gojq, the Go implementation of jq''' (*) |
|||
(*) The C implementation of jq only has support for IEEE 754 64-bit numbers. |
|||
In the first part of this entry, the matrix method is used; the second part uses a recurrence relation for Fib(2^n) that is significantly faster, but still not practical for computing Fib(2^32). |
|||
====Matrix Exponentiation==== |
|||
<syntaxhighlight lang="jq">def mul($m1; $m2): |
|||
($m1|length) as $rows1 |
|||
| ($m1[0]|length) as $cols1 |
|||
| ($m2|length) as $rows2 |
|||
| ($m2[0]|length) as $cols2 |
|||
| if ($cols1 != $rows2) then "Matrices cannot be multiplied."| error |
|||
else reduce range(0; $rows1) as $i (null; |
|||
reduce range(0; $cols2) as $j (.; |
|||
.[$i][$j] = 0 |
|||
| reduce range(0; $rows2) as $k (.; |
|||
.[$i][$j] += $m1[$i][$k] * $m2[$k][$j]) |
|||
) ) |
|||
end ; |
|||
def identityMatrix: |
|||
. as $n |
|||
| [range(0; .) | 0] as $row |
|||
| [range(0; .) | $row] |
|||
| reduce range(0;$n) as $i (.; |
|||
.[$i][$i] = 1 ); |
|||
# . should be a square matrix and $n >= 0 |
|||
def pow( $n ): |
|||
. as $m |
|||
| ($m|length) as $le |
|||
| if $n < 0 then "Negative exponents not supported" | error |
|||
elif $n == 0 then $le|identityMatrix |
|||
elif $n == 1 then $m |
|||
else {pow : ($le | identityMatrix), |
|||
base : $m, |
|||
e : $n } |
|||
| until( .e <= 0.5; |
|||
# debug| |
|||
(.e % 2) as $temp |
|||
| if $temp == 1 then .pow = mul(.pow; .base) else . end |
|||
| .e /= 2 |
|||
| .base = mul(.base; .base) ) |
|||
| .pow |
|||
end; |
|||
def fibonacci: |
|||
. as $n |
|||
| if $n == 0 then 0 |
|||
else {m: [[1, 1], [1, 0]]} |
|||
| .m |= pow($n - 1) |
|||
| .m[0][0] |
|||
end; |
|||
def task1: |
|||
{ i: 10 } |
|||
| while ( .i <= 1e7; |
|||
.n = .i |
|||
| .s = (.n|fibonacci|tostring) |
|||
| .i *= 10) |
|||
| select(.s) |
|||
| "\nThe digits of the \(.n)th Fibonacci number (\(.s|length)) are:", |
|||
(if .s|length > 20 |
|||
then " First 20 : \(.s[0:20])", |
|||
( if (.s|length < 40) |
|||
then "\n Final \(.s|length-20): \(.s[20:])" |
|||
else " Final 20 : \(.s[-20:])" |
|||
end ) |
|||
else " All \(.s|length) : \(.s)" |
|||
end) ; |
|||
def task2: |
|||
pow(2;16) as $n |
|||
| (($n|fibonacci)|tostring) as $s |
|||
| "The digits of the 2^16th Fibonacci number \($s|length) are:", |
|||
" First 20 : \($s[0:20])", |
|||
" Final 20 : \($s[-20:])"; |
|||
task1, "", task2</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 1000th Fibonacci number (209) are: |
|||
First 20 : 43466557686937456435 |
|||
Final 20 : 76137795166849228875 |
|||
The digits of the 10000th Fibonacci number (2090) are: |
|||
First 20 : 33644764876431783266 |
|||
Final 20 : 66073310059947366875 |
|||
The digits of the 100000th Fibonacci number (20899) are: |
|||
First 20 : 25974069347221724166 |
|||
Final 20 : 49895374653428746875 |
|||
The digits of the 1000000th Fibonacci number (208988) are: |
|||
First 20 : 19532821287077577316 |
|||
Final 20 : 68996526838242546875 |
|||
The digits of the 2^16th Fibonacci number 13696 are: |
|||
First 20 : 73199214460290552832 |
|||
Final 20 : 97270190955307463227 |
|||
</pre> |
|||
====Fib(2^n)==== |
|||
<syntaxhighlight lang="jq"># Input: n |
|||
# Output: Fib(2^n) |
|||
def Fib2pN: |
|||
# in: [p, Fn-1, Fn] where n==2^p |
|||
# out: [2p, F(2n-1),F(2n)] |
|||
def fibonacci_recurrence: |
|||
def sq: .*.; |
|||
. as [$p, $fprev, $f] |
|||
| [1+$p, ($f|sq) + ($fprev|sq), (2*$fprev + $f)*$f]; |
|||
. as $n |
|||
| [0,0,1] |
|||
| until( .[0] >= $n; fibonacci_recurrence) |
|||
| .[2] ; |
|||
16, 32 |
|||
| . as $i |
|||
| Fib2pN |
|||
| tostring |
|||
| "The digits of the 2^\($i)th Fibonacci number (with string length \(length)) are:", |
|||
" First 20 : \(.[0:20])", |
|||
" Final 20 : \(.[-20:])", |
|||
""</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
Using gojq to compute Fib(2^32) using this method takes many hours. |
|||
<pre> |
<pre> |
||
The digits of the 2^ |
The digits of the 2^16th Fibonacci number (with string length 13696) are: |
||
First |
First 20 : 73199214460290552832 |
||
Final |
Final 20 : 97270190955307463227 |
||
The digits of the 2^32 Fibonacci number (with string length 897595080) are: |
|||
Took 5m2.722505927s |
|||
First 20 : 61999319689381859818 |
|||
Final 20 : 39623735538208076347 |
|||
</pre> |
</pre> |
||
Line 235: | Line 1,179: | ||
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 |
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. |
Lucas method is used as the alternative method. |
||
<syntaxhighlight 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 |
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) # not used here |
|||
binetfibonacci(n) = ((1+sqrt(big"5"))^n-(1-sqrt(big"5"))^n)/(sqrt(big"5")*big"2"^n) |
|||
# This is derived from the Schönhage-Strassen algorithm using Lucas numbers |
|||
# 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 lucasfibonacci(n) |
||
function inner(n) |
function inner(n) |
||
Line 251: | Line 1,220: | ||
u *= u |
u *= u |
||
v *= v |
v *= v |
||
return isodd(n) ? BigInt(u + v), BigInt(3 * v - 2 * (u - q)) : |
return isodd(n) ? (BigInt(u + v), BigInt(3 * v - 2 * (u - q))) : |
||
BigInt(2 * (v + q) - 3 * u), BigInt(u + v) |
(BigInt(2 * (v + q) - 3 * u), BigInt(u + v)) |
||
end |
end |
||
u, v = inner(n >> 1) |
u, v = inner(n >> 1) |
||
Line 262: | Line 1,231: | ||
return u * l |
return u * l |
||
end |
end |
||
m2s(bits) = string(matrixfibonacci(big"2"^bits)) |
m2s(bits) = string(matrixfibonacci(big"2"^bits)) |
||
l2s(bits) = string(lucasfibonacci(big"2"^bits)) |
l2s(bits) = string(lucasfibonacci(big"2"^bits)) |
||
firstlast(s) = (length(s) < |
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)) |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
N Matrix Lucas Mod |
|||
------------------------------------------------------------------------------------------------------------------------------------------------- |
|||
2^16 73199214460290552832...97270190955307463227 73199214460290552832...97270190955307463227 73199214460290552832...97270190955307463227 |
|||
2^32 61999319689381859818...39623735538208076347 61999319689381859818...39623735538208076347 61999319689381859818...39623735538208076347 |
|||
2^64 11175807536929528424...17529800348089840187 |
|||
</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|Wren}} |
|||
{{trans|Go}} |
|||
Using the Lucas method. |
|||
<syntaxhighlight lang="nim">import strformat, strutils, times |
|||
import bignum |
|||
let |
|||
One = newInt(1) |
|||
Two = newInt(2) |
|||
Three = newInt(3) |
|||
proc lucas(n: Int): Int = |
|||
proc inner(n: Int): (Int, Int) = |
|||
if n.isZero: return (newInt(0), newInt(1)) |
|||
var t = n shr 1 |
|||
var (u, v) = inner(t) |
|||
t = n and Two |
|||
let q = t - One |
|||
var r = newInt(0) |
|||
u *= u |
|||
v *= v |
|||
t = n and One |
|||
if t == One: |
|||
t = (u - q) * Two |
|||
r = v * Three |
|||
result = (u + v, r - t) |
|||
else: |
|||
t = u * Three |
|||
r = v + q |
|||
r *= Two |
|||
result = (r - t, u + v) |
|||
var t = n shr 1 |
|||
let (u, v) = inner(t) |
|||
let l = v * Two - u |
|||
t = n and One |
|||
if t == One: |
|||
let q = n and Two - One |
|||
return v * l + q |
|||
return u * l |
|||
let start = now() |
|||
var n: Int |
|||
var i = 10 |
|||
while i <= 10_000_000: |
|||
n = newInt(i) |
|||
let s = $lucas(n) |
|||
echo &"The digits of the {($i).insertSep}th Fibonacci number ({($s.len).insertSep}) are:" |
|||
if s.len > 20: |
|||
echo &" First 20 : {s[0..19]}" |
|||
if s.len < 40: |
|||
echo &" Final {s.len-20:<2} : {s[20..^1]}" |
|||
else: |
|||
echo &" Final 20 : {s[^20..^1]}" |
|||
else: |
|||
echo &" All {s.len:<2} : {s}" |
|||
echo() |
|||
i *= 10 |
|||
for e in [culong 16, 32]: |
|||
n = One shl e |
|||
let s = $lucas(n) |
|||
echo &"The digits of the 2^{e}th Fibonacci number ({($s.len).insertSep}) are:" |
|||
echo &" First 20 : {s[0..19]}" |
|||
echo &" Final 20 : {s[^20..^1]}" |
|||
echo() |
|||
echo &"Took {now() - start}"</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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^16th 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 6 minutes, 18 seconds, 643 milliseconds, 622 microseconds, and 965 nanoseconds</pre> |
|||
=={{header|Perl}}== |
|||
{{trans|Sidef}} |
|||
<syntaxhighlight 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; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|Phix}}== |
|||
{{libheader|Phix/mpfr}} {{trans|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.<br> |
|||
(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.) |
|||
<!--<syntaxhighlight lang="phix">--> |
|||
<span style="color: #008080;">without</span> <span style="color: #008080;">javascript_semantics</span> <span style="color: #000080;font-style:italic;">-- (no mpfr_log(), mpfr_exp() under pwa/p2js)</span> |
|||
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.0"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (mpfr_set_default_prec[ision] has been renamed)</span> |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
<span style="color: #7060A8;">mpfr_set_default_precision</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">40</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">PHI</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1.25</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">IHP</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(),</span> |
|||
<span style="color: #000000;">HLF</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">L10</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">IHP</span><span style="color: #0000FF;">,</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">HLF</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">HLF</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_neg</span><span style="color: #0000FF;">(</span><span style="color: #000000;">IHP</span><span style="color: #0000FF;">,</span><span style="color: #000000;">IHP</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">L10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">L10</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">binet_approx</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpfr</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #000000;">mpfr_ui_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (n as in 2^n here)</span> |
|||
<span style="color: #000000;">mpfr_log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">PHI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">IHP</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_free</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">nth_fib_first_k_digits</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">mpfr</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">h</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">binet_approx</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">L10</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_add_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">g</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">,</span><span style="color: #000000;">L10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">g</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">g</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_ui_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">h</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpfr_floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%%%d.0Rf"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">mpfr_sprintf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">cache</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">new_dict</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #000080;font-style:italic;">-- key of {n,m} means 2^n+m (where n and m are integers, n>=0, m always 0 or -1)</span> |
|||
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">({</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">),</span><span style="color: #000000;">cache</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- aka fib(0):=0</span> |
|||
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">({</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span><span style="color: #000000;">cache</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- fib(1):=1</span> |
|||
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">({</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span><span style="color: #000000;">cache</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- fib(2):=1</span> |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">mpz_fibmod</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">key</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">getd_index</span><span style="color: #0000FF;">(</span><span style="color: #000000;">key</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cache</span><span style="color: #0000FF;">)!=</span><span style="color: #004600;">NULL</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">getd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">key</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cache</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">else</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">n</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">1</span> |
|||
<span style="color: #000000;">mpz_fibmod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpz_fibmod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">=-</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000080;font-style:italic;">-- fib(2n-1) = fib(n)^2 + fib(n-1)^2</span> |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">else</span> |
|||
<span style="color: #000080;font-style:italic;">-- fib(2n) = (2*fib(n-1)+fib(n))*fib(n)</span> |
|||
<span style="color: #7060A8;">mpz_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #7060A8;">mpz_mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">setd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">key</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cache</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">nth_fib_last_k_digits</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_ui_pow_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">mpz_fibmod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">16</span><span style="color: #0000FF;">,</span><span style="color: #000000;">32</span><span style="color: #0000FF;">,</span><span style="color: #000000;">64</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">first20</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nth_fib_first_k_digits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">last20</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nth_fib_last_k_digits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"F(2^%d) = %s ... %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">first20</span><span style="color: #0000FF;">,</span><span style="color: #000000;">last20</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
|||
<pre> |
|||
F(2^16) = 73199214460290552832 ... 97270190955307463227 |
|||
F(2^32) = 61999319689381859818 ... 39623735538208076347 |
|||
F(2^64) = 11175807536929528424 ... 17529800348089840187 |
|||
</pre> |
|||
===matrix exponentiation (2^16)=== |
|||
Somewhat closer to the original specification, but certainly not recommended for 2^32, let alone 2^64...<br> |
|||
Contains copies of routines from [[Matrix-exponentiation_operator#Phix]], but modified to use gmp. |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">ZERO</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">ONE</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">identity</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ZERO</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ONE</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])!=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"matrices cannot be multiplied"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ZERO</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])),</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">cij</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cij</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cij</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cij</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">cij</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">c</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">matrix_exponent</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">!=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"not a square matrix"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">crash</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"negative exponents not supported"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">identity</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">and_bits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #000080;font-style:italic;">-- (avoid unnecessary matrix_mul)</span> |
|||
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">fibonacci</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">ONE</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ONE</span><span style="color: #0000FF;">},</span> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">ONE</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ZERO</span><span style="color: #0000FF;">}}</span> |
|||
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_exponent</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fibonacci</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">16</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fn</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"fibonnaci(2^16) = %s [%s]\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">shorten</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">),</span><span style="color: #000000;">e</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()=</span><span style="color: #004600;">JS</span><span style="color: #0000FF;">?</span><span style="color: #000000;">6</span><span style="color: #0000FF;">:</span><span style="color: #000000;">7</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">fibonacci</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fn</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span> |
|||
<span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">></span><span style="color: #000000;">0.1</span><span style="color: #0000FF;">?</span><span style="color: #008000;">" ["</span><span style="color: #0000FF;">&</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)&</span><span style="color: #008000;">"]"</span><span style="color: #0000FF;">:</span><span style="color: #008000;">""</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"fibonnaci(%,d) = %s%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">shorten</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">),</span><span style="color: #000000;">e</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
|||
<pre> |
|||
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] |
|||
</pre> |
|||
Note that under pwa/p2js 10^6 takes 11.4s so we'll stop there...<br> |
|||
Change that loop to 8 and a 9 year old 3.3GHz i3 also eventually gets: |
|||
<pre> |
|||
fibonnaci(100,000,000) = 47371034734563369625...06082642167760546875 (20,898,764 digits) [14 minutes and 24s] |
|||
</pre> |
|||
Clearly 2^32 (897 million digits, apparently) is a tad out of bounds, let alone 2^64. |
|||
=={{header|Python}}== |
|||
With fancy custom integer classes that are absolutely useless except for this task. |
|||
<syntaxhighlight lang="python">class Head(): |
|||
def __init__(self, lo, hi=None, shift=0): |
|||
if hi is None: hi = lo |
|||
d = hi - lo |
|||
ds, ls, hs = str(d), str(lo), str(hi) |
|||
if d and len(ls) > len(ds): |
|||
assert(len(ls) - len(ds) + 1 > 21) |
|||
lo = int(str(lo)[:len(ls) - len(ds) + 1]) |
|||
hi = int(str(hi)[:len(hs) - len(ds) + 1]) + 1 |
|||
shift += len(ds) - 1 |
|||
elif len(ls) > 100: |
|||
lo = int(str(ls)[:100]) |
|||
hi = lo + 1 |
|||
shift = len(ls) - 100 |
|||
self.lo, self.hi, self.shift = lo, hi, shift |
|||
def __mul__(self, other): |
|||
lo = self.lo*other.lo |
|||
hi = self.hi*other.hi |
|||
shift = self.shift + other.shift |
|||
return Head(lo, hi, shift) |
|||
def __add__(self, other): |
|||
if self.shift < other.shift: |
|||
return other + self |
|||
sh = self.shift - other.shift |
|||
if sh >= len(str(other.hi)): |
|||
return Head(self.lo, self.hi, self.shift) |
|||
ls = str(other.lo) |
|||
hs = str(other.hi) |
|||
lo = self.lo + int(ls[:len(ls)-sh]) |
|||
hi = self.hi + int(hs[:len(hs)-sh]) |
|||
return Head(lo, hi, self.shift) |
|||
def __repr__(self): |
|||
return str(self.hi)[:20] |
|||
class Tail(): |
|||
def __init__(self, v): |
|||
self.v = int(f'{v:020d}'[-20:]) |
|||
def __add__(self, other): |
|||
return Tail(self.v + other.v) |
|||
def __mul__(self, other): |
|||
return Tail(self.v*other.v) |
|||
def __repr__(self): |
|||
return f'{self.v:020d}'[-20:] |
|||
def mul(a, b): |
|||
return a[0]*b[0] + a[1]*b[1], a[0]*b[1] + a[1]*b[2], a[1]*b[1] + a[2]*b[2] |
|||
def fibo(n, cls): |
|||
n -= 1 |
|||
zero, one = cls(0), cls(1) |
|||
m = (one, one, zero) |
|||
e = (one, zero, one) |
|||
while n: |
|||
if n&1: e = mul(m, e) |
|||
m = mul(m, m) |
|||
n >>= 1 |
|||
return f'{e[0]}' |
|||
for i in range(2, 10): |
|||
n = 10**i |
|||
print(f'10^{i} :', fibo(n, Head), '...', fibo(n, Tail)) |
|||
for i in range(3, 8): |
|||
n = 2**i |
|||
s = f'2^{n}' |
|||
print(f'{s:5s}:', fibo(2**n, Head), '...', fibo(2**n, Tail))</syntaxhighlight> |
|||
{{out}} |
|||
<pre>10^2 : 35422484817926191507 ... 54224848179261915075 |
|||
10^3 : 43466557686937456435 ... 76137795166849228875 |
|||
10^4 : 33644764876431783266 ... 66073310059947366875 |
|||
10^5 : 25974069347221724166 ... 49895374653428746875 |
|||
10^6 : 19532821287077577316 ... 68996526838242546875 |
|||
10^7 : 11298343782253997603 ... 86998673686380546875 |
|||
10^8 : 47371034734563369625 ... 06082642167760546875 |
|||
10^9 : 79523178745546834678 ... 03172326981560546875 |
|||
2^8 : 14169381771405651323 ... 19657707794958199867 |
|||
2^16 : 73199214460290552832 ... 97270190955307463227 |
|||
2^32 : 61999319689381859818 ... 39623735538208076347 |
|||
2^64 : 11175807536929528424 ... 17529800348089840187 |
|||
2^128: 15262728879740471565 ... 17229324095882654267</pre> |
|||
=={{header|Ruby}}== |
|||
===Matrix exponentiation by Ruby's fast exponentiation operator duck-typing applied to Ruby's built-in Integer Class=== |
|||
<pre> |
|||
require 'matrix' |
|||
start_time = Time.now |
|||
[0,1,2,3,4,10,100,256, 1_000, 1024, 10_000, 2 ** 16, 100_000, 1_000_000,10_000_000 ].each {|n| |
|||
## |
|||
fib_Num=(Matrix[[0,1],[1,1]] ** (n))[0,1] ## Matrix exponentiation |
|||
## |
|||
fib_Str= fib_Num.to_s() |
|||
if fib_Str.length <= 21 |
|||
p ["Fibonacci(#{n})",fib_Str.length.to_s + ' digits' , fib_Str] |
|||
else |
|||
p ["Fibonacci(#{n})",fib_Str.length.to_s + ' digits' , fib_Str.slice(0,20) + " ... " + fib_Str.slice(-20,20)] |
|||
end |
|||
} |
|||
puts "Took #{Time.now - start_time}s" |
|||
</pre> |
|||
{{out}} |
|||
<pre> |
|||
["Fibonacci(0)", "1 digits", "0"] |
|||
["Fibonacci(1)", "1 digits", "1"] |
|||
["Fibonacci(2)", "1 digits", "1"] |
|||
["Fibonacci(3)", "1 digits", "2"] |
|||
["Fibonacci(4)", "1 digits", "3"] |
|||
["Fibonacci(10)", "2 digits", "55"] |
|||
["Fibonacci(100)", "21 digits", "354224848179261915075"] |
|||
["Fibonacci(256)", "54 digits", "14169381771405651323 ... 19657707794958199867"] |
|||
["Fibonacci(1000)", "209 digits", "43466557686937456435 ... 76137795166849228875"] |
|||
["Fibonacci(1024)", "214 digits", "45066996336778198131 ... 04103631553925405243"] |
|||
["Fibonacci(10000)", "2090 digits", "33644764876431783266 ... 66073310059947366875"] |
|||
["Fibonacci(65536)", "13696 digits", "73199214460290552832 ... 97270190955307463227"] |
|||
["Fibonacci(100000)", "20899 digits", "25974069347221724166 ... 49895374653428746875"] |
|||
["Fibonacci(1000000)", "208988 digits", "19532821287077577316 ... 68996526838242546875"] |
|||
["Fibonacci(10000000)", "2089877 digits", "11298343782253997603 ... 86998673686380546875"] |
|||
Took 1.027948065s |
|||
</pre> |
|||
===Matrix exponentiation by Ruby's Matlix#exponentiation duck-typeing with Head-tail-BigNumber class=== |
|||
Here, the HeadTailBignum class holds the digits of the head part in Ruby's bigDecimal class and tail part in the BigInteger class. |
|||
and HeadTailBignum defines only the add and multply operators and the coerce method. |
|||
Matlix exponentiation operator is fast and can apply duck-typing to HeadTailBignum. |
|||
require 'matrix' |
|||
require 'bigdecimal' |
|||
class HeadTailBignum < Numeric |
|||
println("N", " "^23, "Matrix", " "^40, "Lucas\n", "-"^97) |
|||
attr_accessor :hi, :lo |
|||
println("2^32 ", rpad(firstlast(m2s(32)), 45), rpad(firstlast(l2s(32)), 45)) |
|||
@@degitsLimit = 20 |
|||
# println("2^64 ", rpad(firstlast(m2s(64)), 45), rpad(firstlast(l2s(64)), 45)) |
|||
@@hiDegitsLimit = (@@degitsLimit + 1) * 2 |
|||
</lang>{{out}} |
|||
@@loDegitsBase = 10 ** @@degitsLimit |
|||
BigDecimal::limit(@@hiDegitsLimit) |
|||
def initialize(other, loValue = nil) |
|||
if other.kind_of?(BigDecimal) && loValue.kind_of?(Integer) |
|||
@hi = other |
|||
@lo = loValue % @@loDegitsBase |
|||
elsif other.kind_of?(HeadTailBignum) |
|||
@hi = other.hi |
|||
@lo = other.lo |
|||
elsif other.kind_of?(Integer) |
|||
@hi = BigDecimal(other) |
|||
@lo = other % @@loDegitsBase |
|||
else |
|||
raise StandardError.new("HeadTailBignum initialize Type Error (" + other.class.to_s + "," + loValue.class.to_s + ")") |
|||
end |
|||
end |
|||
def clone |
|||
HeadTailBignum(self.hi, self.lo) |
|||
end |
|||
def integer?() |
|||
true |
|||
end |
|||
def coerce(other) |
|||
if other.kind_of?(Integer) |
|||
return HeadTailBignum.new(other), self |
|||
else |
|||
super |
|||
end |
|||
end |
|||
def +(other) |
|||
if other.kind_of?(Integer) |
|||
return HeadTailBignum.new(self.hi + other, self.lo + other) |
|||
elsif other.kind_of?(HeadTailBignum) |
|||
return HeadTailBignum.new(self.hi + other.hi , self.lo + other.lo) |
|||
else |
|||
raise StandardError.new("HeadTailBignum add Type Error (" + other.class.to_s + ")") |
|||
end |
|||
end |
|||
def *(other) |
|||
if other.kind_of?(Integer) |
|||
return HeadTailBignum.new(self.hi * other, self.lo * other) |
|||
elsif other.kind_of?(HeadTailBignum) |
|||
return HeadTailBignum.new(self.hi * other.hi , self.lo * other.lo) |
|||
else |
|||
raise StandardError.new("HeadTailBignum mulply Type Error (" + other.class.to_s + ")") |
|||
end |
|||
end |
|||
def exponent |
|||
@hi.exponent |
|||
end |
|||
def to_s |
|||
if @hi < @@loDegitsBase |
|||
return @lo.inspect |
|||
else |
|||
return @hi.to_s("E").slice(2,@@degitsLimit ) + " ... " + @lo.to_s |
|||
end |
|||
end |
|||
end |
|||
start_time = Time.now |
|||
[8,16,32,64].each {|h| |
|||
n = (2**h) |
|||
fib_Num=(Matrix[[ HeadTailBignum.new(0),1],[1,1]] ** (n))[0,1] |
|||
puts "Fibonacci(2^#{h.to_s}) = #{fib_Num.to_s} are #{fib_Num.exponent} digits" |
|||
} |
|||
puts "Took #{Time.now - start_time}s" |
|||
{{out}} |
|||
<pre> |
<pre> |
||
Fibonacci(2^8) = 14169381771405651323 ... 19657707794958199867 are 54 digits |
|||
N Matrix Lucas |
|||
Fibonacci(2^16) = 73199214460290552832 ... 97270190955307463227 are 13696 digits |
|||
------------------------------------------------------------------------------------------------- |
|||
2^32 61999319689381859818...39623735538208076347 |
Fibonacci(2^32) = 61999319689381859818 ... 39623735538208076347 are 897595080 digits |
||
Fibonacci(2^64) = 11175807536929528424 ... 17529800348089840187 are 3855141514259838964 digits |
|||
Took 0.009357194s |
|||
</pre> |
|||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
Following the general approach of Sidef, and relying on Perl for <code>fibmod</code> function. Borrowed/adapted routines from [http://rosettacode.org/wiki/Ramanujan%27s_constant Ramanujan's_constant] task to allow <code>FatRat</code> calculations throughout. Does not quite meet task spec, as n^64 results not working yet. |
|||
<syntaxhighlight lang="raku" line>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 |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>F(2^16) = 73199214460290552832 ... 97270190955307463227 |
|||
F(2^32) = 61999319689381859818 ... 39623735538208076347 |
|||
F(2^16) = 73199214460290552832 ... 97270190955307463227 |
|||
F(2^32) = 61999319689381859818 ... 39623735538208076347</pre> |
|||
=={{header|Sidef}}== |
|||
Computing the n-th Fibonacci number, using matrix-exponentiation (this function is also built-in): |
|||
<syntaxhighlight 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]</syntaxhighlight> |
|||
First and last 20 digits of the n-th Fibonacci number: |
|||
<syntaxhighlight 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}" |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
F(2^16) = 73199214460290552832 ... 97270190955307463227 |
|||
F(2^32) = 61999319689381859818 ... 39623735538208076347 |
|||
</pre> |
|||
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: |
|||
<syntaxhighlight 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}" |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
F(2^16) = 73199214460290552832 ... 97270190955307463227 |
|||
F(2^32) = 61999319689381859818 ... 39623735538208076347 |
|||
F(2^64) = 11175807536929528424 ... 17529800348089840187 |
|||
</pre> |
|||
=={{header|Wren}}== |
|||
===Matrix exponentiation=== |
|||
{{trans|Go}} |
|||
{{libheader|Wren-big}} |
|||
{{libheader|Wren-fmt}} |
|||
A tough task for Wren which takes just under 3 minutes to process even the 1 millionth Fibonacci number. Judging by the times for the compiled, statically typed languages using GMP, the 10 millionth would likely take north of 4 hours so I haven't attempted it. |
|||
<syntaxhighlight lang="wren">import "./big" for BigInt |
|||
import "./fmt" for Fmt |
|||
var mul = Fn.new { |m1, m2| |
|||
var rows1 = m1.count |
|||
var cols1 = m1[0].count |
|||
var rows2 = m2.count |
|||
var cols2 = m2[0].count |
|||
if (cols1 != rows2) Fiber.abort("Matrices cannot be multiplied.") |
|||
var result = List.filled(rows1, null) |
|||
for (i in 0...rows1) { |
|||
result[i] = List.filled(cols2, null) |
|||
for (j in 0...cols2) { |
|||
result[i][j] = BigInt.zero |
|||
for (k in 0...rows2) { |
|||
var temp = m1[i][k] * m2[k][j] |
|||
result[i][j] = result[i][j] + temp |
|||
} |
|||
} |
|||
} |
|||
return result |
|||
} |
|||
var identityMatrix = Fn.new { |n| |
|||
if (n < 1) Fiber.abort("Size of identity matrix can't be less than 1") |
|||
var ident = List.filled(n, 0) |
|||
for (i in 0...n) { |
|||
ident[i] = List.filled(n, null) |
|||
for(j in 0...n) ident[i][j] = (i != j) ? BigInt.zero : BigInt.one |
|||
} |
|||
return ident |
|||
} |
|||
var pow = Fn.new { |m, n| |
|||
var le = m.count |
|||
if (le != m[0].count) Fiber.abort("Not a square matrix") |
|||
if (n < 0) Fiber.abort("Negative exponents not supported") |
|||
if (n == 0) return identityMatrix.call(le) |
|||
if (n == 1) return m |
|||
var pow = identityMatrix.call(le) |
|||
var base = m |
|||
var e = n |
|||
while (e > 0) { |
|||
var temp = e & BigInt.one |
|||
if (temp == BigInt.one) pow = mul.call(pow, base) |
|||
e = e >> 1 |
|||
base = mul.call(base, base) |
|||
} |
|||
return pow |
|||
} |
|||
var fibonacci = Fn.new { |n| |
|||
if (n == 0) return BigInt.zero |
|||
var m = [[BigInt.one, BigInt.one], [BigInt.one, BigInt.zero]] |
|||
m = pow.call(m, n - 1) |
|||
return m[0][0] |
|||
} |
|||
var n = BigInt.zero |
|||
var i = 10 |
|||
while (i <= 1e6) { |
|||
n = BigInt.new(i) |
|||
var s = fibonacci.call(n).toString |
|||
Fmt.print("The digits of the $,sth Fibonacci number ($,s) are:", i, s.count) |
|||
if (s.count > 20) { |
|||
Fmt.print(" First 20 : $s", s[0...20]) |
|||
if (s.count < 40) { |
|||
Fmt.print(" Final $-2d : $s", s.count-20, s[20..-1]) |
|||
} else { |
|||
Fmt.print(" Final 20 : $s", s[s.count-20..-1]) |
|||
} |
|||
} else { |
|||
Fmt.print(" All $-2d : $s", s.count, s) |
|||
} |
|||
System.print() |
|||
i = i * 10 |
|||
} |
|||
n = BigInt.one << 16 |
|||
var s = fibonacci.call(n).toString |
|||
Fmt.print("The digits of the 2^16th Fibonacci number ($,s) are:", s.count) |
|||
Fmt.print(" First 20 : $s", s[0...20]) |
|||
Fmt.print(" Final 20 : $s", s[s.count-20..-1])</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 2^16th Fibonacci number (13,696) are: |
|||
First 20 : 73199214460290552832 |
|||
Final 20 : 97270190955307463227 |
|||
</pre> |
|||
<br> |
|||
===Lucas method=== |
|||
{{trans|Go}} |
|||
Much quicker than the matrix exponentiation method taking about 23 seconds to process the 1 millionth Fibonacci number and around 32 minutes to reach the 10 millionth. |
|||
Apart from the 2^16th number, the extra credit is still out of reach using this approach. |
|||
<syntaxhighlight lang="wren">import "./big" for BigInt |
|||
import "./fmt" for Fmt |
|||
var lucas = Fn.new { |n| |
|||
var inner // recursive function |
|||
inner = Fn.new { |n| |
|||
if (n == BigInt.zero) return [BigInt.zero, BigInt.one] |
|||
var t = n >> 1 |
|||
var res = inner.call(t) |
|||
var u = res[0] |
|||
var v = res[1] |
|||
t = n & BigInt.two |
|||
var q = t - BigInt.one |
|||
var r = BigInt.zero |
|||
u = u.square |
|||
v = v.square |
|||
t = n & BigInt.one |
|||
if (t == BigInt.one) { |
|||
t = u - q |
|||
t = BigInt.two * t |
|||
r = BigInt.three * v |
|||
return [u + v, r - t] |
|||
} else { |
|||
t = BigInt.three * u |
|||
r = v + q |
|||
r = BigInt.two * r |
|||
return [r - t, u + v] |
|||
} |
|||
} |
|||
var t = n >> 1 |
|||
var res = inner.call(t) |
|||
var u = res[0] |
|||
var v = res[1] |
|||
var l = BigInt.two * v |
|||
l = l - u // Lucas function |
|||
t = n & BigInt.one |
|||
if (t == BigInt.one) { |
|||
var q = n & BigInt.two |
|||
q = q - BigInt.one |
|||
t = v * l |
|||
return t + q |
|||
} |
|||
return u * l |
|||
} |
|||
var n = BigInt.zero |
|||
var i = 10 |
|||
while (i <= 1e7) { |
|||
n = BigInt.new(i) |
|||
var s = lucas.call(n).toString |
|||
Fmt.print("The digits of the $,sth Fibonacci number ($,s) are:", i, s.count) |
|||
if (s.count > 20) { |
|||
Fmt.print(" First 20 : $s", s[0...20]) |
|||
if (s.count < 40) { |
|||
Fmt.print(" Final $-2d : $s", s.count-20, s[20..-1]) |
|||
} else { |
|||
Fmt.print(" Final 20 : $s", s[s.count-20..-1]) |
|||
} |
|||
} else { |
|||
Fmt.print(" All $-2d : $s", s.count, s) |
|||
} |
|||
System.print() |
|||
i = i * 10 |
|||
} |
|||
n = BigInt.one << 16 |
|||
var s = lucas.call(n).toString |
|||
Fmt.print("The digits of the 2^16th Fibonacci number ($,s) are:", s.count) |
|||
Fmt.print(" First 20 : $s", s[0...20]) |
|||
Fmt.print(" Final 20 : $s", s[s.count-20..-1])</syntaxhighlight> |
|||
{{out}} |
|||
As matrix exponentiation method, plus: |
|||
<pre> |
|||
The digits of the 10,000,000th Fibonacci number (2,089,877) are: |
|||
First 20 : 11298343782253997603 |
|||
Final 20 : 86998673686380546875 |
|||
</pre> |
|||
<br> |
|||
===Fibmod method (embedded)=== |
|||
{{trans|Go}} |
|||
{{libheader|Wren-gmp}} |
|||
Ridiculously fast (under 5ms) thanks to the stuff borrowed from Sidef and Julia combined with the speed of GMP. |
|||
<syntaxhighlight lang="wren">import "./gmp" for Mpz, Mpf |
|||
import "./fmt" for Fmt |
|||
var nd = 20 // number of digits to be displayed at each end |
|||
var pr = 128 // precision to be used |
|||
var one = Mpz.one |
|||
var two = Mpz.two |
|||
var fibmod = Fn.new { |n, nmod| |
|||
if (n < two) return n |
|||
var fibmods = {} |
|||
var f // recursive closure |
|||
f = Fn.new { |n| |
|||
if (n < two) return one |
|||
var ns = n.toString |
|||
var v = fibmods[ns] |
|||
if (v) return v |
|||
v = Mpz.zero |
|||
var k = n / two |
|||
var t = n & one |
|||
var u = Mpz.zero |
|||
if (t != one) { |
|||
t.set(f.call(k)).square |
|||
v.sub(k, one) |
|||
u.set(f.call(v)).square |
|||
} else { |
|||
t.set(f.call(k)) |
|||
v.add(k, one) |
|||
v.set(f.call(v)) |
|||
u.sub(k, one) |
|||
u.set(f.call(u)).mul(t) |
|||
t.mul(v) |
|||
} |
|||
t.add(u) |
|||
fibmods[ns] = t.rem(nmod) |
|||
return fibmods[ns] |
|||
} |
|||
var w = n - one |
|||
return f.call(w) |
|||
} |
|||
var binetApprox = Fn.new { |n| |
|||
var phi = Mpf.from(0.5, pr) |
|||
var ihp = phi.copy() |
|||
var root = Mpf.from(1.25, pr).sqrt |
|||
phi.add(root) |
|||
ihp.sub(root, ihp).neg |
|||
ihp.sub(phi, ihp).log |
|||
phi.log |
|||
var nn = Mpf.from(n, pr) |
|||
return phi.mul(nn).sub(ihp) |
|||
} |
|||
var firstFibDigits = Fn.new { |n, k| |
|||
var f = binetApprox.call(n) |
|||
var g = Mpf.new(pr) |
|||
g.div(f, Mpf.ln10(pr)).inc |
|||
g.floor.mul(Mpf.ln10(pr)) |
|||
f.sub(g).exp |
|||
var p = Mpz.from(10).pow(k) |
|||
g.set(p) |
|||
f.mul(g) |
|||
return f.floor.toString[0...k] |
|||
} |
|||
var lastFibDigits = Fn.new { |n, k| |
|||
var p = Mpz.from(10).pow(k) |
|||
return fibmod.call(n, p).toString[0...k] |
|||
} |
|||
var start = System.clock |
|||
var n = Mpz.zero |
|||
var i = 10 |
|||
while (i <= 1e7) { |
|||
n.set(i) |
|||
var nn = Mpz.from(i) |
|||
Fmt.print("\nThe digits of the $,r Fibonacci number are:", i) |
|||
var nd2 = nd |
|||
var nd3 = 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 |
|||
} |
|||
var s1 = firstFibDigits.call(n, nd2) |
|||
if (s1.count < 20) { |
|||
Fmt.print(" All $-2d : $s", s1.count, s1) |
|||
} else { |
|||
Fmt.print(" First 20 : $s", s1) |
|||
var s2 = lastFibDigits.call(nn, nd3) |
|||
if (s2.count < 20) { |
|||
Fmt.print(" Final $-2d : $s", s2.count, s2) |
|||
} else { |
|||
Fmt.print(" Final 20 : $s", s2) |
|||
} |
|||
} |
|||
i = i * 10 |
|||
} |
|||
var ord = ["th", "nd", "th"] |
|||
i = 0 |
|||
for (p in [16, 32, 64]) { |
|||
n.lsh(one, p) |
|||
var nn = one << p |
|||
Fmt.print("\nThe digits of the 2^$d$s Fibonacci number are:", p, ord[i]) |
|||
Fmt.print(" First $d : $s", nd, firstFibDigits.call(n, nd)) |
|||
Fmt.print(" Final $d : $s", nd, lastFibDigits.call(nn, nd)) |
|||
i = i + 1 |
|||
} |
|||
System.print("\nTook %(System.clock-start) seconds.")</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 0.004516 seconds. |
|||
</pre> |
</pre> |
Latest revision as of 12:12, 4 December 2023
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
- 10,000,000
Only display the first 20 decimal digits and the last 20 decimal digits of each Fibonacci number.
- Extra
Generate Fibonacci(216 ), Fibonacci(232) and Fibonacci(264) using the same method or another one.
- Related tasks
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.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
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.
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!
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
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) -> 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
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
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
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
jq
Adapted from #Wren
Works with gojq, the Go implementation of jq (*)
(*) The C implementation of jq only has support for IEEE 754 64-bit numbers.
In the first part of this entry, the matrix method is used; the second part uses a recurrence relation for Fib(2^n) that is significantly faster, but still not practical for computing Fib(2^32).
Matrix Exponentiation
def mul($m1; $m2):
($m1|length) as $rows1
| ($m1[0]|length) as $cols1
| ($m2|length) as $rows2
| ($m2[0]|length) as $cols2
| if ($cols1 != $rows2) then "Matrices cannot be multiplied."| error
else reduce range(0; $rows1) as $i (null;
reduce range(0; $cols2) as $j (.;
.[$i][$j] = 0
| reduce range(0; $rows2) as $k (.;
.[$i][$j] += $m1[$i][$k] * $m2[$k][$j])
) )
end ;
def identityMatrix:
. as $n
| [range(0; .) | 0] as $row
| [range(0; .) | $row]
| reduce range(0;$n) as $i (.;
.[$i][$i] = 1 );
# . should be a square matrix and $n >= 0
def pow( $n ):
. as $m
| ($m|length) as $le
| if $n < 0 then "Negative exponents not supported" | error
elif $n == 0 then $le|identityMatrix
elif $n == 1 then $m
else {pow : ($le | identityMatrix),
base : $m,
e : $n }
| until( .e <= 0.5;
# debug|
(.e % 2) as $temp
| if $temp == 1 then .pow = mul(.pow; .base) else . end
| .e /= 2
| .base = mul(.base; .base) )
| .pow
end;
def fibonacci:
. as $n
| if $n == 0 then 0
else {m: [[1, 1], [1, 0]]}
| .m |= pow($n - 1)
| .m[0][0]
end;
def task1:
{ i: 10 }
| while ( .i <= 1e7;
.n = .i
| .s = (.n|fibonacci|tostring)
| .i *= 10)
| select(.s)
| "\nThe digits of the \(.n)th Fibonacci number (\(.s|length)) are:",
(if .s|length > 20
then " First 20 : \(.s[0:20])",
( if (.s|length < 40)
then "\n Final \(.s|length-20): \(.s[20:])"
else " Final 20 : \(.s[-20:])"
end )
else " All \(.s|length) : \(.s)"
end) ;
def task2:
pow(2;16) as $n
| (($n|fibonacci)|tostring) as $s
| "The digits of the 2^16th Fibonacci number \($s|length) are:",
" First 20 : \($s[0:20])",
" Final 20 : \($s[-20:])";
task1, "", task2
- Output:
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 1000th Fibonacci number (209) are: First 20 : 43466557686937456435 Final 20 : 76137795166849228875 The digits of the 10000th Fibonacci number (2090) are: First 20 : 33644764876431783266 Final 20 : 66073310059947366875 The digits of the 100000th Fibonacci number (20899) are: First 20 : 25974069347221724166 Final 20 : 49895374653428746875 The digits of the 1000000th Fibonacci number (208988) are: First 20 : 19532821287077577316 Final 20 : 68996526838242546875 The digits of the 2^16th Fibonacci number 13696 are: First 20 : 73199214460290552832 Final 20 : 97270190955307463227
Fib(2^n)
# Input: n
# Output: Fib(2^n)
def Fib2pN:
# in: [p, Fn-1, Fn] where n==2^p
# out: [2p, F(2n-1),F(2n)]
def fibonacci_recurrence:
def sq: .*.;
. as [$p, $fprev, $f]
| [1+$p, ($f|sq) + ($fprev|sq), (2*$fprev + $f)*$f];
. as $n
| [0,0,1]
| until( .[0] >= $n; fibonacci_recurrence)
| .[2] ;
16, 32
| . as $i
| Fib2pN
| tostring
| "The digits of the 2^\($i)th Fibonacci number (with string length \(length)) are:",
" First 20 : \(.[0:20])",
" Final 20 : \(.[-20:])",
""
- Output:
Using gojq to compute Fib(2^32) using this method takes many hours.
The digits of the 2^16th Fibonacci number (with string length 13696) are: First 20 : 73199214460290552832 Final 20 : 97270190955307463227 The digits of the 2^32 Fibonacci number (with string length 897595080) are: First 20 : 61999319689381859818 Final 20 : 39623735538208076347
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)
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
Nim
Using the Lucas method.
import strformat, strutils, times
import bignum
let
One = newInt(1)
Two = newInt(2)
Three = newInt(3)
proc lucas(n: Int): Int =
proc inner(n: Int): (Int, Int) =
if n.isZero: return (newInt(0), newInt(1))
var t = n shr 1
var (u, v) = inner(t)
t = n and Two
let q = t - One
var r = newInt(0)
u *= u
v *= v
t = n and One
if t == One:
t = (u - q) * Two
r = v * Three
result = (u + v, r - t)
else:
t = u * Three
r = v + q
r *= Two
result = (r - t, u + v)
var t = n shr 1
let (u, v) = inner(t)
let l = v * Two - u
t = n and One
if t == One:
let q = n and Two - One
return v * l + q
return u * l
let start = now()
var n: Int
var i = 10
while i <= 10_000_000:
n = newInt(i)
let s = $lucas(n)
echo &"The digits of the {($i).insertSep}th Fibonacci number ({($s.len).insertSep}) are:"
if s.len > 20:
echo &" First 20 : {s[0..19]}"
if s.len < 40:
echo &" Final {s.len-20:<2} : {s[20..^1]}"
else:
echo &" Final 20 : {s[^20..^1]}"
else:
echo &" All {s.len:<2} : {s}"
echo()
i *= 10
for e in [culong 16, 32]:
n = One shl e
let s = $lucas(n)
echo &"The digits of the 2^{e}th Fibonacci number ({($s.len).insertSep}) are:"
echo &" First 20 : {s[0..19]}"
echo &" Final 20 : {s[^20..^1]}"
echo()
echo &"Took {now() - start}"
- Output:
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^16th 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 6 minutes, 18 seconds, 643 milliseconds, 622 microseconds, and 965 nanoseconds
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;
}
- 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
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.)
without javascript_semantics -- (no mpfr_log(), mpfr_exp() under pwa/p2js) requires("1.0.0") -- (mpfr_set_default_prec[ision] has been renamed) include mpfr.e mpfr_set_default_precision(-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)
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.
with javascript_semantics 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 iff(platform()=JS?6: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]
Note that under pwa/p2js 10^6 takes 11.4s so we'll stop there...
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.
Python
With fancy custom integer classes that are absolutely useless except for this task.
class Head():
def __init__(self, lo, hi=None, shift=0):
if hi is None: hi = lo
d = hi - lo
ds, ls, hs = str(d), str(lo), str(hi)
if d and len(ls) > len(ds):
assert(len(ls) - len(ds) + 1 > 21)
lo = int(str(lo)[:len(ls) - len(ds) + 1])
hi = int(str(hi)[:len(hs) - len(ds) + 1]) + 1
shift += len(ds) - 1
elif len(ls) > 100:
lo = int(str(ls)[:100])
hi = lo + 1
shift = len(ls) - 100
self.lo, self.hi, self.shift = lo, hi, shift
def __mul__(self, other):
lo = self.lo*other.lo
hi = self.hi*other.hi
shift = self.shift + other.shift
return Head(lo, hi, shift)
def __add__(self, other):
if self.shift < other.shift:
return other + self
sh = self.shift - other.shift
if sh >= len(str(other.hi)):
return Head(self.lo, self.hi, self.shift)
ls = str(other.lo)
hs = str(other.hi)
lo = self.lo + int(ls[:len(ls)-sh])
hi = self.hi + int(hs[:len(hs)-sh])
return Head(lo, hi, self.shift)
def __repr__(self):
return str(self.hi)[:20]
class Tail():
def __init__(self, v):
self.v = int(f'{v:020d}'[-20:])
def __add__(self, other):
return Tail(self.v + other.v)
def __mul__(self, other):
return Tail(self.v*other.v)
def __repr__(self):
return f'{self.v:020d}'[-20:]
def mul(a, b):
return a[0]*b[0] + a[1]*b[1], a[0]*b[1] + a[1]*b[2], a[1]*b[1] + a[2]*b[2]
def fibo(n, cls):
n -= 1
zero, one = cls(0), cls(1)
m = (one, one, zero)
e = (one, zero, one)
while n:
if n&1: e = mul(m, e)
m = mul(m, m)
n >>= 1
return f'{e[0]}'
for i in range(2, 10):
n = 10**i
print(f'10^{i} :', fibo(n, Head), '...', fibo(n, Tail))
for i in range(3, 8):
n = 2**i
s = f'2^{n}'
print(f'{s:5s}:', fibo(2**n, Head), '...', fibo(2**n, Tail))
- Output:
10^2 : 35422484817926191507 ... 54224848179261915075 10^3 : 43466557686937456435 ... 76137795166849228875 10^4 : 33644764876431783266 ... 66073310059947366875 10^5 : 25974069347221724166 ... 49895374653428746875 10^6 : 19532821287077577316 ... 68996526838242546875 10^7 : 11298343782253997603 ... 86998673686380546875 10^8 : 47371034734563369625 ... 06082642167760546875 10^9 : 79523178745546834678 ... 03172326981560546875 2^8 : 14169381771405651323 ... 19657707794958199867 2^16 : 73199214460290552832 ... 97270190955307463227 2^32 : 61999319689381859818 ... 39623735538208076347 2^64 : 11175807536929528424 ... 17529800348089840187 2^128: 15262728879740471565 ... 17229324095882654267
Ruby
Matrix exponentiation by Ruby's fast exponentiation operator duck-typing applied to Ruby's built-in Integer Class
require 'matrix' start_time = Time.now [0,1,2,3,4,10,100,256, 1_000, 1024, 10_000, 2 ** 16, 100_000, 1_000_000,10_000_000 ].each {|n| ## fib_Num=(Matrix[[0,1],[1,1]] ** (n))[0,1] ## Matrix exponentiation ## fib_Str= fib_Num.to_s() if fib_Str.length <= 21 p ["Fibonacci(#{n})",fib_Str.length.to_s + ' digits' , fib_Str] else p ["Fibonacci(#{n})",fib_Str.length.to_s + ' digits' , fib_Str.slice(0,20) + " ... " + fib_Str.slice(-20,20)] end } puts "Took #{Time.now - start_time}s"
- Output:
["Fibonacci(0)", "1 digits", "0"] ["Fibonacci(1)", "1 digits", "1"] ["Fibonacci(2)", "1 digits", "1"] ["Fibonacci(3)", "1 digits", "2"] ["Fibonacci(4)", "1 digits", "3"] ["Fibonacci(10)", "2 digits", "55"] ["Fibonacci(100)", "21 digits", "354224848179261915075"] ["Fibonacci(256)", "54 digits", "14169381771405651323 ... 19657707794958199867"] ["Fibonacci(1000)", "209 digits", "43466557686937456435 ... 76137795166849228875"] ["Fibonacci(1024)", "214 digits", "45066996336778198131 ... 04103631553925405243"] ["Fibonacci(10000)", "2090 digits", "33644764876431783266 ... 66073310059947366875"] ["Fibonacci(65536)", "13696 digits", "73199214460290552832 ... 97270190955307463227"] ["Fibonacci(100000)", "20899 digits", "25974069347221724166 ... 49895374653428746875"] ["Fibonacci(1000000)", "208988 digits", "19532821287077577316 ... 68996526838242546875"] ["Fibonacci(10000000)", "2089877 digits", "11298343782253997603 ... 86998673686380546875"] Took 1.027948065s
Matrix exponentiation by Ruby's Matlix#exponentiation duck-typeing with Head-tail-BigNumber class
Here, the HeadTailBignum class holds the digits of the head part in Ruby's bigDecimal class and tail part in the BigInteger class. and HeadTailBignum defines only the add and multply operators and the coerce method.
Matlix exponentiation operator is fast and can apply duck-typing to HeadTailBignum.
require 'matrix' require 'bigdecimal' class HeadTailBignum < Numeric attr_accessor :hi, :lo @@degitsLimit = 20 @@hiDegitsLimit = (@@degitsLimit + 1) * 2 @@loDegitsBase = 10 ** @@degitsLimit BigDecimal::limit(@@hiDegitsLimit) def initialize(other, loValue = nil) if other.kind_of?(BigDecimal) && loValue.kind_of?(Integer) @hi = other @lo = loValue % @@loDegitsBase elsif other.kind_of?(HeadTailBignum) @hi = other.hi @lo = other.lo elsif other.kind_of?(Integer) @hi = BigDecimal(other) @lo = other % @@loDegitsBase else raise StandardError.new("HeadTailBignum initialize Type Error (" + other.class.to_s + "," + loValue.class.to_s + ")") end end def clone HeadTailBignum(self.hi, self.lo) end def integer?() true end def coerce(other) if other.kind_of?(Integer) return HeadTailBignum.new(other), self else super end end def +(other) if other.kind_of?(Integer) return HeadTailBignum.new(self.hi + other, self.lo + other) elsif other.kind_of?(HeadTailBignum) return HeadTailBignum.new(self.hi + other.hi , self.lo + other.lo) else raise StandardError.new("HeadTailBignum add Type Error (" + other.class.to_s + ")") end end def *(other) if other.kind_of?(Integer) return HeadTailBignum.new(self.hi * other, self.lo * other) elsif other.kind_of?(HeadTailBignum) return HeadTailBignum.new(self.hi * other.hi , self.lo * other.lo) else raise StandardError.new("HeadTailBignum mulply Type Error (" + other.class.to_s + ")") end end def exponent @hi.exponent end def to_s if @hi < @@loDegitsBase return @lo.inspect else return @hi.to_s("E").slice(2,@@degitsLimit ) + " ... " + @lo.to_s end end end start_time = Time.now [8,16,32,64].each {|h| n = (2**h) fib_Num=(Matrix[[ HeadTailBignum.new(0),1],[1,1]] ** (n))[0,1] puts "Fibonacci(2^#{h.to_s}) = #{fib_Num.to_s} are #{fib_Num.exponent} digits" } puts "Took #{Time.now - start_time}s"
- Output:
Fibonacci(2^8) = 14169381771405651323 ... 19657707794958199867 are 54 digits Fibonacci(2^16) = 73199214460290552832 ... 97270190955307463227 are 13696 digits Fibonacci(2^32) = 61999319689381859818 ... 39623735538208076347 are 897595080 digits Fibonacci(2^64) = 11175807536929528424 ... 17529800348089840187 are 3855141514259838964 digits Took 0.009357194s
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 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
}
- 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
Wren
Matrix exponentiation
A tough task for Wren which takes just under 3 minutes to process even the 1 millionth Fibonacci number. Judging by the times for the compiled, statically typed languages using GMP, the 10 millionth would likely take north of 4 hours so I haven't attempted it.
import "./big" for BigInt
import "./fmt" for Fmt
var mul = Fn.new { |m1, m2|
var rows1 = m1.count
var cols1 = m1[0].count
var rows2 = m2.count
var cols2 = m2[0].count
if (cols1 != rows2) Fiber.abort("Matrices cannot be multiplied.")
var result = List.filled(rows1, null)
for (i in 0...rows1) {
result[i] = List.filled(cols2, null)
for (j in 0...cols2) {
result[i][j] = BigInt.zero
for (k in 0...rows2) {
var temp = m1[i][k] * m2[k][j]
result[i][j] = result[i][j] + temp
}
}
}
return result
}
var identityMatrix = Fn.new { |n|
if (n < 1) Fiber.abort("Size of identity matrix can't be less than 1")
var ident = List.filled(n, 0)
for (i in 0...n) {
ident[i] = List.filled(n, null)
for(j in 0...n) ident[i][j] = (i != j) ? BigInt.zero : BigInt.one
}
return ident
}
var pow = Fn.new { |m, n|
var le = m.count
if (le != m[0].count) Fiber.abort("Not a square matrix")
if (n < 0) Fiber.abort("Negative exponents not supported")
if (n == 0) return identityMatrix.call(le)
if (n == 1) return m
var pow = identityMatrix.call(le)
var base = m
var e = n
while (e > 0) {
var temp = e & BigInt.one
if (temp == BigInt.one) pow = mul.call(pow, base)
e = e >> 1
base = mul.call(base, base)
}
return pow
}
var fibonacci = Fn.new { |n|
if (n == 0) return BigInt.zero
var m = [[BigInt.one, BigInt.one], [BigInt.one, BigInt.zero]]
m = pow.call(m, n - 1)
return m[0][0]
}
var n = BigInt.zero
var i = 10
while (i <= 1e6) {
n = BigInt.new(i)
var s = fibonacci.call(n).toString
Fmt.print("The digits of the $,sth Fibonacci number ($,s) are:", i, s.count)
if (s.count > 20) {
Fmt.print(" First 20 : $s", s[0...20])
if (s.count < 40) {
Fmt.print(" Final $-2d : $s", s.count-20, s[20..-1])
} else {
Fmt.print(" Final 20 : $s", s[s.count-20..-1])
}
} else {
Fmt.print(" All $-2d : $s", s.count, s)
}
System.print()
i = i * 10
}
n = BigInt.one << 16
var s = fibonacci.call(n).toString
Fmt.print("The digits of the 2^16th Fibonacci number ($,s) are:", s.count)
Fmt.print(" First 20 : $s", s[0...20])
Fmt.print(" Final 20 : $s", s[s.count-20..-1])
- Output:
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 2^16th Fibonacci number (13,696) are: First 20 : 73199214460290552832 Final 20 : 97270190955307463227
Lucas method
Much quicker than the matrix exponentiation method taking about 23 seconds to process the 1 millionth Fibonacci number and around 32 minutes to reach the 10 millionth.
Apart from the 2^16th number, the extra credit is still out of reach using this approach.
import "./big" for BigInt
import "./fmt" for Fmt
var lucas = Fn.new { |n|
var inner // recursive function
inner = Fn.new { |n|
if (n == BigInt.zero) return [BigInt.zero, BigInt.one]
var t = n >> 1
var res = inner.call(t)
var u = res[0]
var v = res[1]
t = n & BigInt.two
var q = t - BigInt.one
var r = BigInt.zero
u = u.square
v = v.square
t = n & BigInt.one
if (t == BigInt.one) {
t = u - q
t = BigInt.two * t
r = BigInt.three * v
return [u + v, r - t]
} else {
t = BigInt.three * u
r = v + q
r = BigInt.two * r
return [r - t, u + v]
}
}
var t = n >> 1
var res = inner.call(t)
var u = res[0]
var v = res[1]
var l = BigInt.two * v
l = l - u // Lucas function
t = n & BigInt.one
if (t == BigInt.one) {
var q = n & BigInt.two
q = q - BigInt.one
t = v * l
return t + q
}
return u * l
}
var n = BigInt.zero
var i = 10
while (i <= 1e7) {
n = BigInt.new(i)
var s = lucas.call(n).toString
Fmt.print("The digits of the $,sth Fibonacci number ($,s) are:", i, s.count)
if (s.count > 20) {
Fmt.print(" First 20 : $s", s[0...20])
if (s.count < 40) {
Fmt.print(" Final $-2d : $s", s.count-20, s[20..-1])
} else {
Fmt.print(" Final 20 : $s", s[s.count-20..-1])
}
} else {
Fmt.print(" All $-2d : $s", s.count, s)
}
System.print()
i = i * 10
}
n = BigInt.one << 16
var s = lucas.call(n).toString
Fmt.print("The digits of the 2^16th Fibonacci number ($,s) are:", s.count)
Fmt.print(" First 20 : $s", s[0...20])
Fmt.print(" Final 20 : $s", s[s.count-20..-1])
- Output:
As matrix exponentiation method, plus:
The digits of the 10,000,000th Fibonacci number (2,089,877) are: First 20 : 11298343782253997603 Final 20 : 86998673686380546875
Fibmod method (embedded)
Ridiculously fast (under 5ms) thanks to the stuff borrowed from Sidef and Julia combined with the speed of GMP.
import "./gmp" for Mpz, Mpf
import "./fmt" for Fmt
var nd = 20 // number of digits to be displayed at each end
var pr = 128 // precision to be used
var one = Mpz.one
var two = Mpz.two
var fibmod = Fn.new { |n, nmod|
if (n < two) return n
var fibmods = {}
var f // recursive closure
f = Fn.new { |n|
if (n < two) return one
var ns = n.toString
var v = fibmods[ns]
if (v) return v
v = Mpz.zero
var k = n / two
var t = n & one
var u = Mpz.zero
if (t != one) {
t.set(f.call(k)).square
v.sub(k, one)
u.set(f.call(v)).square
} else {
t.set(f.call(k))
v.add(k, one)
v.set(f.call(v))
u.sub(k, one)
u.set(f.call(u)).mul(t)
t.mul(v)
}
t.add(u)
fibmods[ns] = t.rem(nmod)
return fibmods[ns]
}
var w = n - one
return f.call(w)
}
var binetApprox = Fn.new { |n|
var phi = Mpf.from(0.5, pr)
var ihp = phi.copy()
var root = Mpf.from(1.25, pr).sqrt
phi.add(root)
ihp.sub(root, ihp).neg
ihp.sub(phi, ihp).log
phi.log
var nn = Mpf.from(n, pr)
return phi.mul(nn).sub(ihp)
}
var firstFibDigits = Fn.new { |n, k|
var f = binetApprox.call(n)
var g = Mpf.new(pr)
g.div(f, Mpf.ln10(pr)).inc
g.floor.mul(Mpf.ln10(pr))
f.sub(g).exp
var p = Mpz.from(10).pow(k)
g.set(p)
f.mul(g)
return f.floor.toString[0...k]
}
var lastFibDigits = Fn.new { |n, k|
var p = Mpz.from(10).pow(k)
return fibmod.call(n, p).toString[0...k]
}
var start = System.clock
var n = Mpz.zero
var i = 10
while (i <= 1e7) {
n.set(i)
var nn = Mpz.from(i)
Fmt.print("\nThe digits of the $,r Fibonacci number are:", i)
var nd2 = nd
var nd3 = 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
}
var s1 = firstFibDigits.call(n, nd2)
if (s1.count < 20) {
Fmt.print(" All $-2d : $s", s1.count, s1)
} else {
Fmt.print(" First 20 : $s", s1)
var s2 = lastFibDigits.call(nn, nd3)
if (s2.count < 20) {
Fmt.print(" Final $-2d : $s", s2.count, s2)
} else {
Fmt.print(" Final 20 : $s", s2)
}
}
i = i * 10
}
var ord = ["th", "nd", "th"]
i = 0
for (p in [16, 32, 64]) {
n.lsh(one, p)
var nn = one << p
Fmt.print("\nThe digits of the 2^$d$s Fibonacci number are:", p, ord[i])
Fmt.print(" First $d : $s", nd, firstFibDigits.call(n, nd))
Fmt.print(" Final $d : $s", nd, lastFibDigits.call(nn, nd))
i = i + 1
}
System.print("\nTook %(System.clock-start) seconds.")
- 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 0.004516 seconds.