ejolson wrote: ↑Thu Jun 20, 2019 12:01 am
I just rounded up four distinctly different algorithms in C using the GNU Multi-precision Library.
This post presents four different C codes implementing four different algorithms for quickly computing the n-th Fibonacci number. The codes are
- method1.c -- optimized doubling: The optimized form of the doubling formula that was used in classic.bas, fibo.bas, fibompi.c, parallel.c, visual.bas and the PL/I code posted here.
- method2.c -- recursive doubling: The recursive implementation of the doubling formula implemented in hfibo.sb and the Python code used for the Dramatic Fibonacci Comedy of Python Interpreters. Note that the Python code in the GitHub repository is similar but has been memoized. While memoization makes the speed of the algorithm roughly on par with the other methods presented here, for simplicity the memoized version has not been considered in this roundup. Note that slightly modified versions of recursive doubling--some with and some without memoization--have been written in Haskell, Java, JavaScript and Scheme.
- method3.c -- golden ratio: An algebraic method based on taking powers of the golden ratio in the ring Q[sqrt(5)] that was used here in the pfibo.pas.
- method4.c -- mystery while loop: The iterative method used here in fib2.sb.
For reference and to check for possible memory leaks, the exact C codes used for the roundup were
Code: Select all
/* method1.c -- Compute the nth Fibonacci Number
Written June 19, 2019 by Eric Olson
This program uses the doubling formulas
F(2k) = F(k)[2F(k+1)-F(k)]
F(2k+1) = F(k+1)[2F(k+1)-F(k)]+(-1)^(k+1)
and
F(2k+1) = F(k)[2F(k)+F(k+1)]+(-1)^(k)
F(2k+2) = F(k+1)[2F(k)+F(k+1)]
The results of this code can also serve as a point of comparison
to implementations of the doubling formulas using other languages
that leverage GMP big numbers behind the scenes. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
static mpz_t a,b,t1;
void fibowork(int n){
if(n==0){
mpz_set_ui(a,0);
mpz_set_ui(b,1);
return;
}
fibowork(n/2);
if(n%2==0){
// [a,b]=[a*(2*b-a),b*(2*b-a)-(-1)^k]
mpz_add(t1,b,b); mpz_sub(t1,t1,a);
mpz_mul(a,a,t1); mpz_mul(b,b,t1);
if(n%4==0) mpz_sub_ui(b,b,1);
else mpz_add_ui(b,b,1);
} else {
// [a,b]=[a*(2*a+b)+(-1)^k,b*(2*a+b)]
mpz_add(t1,a,a); mpz_add(t1,t1,b);
mpz_mul(b,b,t1); mpz_mul(a,a,t1);
if(n%4==1) mpz_add_ui(a,a,1);
else mpz_sub_ui(a,a,1);
}
return;
}
void fibo(int n){
if(n<2){
mpz_set_ui(b,n);
return;
}
fibowork((n-1)/2);
if(n%2==0){
// b=b*(2*a+b)
mpz_add(t1,a,a); mpz_add(t1,t1,b); mpz_mul(b,b,t1);
} else {
// b=b*(2*b-a)-(-1)^k
mpz_add(t1,b,b); mpz_sub(t1,t1,a); mpz_mul(b,b,t1);
if(n%4==1) mpz_sub_ui(b,b,1);
else mpz_add_ui(b,b,1);
}
return;
}
#define bits 4000000
int main(){
mpz_init2(a,bits); mpz_init2(b,bits); mpz_init2(t1,bits);
for(;;){
printf("? ");
char buf[128];
fgets(buf,sizeof(buf),stdin);
int n=atoi(buf);
if(n==0) break;
fibo(n);
mpz_out_str(stdout,10,b);
printf("\n");
}
mpz_clear(t1); mpz_clear(b); mpz_clear(a);
return 0;
}
Code: Select all
/* method2.c -- Compute the nth Fibonacci Number
Written June 19, 2019 by Eric Olson
This program recursively uses the doubling formulas
F(2k) = F(k)[2F(k+1)-F(k)]
F(2k+1) = F(k+1)F(k+1)+F(k)F(k)
The results of this code can also serve as a point of comparison
to implementations of the doubling formulas using other languages
that leverage GMP big numbers behind the scenes. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
void hfibo(int n,mpz_t r){
mpz_t fk,fk1;
mpz_init(r);
if(n<=2){
mpz_set_ui(r,1);
} else {
int k=n/2;
hfibo(k,fk);
hfibo(k+1,fk1);
if(n%2==0){
mpz_add(r,fk1,fk1); mpz_sub(r,r,fk); mpz_mul(r,r,fk);
} else {
mpz_mul(r,fk,fk); mpz_mul(fk1,fk1,fk1); mpz_add(r,r,fk1);
}
mpz_clear(fk); mpz_clear(fk1);
}
}
int main(){
for(;;){
printf("? ");
char buf[128];
fgets(buf,sizeof(buf),stdin);
int n=atoi(buf);
if(n==0) break;
mpz_t b; hfibo(n,b);
mpz_out_str(stdout,10,b);
mpz_clear(b);
printf("\n");
}
return 0;
}
Code: Select all
/* method3.c -- Compute the nth Fibonacci Number
Written June 19, 2019 by Eric Olson
This program uses the golden ratio and the formula
F(n) = (phi^n - (-phi)^(-n))/sqrt(5)
The results of this code can also serve as a point of comparison
to implementations of the doubling formulas using other languages
that leverage GMP big numbers behind the scenes. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
typedef struct { mpz_t a,b; } ring5;
mpz_t t1,t2;
ring5 ring5mul(ring5 x,ring5 y){
ring5 r;
mpz_init(r.a); mpz_init(r.b);
mpz_mul(t1,x.a,y.a); mpz_mul(t2,x.b,y.b); mpz_mul_ui(t2,t2,5);
mpz_add(r.a,t1,t2);
mpz_mul(t1,x.a,y.b); mpz_mul(t2,x.b,y.a);
mpz_add(r.b,t1,t2);
return r;
}
ring5 goldmul(ring5 x,ring5 y){
ring5 r=ring5mul(x,y);
mpz_tdiv_q_2exp(r.a,r.a,1);
mpz_tdiv_q_2exp(r.b,r.b,1);
return r;
}
ring5 goldpow(ring5 x, int n){
if(n<=1){
ring5 r;
mpz_init(r.a); mpz_set(r.a,x.a);
mpz_init(r.b); mpz_set(r.b,x.b);
return r;
} else {
ring5 g1=goldpow(x,n/2);
ring5 g2=goldmul(g1,g1);
mpz_clear(g1.a); mpz_clear(g1.b);
if(n%2==0) {
return g2;
} else {
ring5 r=goldmul(x,g2);
mpz_clear(g2.a); mpz_clear(g2.b);
return r;
}
}
}
int main(){
ring5 phinum;
mpz_init(t1); mpz_init(t2);
mpz_init(phinum.a); mpz_set_ui(phinum.a,1);
mpz_init(phinum.b); mpz_set_ui(phinum.b,1);
for(;;){
printf("? ");
char buf[128];
fgets(buf,sizeof(buf),stdin);
int n=atoi(buf);
if(n==0) break;
ring5 z=goldpow(phinum,n);
mpz_out_str(stdout,10,z.b);
mpz_clear(z.a); mpz_clear(z.b);
printf("\n");
}
mpz_clear(phinum.b); mpz_clear(phinum.a);
mpz_clear(t2); mpz_clear(t1);
return 0;
}
Code: Select all
/* method4.c -- Compute the nth Fibonacci Number
Written June 19, 2019 by Eric Olson
This program uses a mystery iterative method, probably based on
taking powers of a matrix such as
F(n) = [1 1; 1 0]^n [1,2]
The results of this code can also serve as a point of comparison
to implementations of the doubling formulas using other languages
that leverage GMP big numbers behind the scenes. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
static mpz_t a,b,p,q,psq,qsq,twopq,ap,aq,bp,bq;
void mystery(int n){
int count=0;
mpz_set_ui(a,1);
mpz_set_ui(b,0);
mpz_set_ui(p,0);
mpz_set_ui(q,1);
while(n>0){
if(n%2==0){
mpz_mul(psq,p,p);
mpz_mul(qsq,q,q);
mpz_mul(twopq,p,q);
mpz_mul_si(twopq,twopq,2);
mpz_add(p,psq,qsq);
mpz_add(q,twopq,qsq);
n=n/2;
} else {
mpz_mul(bq,b,q);
mpz_mul(aq,a,q);
mpz_mul(ap,a,p);
mpz_mul(bp,b,p);
mpz_add(a,bq,aq);
mpz_add(a,a,ap);
mpz_add(b,bp,aq);
n=n-1;
}
}
}
#define bits 4000000
int main(){
mpz_init2(a,bits); mpz_init2(b,bits);
mpz_init2(p,bits); mpz_init2(q,bits);
mpz_init2(psq,bits); mpz_init2(qsq,bits);
mpz_init2(twopq,bits);
mpz_init2(ap,bits); mpz_init2(aq,bits);
mpz_init2(bp,bits); mpz_init2(bq,bits);
for(;;){
printf("? ");
char buf[128];
fgets(buf,sizeof(buf),stdin);
int n=atoi(buf);
if(n==0) break;
mystery(n);
mpz_out_str(stdout,10,b);
printf("\n");
}
mpz_clear(bq); mpz_clear(bp);
mpz_clear(aq); mpz_clear(ap);
mpz_clear(twopq);
mpz_clear(qsq); mpz_clear(psq);
mpz_clear(q); mpz_clear(p);
mpz_clear(b); mpz_clear(a);
return 0;
}
