使用C中的黄金比率计算模m的第n个斐波那契数

Question

我正在尝试计算第n个fibonacci数，模10^9+7，其中n由用户输入。

我用黄金比例来计算斐波那契数。

以下代码生成正确的结果，直到n=43。但对于n>=44，phi超过10^9+7，我开始得到意想不到的结果。此外，如果去除模量，则n>=44给出了正确的结果。

#include <stdio.h>
#include <math.h>
long double mod=1000000007;
long double power(long double base, long long int expo)
{
    if(base==1  ||  expo==0)
        return 1;
    if(expo&1)
    {
        long double temp = power(base, expo>>1);
        return fmodl(base * fmodl(temp*temp, mod), mod);
    }
    else
    {
        long double temp=power(base, expo>>1);
        return fmodl(temp*temp,mod);
    }
}

int main(void) {
    // your code goes here
    long double phi = (1+powl(5, 0.5))/2;
    long double phi_cap = (1 - powl(5, 0.5))/2;
    long double root5 = powl(5, 0.5);
    long long int n;
    scanf("%lld",&n);
    long double ans = fmodl(  (power(phi, n) - power(phi_cap, n)) * power(root5,mod-2), mod);
    printf("%.0Lfn", ans);
    return 0;
}

为什么会发生这种情况？使用长双精度来存储无理数是错误的吗？

感谢

Answer 1

第n个Fibonacci数的闭式表达式是

F_n= φⁿ/√5 - ψⁿ/√5

其中

φ=1/2+√5/2Ş1.6180339887
ψ=1/2-√5/2=φ-1Ş0.6180339887

由于F_n中的被减法部分总是小于一半，我们可以通过向零取整（或者，向负无穷大取整，或者floor()取整，因为F_n对于n≥0是非负的）来计算F₂：

F_n=⌊φⁿ/√5⌋=地板（φⁿ>/√5）

如果我们对F_n模M，M≥2感兴趣，我们需要观察模运算对上述公式的影响。注意，表达式"expr MOD M"通常使用C.中的fmod(expr, M)计算

我们可以应用模乘法的性质：对于正实数a、b和m，（amodm）（bmodm-。地板操作不受影响。这里，a=φⁿ，b=1/√5。这意味着我们可以简化表达式

F_nMOD M=⌊φⁿ/√5⌋MOD M

进入

F_nMOD M=⌊（φⁿMOD（M√5））/√5⌋

在这里，我们可以应用模幂运算，注意这一点上的模不是通过整数M，而是通过M·√5。

换句话说，如果我们有一个函数，它使用浮点模的模幂运算来计算浮点值的整数幂，比如

double pow_mod(double base, unsigned int exponent, double modulus);

我们可以使用计算模m的第n个Fibonacci数F_n

#define PHI    1.61803398874989484820458683436563811772
#define SQRT5  2.236067977499789696409173668731276235441;
double fn_mod_m;
unsigned int n, m;
fn_mod_m = floor(pow_mod(PHI, n, SQRT5 * (double)m) / SQRT5);

从右到左的方法是模幂运算的一个很好的候选者。