The Billionth Fibonacci Number
A few weeks ago, rockstarling wrote a post about calculating large Fibonacci numbers. Never one to resist a program challenge, I really wanted to try my hand at it and see if I could write an even faster program for calculating Fibonacci numbers. At the time, I was rather busy working on a paper but that's done now, so I gave it a shot today. Here's what I came up with:
#include
#include
#include
int main(
int argc,
char **argv )
{
int nth = atoi( argv[ 1 ] );
mpz_t a, b, temp1, temp2;
mpz_init2( a, 800000000 );
mpz_init2( b, 800000000 );
mpz_init2( temp1, 800000000 );
mpz_init2( temp2, 800000000 );
mpz_set_ui( a, 0 );
mpz_set_ui( b, 1 );
int bit = 1;
for ( ; bit <= nth; bit <<= 1 );
for ( bit >>= 1; bit; bit >>= 1 )
{
mpz_mul( temp1, a, a );
mpz_mul( temp2, a, b );
mpz_add( temp2, temp2, temp2 );
mpz_add( a, temp1, temp2 );
mpz_mul( temp2, b, b );
mpz_add( b, temp1, temp2 );
if ( nth & bit )
{
mpz_add( a, a, b );
mpz_sub( b, a, b );
}
}
mpz_out_str( stdout, 10, a );
mpz_clear( a );
mpz_clear( b );
mpz_clear( temp1 );
mpz_clear( temp2 );
return 0;
}
Well, okay, that's not really C++ -- it should also compile as C99. Most of the real work is done by the GNU multiple precision arithmetic library, anyway. (Yes, I know they already have mpz_fib_ui() -- that's cheating!) This also uses a different algorithm than rockstarling's though they're both logarithmic in the number of basic steps. Whereas the algorithm that she implemented is usually attributed to Dijkstra, my program implements the algorithm outlined by Gosper and Salamin. The loop essentially implements repeated squaring, but implemented slightly unusually with bit-logic to iterate from the most to least significant bits instead of doing the standard recursive odd/even thing. In the end, it's the same result but this approach lets it do all the arithmetic in-place. The other nice thing about the Gosper and Salamin algorithm is that it only requires three multiplications per step, instead of the four that Dijkstra's seems to need.
Giving this all a whirl on a big 3GHz Opteron machine:
[hex:~/Temp] aek% time ./fib 100000000 | head -c 75
473710347345633696254897131335103865754868289377201030193412163431081491642
real 0m29.385s
user 0m29.238s
sys 0m0.156s
Not bad. It all checks out with the results that rockstarling reported. For some real fun, though, I've tried calculating the billionth Fibonacci number!
[hex:~/Temp] aek% time ./fib 1000000000 > billionth.txt
real 9m40.174s
user 9m34.304s
sys 0m5.888s
This number has 208,987,640 decimal digits with the first and last hundred being:
7952317874554683467829385196197148189255542185234398913453039937343246686182519370050999626136556779...
6320554356934532581033548032811298867411661975017768296198982902073319952559425703172326981560546875
It is obviously composite, as it passes the divisible-by-5 test.
The sequence "12345678" appears once, beginning at the 159,774,853 digit. The number "1248163264" appears also, beginning at the 169,057,940 digit. And so does "31415927", beginning with the 140,853,142 digit. The longest contiguous sequences of any single digit are "555555555" and "777777777" at nine digits each. The sequence "1337" occurs 20,787 times.
All digits are fairly equally represented, but zero is the most common digit, appearing 20,906,722 times. One is the least common digit, appearing only 20,889,487 times. By descending order of frequency, the digits are 0,5,2,3,8,4,7,6,9,1.
One other interesting tidbit: it's actually the conversion of the final result to decimal for output that takes the most time. If I comment out the call to mpz_out_str(), it finishes in 69 seconds.
This file is obviously a little large to post somewhere for download but if anybody would like to try confirming my results, the MD5 sum on the file is 10f1f75f783b7ae9802f38a6921ac561.
#include
#include
#include
int main(
int argc,
char **argv )
{
int nth = atoi( argv[ 1 ] );
mpz_t a, b, temp1, temp2;
mpz_init2( a, 800000000 );
mpz_init2( b, 800000000 );
mpz_init2( temp1, 800000000 );
mpz_init2( temp2, 800000000 );
mpz_set_ui( a, 0 );
mpz_set_ui( b, 1 );
int bit = 1;
for ( ; bit <= nth; bit <<= 1 );
for ( bit >>= 1; bit; bit >>= 1 )
{
mpz_mul( temp1, a, a );
mpz_mul( temp2, a, b );
mpz_add( temp2, temp2, temp2 );
mpz_add( a, temp1, temp2 );
mpz_mul( temp2, b, b );
mpz_add( b, temp1, temp2 );
if ( nth & bit )
{
mpz_add( a, a, b );
mpz_sub( b, a, b );
}
}
mpz_out_str( stdout, 10, a );
mpz_clear( a );
mpz_clear( b );
mpz_clear( temp1 );
mpz_clear( temp2 );
return 0;
}
Well, okay, that's not really C++ -- it should also compile as C99. Most of the real work is done by the GNU multiple precision arithmetic library, anyway. (Yes, I know they already have mpz_fib_ui() -- that's cheating!) This also uses a different algorithm than rockstarling's though they're both logarithmic in the number of basic steps. Whereas the algorithm that she implemented is usually attributed to Dijkstra, my program implements the algorithm outlined by Gosper and Salamin. The loop essentially implements repeated squaring, but implemented slightly unusually with bit-logic to iterate from the most to least significant bits instead of doing the standard recursive odd/even thing. In the end, it's the same result but this approach lets it do all the arithmetic in-place. The other nice thing about the Gosper and Salamin algorithm is that it only requires three multiplications per step, instead of the four that Dijkstra's seems to need.
Giving this all a whirl on a big 3GHz Opteron machine:
[hex:~/Temp] aek% time ./fib 100000000 | head -c 75
473710347345633696254897131335103865754868289377201030193412163431081491642
real 0m29.385s
user 0m29.238s
sys 0m0.156s
Not bad. It all checks out with the results that rockstarling reported. For some real fun, though, I've tried calculating the billionth Fibonacci number!
[hex:~/Temp] aek% time ./fib 1000000000 > billionth.txt
real 9m40.174s
user 9m34.304s
sys 0m5.888s
This number has 208,987,640 decimal digits with the first and last hundred being:
7952317874554683467829385196197148189255542185234398913453039937343246686182519370050999626136556779...
6320554356934532581033548032811298867411661975017768296198982902073319952559425703172326981560546875
It is obviously composite, as it passes the divisible-by-5 test.
The sequence "12345678" appears once, beginning at the 159,774,853 digit. The number "1248163264" appears also, beginning at the 169,057,940 digit. And so does "31415927", beginning with the 140,853,142 digit. The longest contiguous sequences of any single digit are "555555555" and "777777777" at nine digits each. The sequence "1337" occurs 20,787 times.
All digits are fairly equally represented, but zero is the most common digit, appearing 20,906,722 times. One is the least common digit, appearing only 20,889,487 times. By descending order of frequency, the digits are 0,5,2,3,8,4,7,6,9,1.
One other interesting tidbit: it's actually the conversion of the final result to decimal for output that takes the most time. If I comment out the call to mpz_out_str(), it finishes in 69 seconds.
This file is obviously a little large to post somewhere for download but if anybody would like to try confirming my results, the MD5 sum on the file is 10f1f75f783b7ae9802f38a6921ac561.