Exploring Fibonacci numbers with math/big
1,095 views
Skip to first unread message
Michael Jones
Mar 26, 2012, 6:01:59 PM3/26/12
to golang-nuts
Following the recent discussions about using math/big I decided to write a quick example based around Fibonacci numbers. The idea here is a fast logarithmic time approach that recurses to compute F(n) from F(n/2). The code is optimized to find an individual F(n) rather than the series from F(0)..F(n), however, it is very fast even for F(1000000) and F(10000000) so I used it to scan through the first few hundred numbers with math/big's probable primality code to replicate sequence A001605 in OEIS, the first prime Fibonacci numbers. (Which also reminds that all known prime Fibonacci numbers have a prime index, except F(4).)
Go Playground: http://play.golang.org/p/h1PrDY_1u-
Michael
P.S. (For posterity ;-) The reason the else-clause is not coded normally at the end of fib_half() is because of present Go compiler limitations. When this limit is removed the code should change.)
CODE:
package main
import (
"fmt"
"math/big"
)
// recursive splitting Fibonacci evaluation
func fib_half(n int) (*big.Int, *big.Int) {
if n <= 0 {
zero := big.NewInt(0)
return zero, zero
}
if n == 1 {
return big.NewInt(0), big.NewInt(1)
}
h := n / 2 // compute F(n) from F(n/2) for O(log2 n) complexity
f1, f2 := fib_half(h) // f1 = F(h-1), f2 = F(h)
f3 := big.NewInt(0).Add(f1, f2) // f3 = F(h+1)
t := big.NewInt(0).Add(f1, f3)
t = t.Mul(t, f2)
s2 := big.NewInt(0).Mul(f2, f2) // F(h)**2
if n&1 == 1 { // handle odd n
s3 := big.NewInt(0).Mul(f3, f3) // F(h+1)**2
return t, s2.Add(s2, s3) // (f1 + f3)*f2, f2*f2 + f3*f3
} // else // handle even n
s1 := big.NewInt(0).Mul(f1, f1) // F(h-1)**2
return s1.Add(s1, s2), t // f1*f1 + f2*f2, (f1 + f3)*f2
}
func fib(n int) *big.Int {
_, f := fib_half(n)
return f
}
func trim(n *big.Int) string {
s := fmt.Sprintf("%d", n)
width := 91
length := len(s)
if length > width {
s = s[0:(width-3)/2] + "..." + s[length-(width-3)/2:]
}
s = fmt.Sprintf("%*s [%7d digits]", width, s, length)
return s
}
var fibTests = []struct {
n int
f string
}{
{0, "0"},
{1, "1"},
{2, "1"},
{3, "2"},
{4, "3"},
{480, "9216845717656874712980450562726202415567360565980794777111390850331644813674856981646960226192287360"},
{485, "102216385354134324778078911974689347564945335253989394162085272183728832930964492342791374522266860485"},
{490, "1133597084613134447271848482284309025629966048359864130560049384871348807054284272752352079971127752695"},
{495, "12571784316098613244768412217102088629494571867212494830322628505768565710528091492618664254204672140130"},
{500, "139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125"},
}
func testFib() int {
failures := 0
for i, test := range fibTests {
s := fmt.Sprintf("%d", fib(test.n))
if s != test.f {
fmt.Printf("error: #%d got %v want %v\n", i, s, test.f)
failures++
}
}
if failures == 0 {
fmt.Printf("passed\n")
} else {
fmt.Printf("failed: %d failures\n", failures)
}
return failures
}
func main() {
fmt.Printf("testing... ")
testFib()
fmt.Printf("\n")
fmt.Printf("first\n")
for n := 0; n <= 16; n++ {
fmt.Printf(" F(%2d) = %d\n", n, fib(n))
}
fmt.Printf("\n")
fmt.Printf("F(2**n)\n")
for n := 1; n <= 1<<20; n <<= 1 {
fmt.Printf(" F(%7d) = %s\n", n, trim(fib(n)))
}
fmt.Printf("\n")
fmt.Printf("F(10**n)\n")
for n := 1; n <= 1000000; n *= 10 {
fmt.Printf(" F(%7d) = %s\n", n, trim(fib(n)))
}
fmt.Printf("\n")
fmt.Printf("probable prime Fibonacci numbers (sequence A001605 in OEIS)\n")
for n := 0; n <= 600; n++ {
f := fib(n)
if f.ProbablyPrime(64) {
fmt.Printf(" F(%7d) = %s\n", n, trim(f))
}
}
fmt.Printf("\n")
}
OUTPUT:
mtj$ go run fib.go
testing... passed
first
F( 0) = 0
F( 1) = 1
F( 2) = 1
F( 3) = 2
F( 4) = 3
F( 5) = 5
F( 6) = 8
F( 7) = 13
F( 8) = 21
F( 9) = 34
F(10) = 55
F(11) = 89
F(12) = 144
F(13) = 233
F(14) = 377
F(15) = 610
F(16) = 987
F(2**n)
F( 1) = 1 [ 1 digits]
F( 2) = 1 [ 1 digits]
F( 4) = 3 [ 1 digits]
F( 8) = 21 [ 2 digits]
F( 16) = 987 [ 3 digits]
F( 32) = 2178309 [ 7 digits]
F( 64) = 10610209857723 [ 14 digits]
F( 128) = 251728825683549488150424261 [ 27 digits]
F( 256) = 141693817714056513234709965875411919657707794958199867 [ 54 digits]
F( 512) = 44893845313309942978077298160660626646181883...66661322268805744081870933775586567858979269 [ 107 digits]
F( 1024) = 45066996336778198131043832357288860493678605...68899711089246778413354004103631553925405243 [ 214 digits]
F( 2048) = 45415304437437894250455714462906892027009082...10145547104337536614028338983073680262684101 [ 428 digits]
F( 4096) = 46120017322804312474564457085636141271732249...85050077521931600327064840520442470066917947 [ 856 digits]
F( 8192) = 47562418031541483249676431926694568862149480...34121942551466950077266295376251851089202629 [ 1712 digits]
F( 16384) = 50583963273256865512698275499738614264772589...16487934517684664695772372793939867475669563 [ 3424 digits]
F( 32768) = 57215106297689468939590948384261448582742164...43607205126221971914921497265149891505753541 [ 6848 digits]
F( 65536) = 73199214460290552832639915812598454538952425...17531793911596579203414697270190955307463227 [ 13696 digits]
F( 131072) = 11981131726582568764552857855321239080658488...44491510870363731782250304265906436130344389 [ 27393 digits]
F( 262144) = 32098200701891878424611639413278177386990176...46534756609059187480173020581102949038682683 [ 54785 digits]
F( 524288) = 23038085126797742368512278649433257358762739...35892895641064341683333020246122520212016581 [ 109570 digits]
F(1048576) = 11868006063550661146017510449861600696082418...24222218017383101627183453200033680691163707 [ 219140 digits]
F(10**n)
F( 1) = 1 [ 1 digits]
F( 10) = 55 [ 2 digits]
F( 100) = 354224848179261915075 [ 21 digits]
F( 1000) = 43466557686937456435688527675040625802564660...38298969649928516003704476137795166849228875 [ 209 digits]
F( 10000) = 33644764876431783266621612005107543310302148...62430701794976171121233066073310059947366875 [ 2090 digits]
F( 100000) = 25974069347221724166155034021275915414880485...09970015699780289236362349895374653428746875 [ 20899 digits]
F(1000000) = 19532821287077577316320149475962563324435429...19033684684301719893411568996526838242546875 [ 208988 digits]
probable prime Fibonacci numbers (sequence A001605 in OEIS)
F( 3) = 2 [ 1 digits]
F( 4) = 3 [ 1 digits]
F( 5) = 5 [ 1 digits]
F( 7) = 13 [ 2 digits]
F( 11) = 89 [ 2 digits]
F( 13) = 233 [ 3 digits]
F( 17) = 1597 [ 4 digits]
F( 23) = 28657 [ 5 digits]
F( 29) = 514229 [ 6 digits]
F( 43) = 433494437 [ 9 digits]
F( 47) = 2971215073 [ 10 digits]
F( 83) = 99194853094755497 [ 17 digits]
F( 131) = 1066340417491710595814572169 [ 28 digits]
F( 137) = 19134702400093278081449423917 [ 29 digits]
F( 359) = 475420437734698220747368027166749382927701417016557193662268716376935476241 [ 75 digits]
F( 431) = 529892711006095621792039556787784670197112759029534506620905162834769955134424689676262369 [ 90 digits]
F( 433) = 1387277127804783827114186103186246392258450358171783690079918032136025225954602593712568353 [ 91 digits]
F( 449) = 30617199924845450305543138480837172081112854...38497131674799321571238149015933442805665949 [ 94 digits]
F( 509) = 10597999265301490732599643671505003412515860...80142974347195483140293254396195769876129909 [ 107 digits]
F( 569) = 36684474316080978061473613646275630451100586...93560857904335017879540515228143777781065869 [ 119 digits]
F( 571) = 96041200618922553823942883360924865026104917...08478864192589084185254331637646183008074629 [ 119 digits]
Steven Blenkinsop
Mar 26, 2012, 10:56:02 PM3/26/12
to golang-nuts
I loopified the code <http://play.golang.org/p/SGS6a2LTb6>, and discovered the obvious: any gains by consolidating big.Int instances and constant stack space are imperceptible next to the allocation and computation involved in big arithmetic. Who'da thunk it? ;)
Steven Blenkinsop
Mar 26, 2012, 11:16:10 PM3/26/12
to golang-nuts
On Mon, Mar 26, 2012 at 10:56 PM, Steven Blenkinsop <stev...@gmail.com> wrote:
I loopified the code <http://play.golang.org/p/SGS6a2LTb6>, and discovered the obvious: any gains by consolidating big.Int instances and constant stack space are imperceptible next to the allocation and computation involved in big arithmetic. Who'da thunk it? ;)
Scratch that: the conversion to string is by far the limiting factor, it would seem. Though I'm sure it's a lot faster than it used to be :).
Michael Jones
Mar 27, 2012, 9:08:30 AM3/27/12
to Steven Blenkinsop, golang-nuts
Nice! Yes, string conversion is expensive, even though it's > 300x faster than before. Hexadecimal, octal, or binary output formatting are much faster ("%x" > "%o" > "%b".) Memory allocation is very expensive too. Best is when operators are provided with temporaries...
a.op(a, b) means { new t; t = a op b; a = t; free t }
t.op(a, b) means { t = a op b }
--
Steven Blenkinsop
Mar 27, 2012, 12:10:57 PM3/27/12
to Michael Jones, golang-nuts
On Tue, Mar 27, 2012 at 9:08 AM, Michael Jones <m...@google.com> wrote:
Nice! Yes, string conversion is expensive, even though it's > 300x faster than before.
Wow! I new it had gotten faster, but I didn't realize by how much :)
Hexadecimal, octal, or binary output formatting are much faster ("%x" > "%o" > "%b".) Memory allocation is very expensive too. Best is when operators are provided with temporaries...a.op(a, b) means { new t; t = a op b; a = t; free t }t.op(a, b) means { t = a op b }
Thanks, it didn't occur to me that big was making temporaries for me. So this would be a better way to do it?
Michael Jones
Mar 27, 2012, 1:02:01 PM3/27/12
to Steven Blenkinsop, golang-nuts
Benchmark. I'm on an airplane and can't do it at the moment.
Michael Jones
Mar 27, 2012, 1:09:07 PM3/27/12
to Steven Blenkinsop, golang-nuts
Certainly looks much better!
Of course, in your programs and mine, if we really wanted a dense range of solutions then the simple summation would be the best answer after the first two. ;-)
There is a minimal work algorithm involving Lucas multipliers that I'll code in go for comparison once I take off.
Steven Blenkinsop
Mar 27, 2012, 8:53:45 PM3/27/12
to Michael Jones, golang-nuts
On Tue, Mar 27, 2012 at 1:02 PM, Michael Jones <m...@google.com> wrote:
Benchmark. I'm on an airplane and can't do it at the moment
Not sure the best benchmark strategy. I just computed the same fib(N) b.N times for each function, and got these results. Fib1 is your original one. Fib2 is my original one. Fib3 is my modified one.
N = 100
BenchmarkFib1 200000 7281 ns/op
BenchmarkFib2 1000000 2049 ns/op
BenchmarkFib3 1000000 2448 ns/op
N = 1,000
BenchmarkFib1 200000 11945 ns/op
BenchmarkFib2 500000 4711 ns/op
BenchmarkFib3 500000 5192 ns/op
N = 10,000
BenchmarkFib1 50000 47952 ns/op
BenchmarkFib2 50000 37102 ns/op
BenchmarkFib3 50000 39687 ns/op
N = 100,000
BenchmarkFib1 1000 1458540 ns/op
BenchmarkFib2 1000 1435098 ns/op
BenchmarkFib3 1000 1442553 ns/op
N = 1,000,000
BenchmarkFib1 10 132721800 ns/op
BenchmarkFib2 10 133016900 ns/op
BenchmarkFib3 10 133108000 ns/op
Fib2 and Fib3 are faster than Fib1, though the difference becomes less significant as N increases. However, it appears Fib3 is actually slower than Fib2, so the changes I made had opposite the desired effect.
Michael Jones
Mar 28, 2012, 12:54:25 AM3/28/12
to Steven Blenkinsop, golang-nuts
Very nice! I show something similar here.
Old way (my code, EWD way, recursive):
testing...
passed
timing...
1.786822 seconds to compute F(0)..F(100) 2000 times
2.075858 seconds to compute F(1000)..F(1100) 1000 times
1.599733 seconds to compute F(10000)..F(10100) 200 times
2.515566 seconds to compute F(100000)..F(100100) 10 times
2.679213 seconds to compute F(1000000)..F(1000000) 10 times
New way (your code, iterative using simple/elegant formulation and
minimal temporaries):
testing...
passed
timing...
0.502829 seconds to compute F(0)..F(100) 2000 times
0.814404 seconds to compute F(1000)..F(1100) 1000 times
1.297491 seconds to compute F(10000)..F(10100) 200 times
2.444358 seconds to compute F(100000)..F(100100) 10 times
2.719326 seconds to compute F(1000000)..F(1000000) 10 times
New way modified to reuse temporaries:
testing...
passed
timing...
0.633315 seconds to compute F(0)..F(100) 2000 times
0.889963 seconds to compute F(1000)..F(1100) 1000 times
1.358165 seconds to compute F(10000)..F(10100) 200 times
2.498833 seconds to compute F(100000)..F(100100) 10 times
2.82191 seconds to compute F(1000000)..F(1000000) 10 times
Your original iterative version is the fastest so far. Bravo!
Michael
--
Reply all
Reply to author
Forward
0 new messages