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

Question

我正在尝试计算第 n 个斐波那契数模 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("%.0Lf\n", ans);
    return 0;
}

为什么会这样？用long double存储无理数有错吗？

谢谢

Answer 1

closed form expression for the n'th Fibonacci number 是

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

在哪里

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

由于 F_n 中的减法部分总是小于二分之一，我们可以计算 F_n via rounding 趋向于零（或者，趋向于负无穷大，或 floor()，因为 F_n 对于 n ≥ 0 是非负数）：

F_n = ⌊ φⁿ / √5 ⌋ = floor( φⁿ / √5 )

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

我们可以应用property of modular multiplication：对于正数a、b和m，（a mod m)(b mod m) = (a b mod m）。地板操作不受影响。这里，a = φⁿ，b = 1/√5。这意味着我们可以简化表达式

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

进入

F_n MOD M = ⌊ ( φⁿ MOD (M√5) ) / √5 ⌋

在这里，我们可以申请modular exponentiation，注意此时的模不是整数M，而是M·√5。

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

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

我们可以计算第n个斐波那契数F_n模m使用

#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);

right-to-left method 是此处模幂运算的绝佳候选者。