超快速计算斐波那契数列的某项的程序

傻蛋

程序原理是用矩阵乘法

@团长大人的团 @忆_碎碎念 @LarryGoode
可以删帖了

源代码
#include <stdio.h>
#include <gmp.h>
typedef mpz_t matrix4[2][2];
void initMatrix4(matrix4 m)
{
mpz_init( m[0][0] );
mpz_init( m[0][1] );
mpz_init( m[1][0] );
mpz_init( m[1][1] );
}
void freeMatrix4(matrix4 m)
{
mpz_clear( m[0][0] );
mpz_clear( m[0][1] );
mpz_clear( m[1][0] );
mpz_clear( m[1][1] );
}
void mulMatrix4(matrix4 &result, matrix4 left, matrix4 right)
{
mpz_t tmp0, tmp1;
mpz_init( tmp0 );
mpz_init( tmp1 );
// result[0][0] = left[0][0] * right[0][0] + left[0][1] * right[1][0];
mpz_mul(tmp0, left[0][0], right[0][0]);
mpz_mul(tmp1, left[0][1], right[1][0]);
mpz_add(result[0][0], tmp0, tmp1);
// result[0][1] = left[0][0] * right[0][1] + left[0][1] * right[1][1];
mpz_mul(tmp0, left[0][0], right[0][1]);
mpz_mul(tmp1, left[0][1], right[1][1]);
mpz_add(result[0][1], tmp0, tmp1);
// result[1][0] = left[1][0] * right[0][0] + left[1][1] * right[1][0];
mpz_mul(tmp0, left[1][0], right[0][0]);
mpz_mul(tmp1, left[1][1], right[1][0]);
mpz_add(result[1][0], tmp0, tmp1);
// result[1][1] = left[1][0] * right[0][1] + left[1][1] * right[1][1];
mpz_mul(tmp0, left[1][0], right[0][1]);
mpz_mul(tmp1, left[1][1], right[1][1]);
mpz_add(result[1][1], tmp0, tmp1);
mpz_clear( tmp0 );
mpz_clear( tmp1 );
}
void copyMatrix4(matrix4 &result, matrix4 source)
{
mpz_set(result[0][0], source[0][0]);
mpz_set(result[0][1], source[0][1]);
mpz_set(result[1][0], source[1][0]);
mpz_set(result[1][1], source[1][1]);
}
void powerMatrix4(matrix4 &result, matrix4 source, unsigned int n)
{
if (n == 0)
{
return;
}
if (n == 1)
{
copyMatrix4(result, source);
}
else
{
matrix4 tmp3, tmp4;
initMatrix4(tmp3);
initMatrix4(tmp4);
powerMatrix4(result, source, n / 2);
mulMatrix4(tmp3, result, result);
if (n & 1)
{
mulMatrix4(tmp4, tmp3, source);
copyMatrix4(result, tmp4);
}
else
copyMatrix4(result, tmp3);
freeMatrix4(tmp3);
freeMatrix4(tmp4);
}
}
void fib(mpz_t &result, unsigned int n)
{
matrix4 one, tmp;
if (n == 0)
{
mpz_set_ui(result, 0);
return;
}
initMatrix4( one );
initMatrix4( tmp );
mpz_set_ui(one[0][0], 1);
mpz_set_ui(one[0][1], 1);
mpz_set_ui(one[1][0], 1);
mpz_set_ui(one[1][1], 0);
powerMatrix4(tmp, one, n);
mpz_set(result, tmp[0][1]);
freeMatrix4( one );
freeMatrix4( tmp );
}
int main( void )
{
int n;
mpz_t r;
printf("Input number n:(0<=n<=100000000): ");
scanf("%u", &n);
if (n > 100000000)
return 0;
mpz_init(r);
fib(r, n);
gmp_printf("\nfib(%d) = %Zd\n", n, r);
mpz_clear(r);
return 0;
}

删帖准备

其实gmp库里就有计算斐波那契数列某项的函数，他使用的方法是
对于小数字，直接列表法，对于大数字用递归公式
F[2k] = F[k]*(F[k]+2F[k-1])
F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k
来计算，相对比矩阵法乘法数要少很多，
当然，大家可以看到，其实是从矩阵推导出来的

去找梅森素数吧，找到你发财了！