Fibonacci matrix-exponentiation: Difference between revisions
m
→{{header|Wren}}: Minor tidy
(→{{header|jq}}: Fib(2^n)) |
m (→{{header|Wren}}: Minor tidy) |
||
| (12 intermediate revisions by 4 users not shown) | |||
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,035:
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====
<
($m1|length) as $rows1
| ($m1[0]|length) as $cols1
Line 1,108:
" Final 20 : \($s[-20:])";
task1, "", task2</
{{out}}
<pre>
Line 1,140:
====Fib(2^n)====
<
# Output: Fib(2^n)
def Fib2pN:
Line 1,162:
" First 20 : \(.[0:20])",
" Final 20 : \(.[-20:])",
""</
{{out}}
Using gojq to compute Fib(2^32) using this method takes many hours.
Line 1,171:
The digits of the 2^32 Fibonacci number (with string length 897595080) are:
First 20 :
Final 20 :
</pre>
Line 1,181 ⟶ 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,244 ⟶ 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,257 ⟶ 1,255:
{{trans|Go}}
Using the Lucas method.
<
import bignum
Line 1,325 ⟶ 1,323:
echo()
echo &"Took {now() - start}"</
{{out}}
Line 1,367 ⟶ 1,365:
=={{header|Perl}}==
{{trans|Sidef}}
<
use warnings;
Line 1,406 ⟶ 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,418 ⟶ 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,507 ⟶ 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,517 ⟶ 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,589 ⟶ 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,607 ⟶ 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,716 ⟶ 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,725 ⟶ 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,741 ⟶ 1,970:
say "F(2^#{n}) = #{a} ... #{b}"
}
}</
{{out}}
<pre>
Line 1,749 ⟶ 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,768 ⟶ 1,997:
var last20 = nth_fib_last_k_digits(2**n, 20)
say "F(2^#{n}) = #{first20} ... #{last20}"
}</
{{out}}
<pre>
Line 1,782 ⟶ 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,863 ⟶ 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,900 ⟶ 2,129:
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,968 ⟶ 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,982 ⟶ 2,211:
{{libheader|Wren-gmp}}
Ridiculously fast (under 5ms) thanks to the stuff borrowed from Sidef and Julia combined with the speed of GMP.
<
import "./fmt" for Fmt
Line 2,093 ⟶ 2,322:
i = i + 1
}
System.print("\nTook %(System.clock-start) seconds.")</
{{out}}
| |||