EDIT: solved. ans var needed to be set to 0 in every iteration.That was a mistake I overlooked.A mpz_set_ui(ans,0); at the beginning of each loop solves the problem.Thanks to Oli Charlesworth for his effort!
I'm coding a fibonacci function in C using the GMP library,I use the matrix multiplication algorithm,but for some unexpected reason,I get wrong results(I strongly believe that the algorithm is right).It must have something to do with Gmp structs and inits.Any ideas?
/*********************************************************
Find n th Fibonacci number.Matrix algorithm.
+-- --+ +-- --+
|f(n+1) f(n) | | 1 1 | ^n
| | == | |
|f(n) f(n-1) | | 1 0 |
+-- --+ +-- --+
*********************************************************/
int
fastfib(mpz_t result,int n)
{
int cons[2][2]={{1,1},{1,0}}; //"constant" matrix
mpz_t mat[2][2]; //actual matrix holding fibonacci numbers
mpz_t ans; //holds f(n+1)
mpz_t return_val; //used for calculations
mpz_t temp; //used for calculations as well
//initialize (automatically set to zero)
mpz_init(ans);
mpz_init(return_val);
mpz_init(temp);
mpz_init(mat[0][0]);
mpz_init(mat[0][1]);
mpz_init(mat[1][0]);
mpz_init(mat[1][1]);
//start with n=1
mpz_set_ui(mat[0][0],1);
mpz_set_ui(mat[1][0],1);
mpz_set_ui(mat[0][1],1);
mpz_set_ui(mat[1][1],0);
//some trivial cases
if(n==0){
mpz_set_ui(result,0);
return 0;
}
if(n==1){
mpz_set_ui(result,1);
return 0;
}
if(n==2){
mpz_set_ui(result,1);
return 0;
}
n--;
while(n>1){
//fib[n+1]
//ans=mat[0][0]*cons[0][0]+mat[0][1]*cons[1][0];
mpz_set_ui(ans,0);
mpz_mul_ui(temp,mat[0][0],cons[0][0]);
mpz_add(ans,ans,temp);
mpz_mul_ui(temp,mat[0][1],cons[1][0]);
mpz_add(ans,ans,temp);
//update matrix
mpz_set(mat[1][1],mat[1][0]); //mat[1][1]=mat[1][0];
mpz_set(mat[1][0],mat[0][0]); //mat[0][1]=mat[1][0]=mat[0][0];
mpz_set(mat[0][1],mat[0][0]);
mpz_set(mat[0][0],ans); //mat[0][0]=ans;
n--;
}
//clear vars
mpz_clear(ans);
mpz_clear(return_val);
mpz_clear(temp);
mpz_set(result,mat[0][0]);
return 0;
}
Some debug info
results for n=2 // expected
2 1 // 2 1
1 1 // 1 1
results for n=3
5 2 // 3 2
2 1 // 2 1
results for n=4
12 5 // 5 3
5 2 // 3 2
results for n=5
29 12 // 8 5
12 5 // 5 3
results for n=6
70 29 // 13 8
29 12 // 8 5
70
solved. ans var needed to be set to 0 in every iteration.That was a mistake I overlooked.A mpz_set_ui(ans,0); at the beginning of each loop solves the problem.Thanks to Oli Charlesworth for his effort!