Fibonacci matrix-exponentiation: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (→{{header|Wren}}: Typo) |
m (→{{header|Wren}}: Minor tidy) |
||
| (19 intermediate revisions by 8 users not shown) | |||
Line 7:
;Task:
Write a program using ''matrix exponentiation'' to generate Fibonacci(n) for '''n''' equal to:
:::* 10
:::* 100
Line 30:
=={{header|C sharp|C#}}==
===Matrix exponentiation===
<
using System.IO;
using System.Numerics;
Line 124:
}
</syntaxhighlight>
{{out}}
<pre>
Line 143:
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)}}
<
import (
Line 276:
fmt.Printf("Took %s\n\n", time.Since(start))
}</
{{out}}
Line 323:
<br>
{{trans|Julia}}
<
import (
Line 419:
fmt.Printf("Took %s\n\n", time.Since(start))
}</
{{out}}
Line 433:
{{trans|Sidef}}
{{trans|Julia}}
<
import (
Line 587:
}
fmt.Printf("\nTook %s\n", time.Since(start))
}</
{{out}}
Line 635:
=={{header|Haskell}}==
===Matrix exponentiation===
<
import Data.List
Line 702:
fromInteger n = Mat [[fromInteger n]]
abs = undefined
signum = undefined</
{{out}}
Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5:
Line 718:
</pre>
===Matrix exponentiation - printing alternative ===
<
import Data.List
Line 796:
fromInteger n = Mat [[fromInteger n]]
abs = undefined
signum = undefined</
{{out}}
Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5:
Line 820:
At each step of the exponentiation of a symmetric matric, we multiply 2 symmetric matrices which commute.
<
import System.CPUTime (getCPUTime)
import Data.List
Line 908:
power m 1 = m
power m n = if even n then w else mult w m
where w = square.power m.quot n $ 2</
{{out}}
Timing is on Intel Core i5-4300U CPU, Windows 10 Professional, using GHCi Version 8.6.5:
Line 927:
Implemented fib and fibMod.
<
import java.math.BigInteger;
import java.util.Arrays;
Line 1,010:
}
</syntaxhighlight>
{{out}}
Line 1,024:
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}}
Using gojq to compute Fib(2^32) using this method takes many hours.
<pre>
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
</pre>
Line 1,030 ⟶ 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
Lucas method is used as the alternative method.
<
const b = [big"1" 1; 1 0]
matrixfibonacci(n) = n == 0 ? 0 : n == 1 ? 1 : (b^(n+1))[2,2]
Line 1,093 ⟶ 1,242:
rpad(firstlast(firstbinet(32) * lastfibmod(32)), 45))
println("2^64 ", " "^90, rpad(firstlast(firstbinet(64) * lastfibmod(64)), 45))
</
<pre>
N Matrix Lucas Mod
Line 1,101 ⟶ 1,250:
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}}
<
use warnings;
Line 1,143 ⟶ 1,404:
my $last20 = nth_fib_last_k_digits(2**$n, 20);
printf "F(2^$n) = %s ... %s\n", $first20, $last20;
}</
{{out}}
<pre>F(2^16) = 73199214460290552832 ... 97270190955307463227
Line 1,155 ⟶ 1,416:
{{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.)
<!--<
<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>
Line 1,244 ⟶ 1,505:
<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>
<!--</
{{out}}
<pre>
Line 1,254 ⟶ 1,515:
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.
<!--<
<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>
Line 1,326 ⟶ 1,587:
<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>
<!--</
{{out}}
<pre>
Line 1,344 ⟶ 1,605:
</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
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"
{{out}}
<pre>
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
</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"
use Inline::Perl5;
my $p5 = Inline::Perl5.new();
Line 1,453 ⟶ 1,945:
my $last20 = nth_fib_last_k_digits(2**$n, 20);
printf "F(2^$n) = %s ... %s\n", $first20, $last20
}</
{{out}}
<pre>F(2^16) = 73199214460290552832 ... 97270190955307463227
Line 1,462 ⟶ 1,954:
=={{header|Sidef}}==
Computing the n-th Fibonacci number, using matrix-exponentiation (this function is also built-in):
<
([[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:
<
var f = fibonacci(2**n)
Line 1,478 ⟶ 1,970:
say "F(2^#{n}) = #{a} ... #{b}"
}
}</
{{out}}
<pre>
Line 1,486 ⟶ 1,978:
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:
<
const PHI = (1.25.sqrt + 0.5)
const IHP = -(1.25.sqrt - 0.5)
Line 1,505 ⟶ 1,997:
var last20 = nth_fib_last_k_digits(2**n, 20)
say "F(2^#{n}) = #{first20} ... #{last20}"
}</
{{out}}
<pre>
Line 1,519 ⟶ 2,011:
{{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.
<
import "./fmt" for Fmt
var mul = Fn.new { |m1, m2|
Line 1,600 ⟶ 2,092:
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])</
{{out}}
Line 1,636 ⟶ 2,128:
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 "./fmt" for Fmt
var lucas = Fn.new { |n|
Line 1,705 ⟶ 2,197:
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])</
{{out}}
Line 1,713 ⟶ 2,205:
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>
| |||