查找斐波那契数列中第 n 项的最后一位数字的问题

Question

我使用不错的博客 post 编写了一些代码来确定第 n 个 Fibonacci 数，该博客 post 在这个问题的已接受答案中给出：Finding out nth fibonacci number for very large 'n'。我这样做是为了练习在 projecteuler 上给出的更困难的递归问题，但这并不是真正相关的。该方法依赖于将问题更改为

形式的小型线性代数问题

Fn = T^n F1

其中 F1 = (1 1)^t 和 Fn 包含第 n 个和第 (n-1) 个斐波那契数。然后可以在 O(log n) 时间内确定 T^n 项。我成功地实施了这个并且它似乎工作正常。当我执行矩阵求幂时，我使用 %10000 所以我只得到最后 4 位数字，这似乎有效（我检查了一些大的斐波那契数）。但是，我想通过增加数字 10000 来尝试获得更多的最后一位数字。但这似乎不起作用。我不再得到正确的答案。这是我的代码

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

const unsigned long M = 10000;

unsigned long int * matProd(unsigned long int * A, unsigned long int * B){
  unsigned long int * C;
  C = malloc(4*sizeof(unsigned long int));
  C[0] = ((A[0]*B[0]%M) + (A[1]*B[2]%M)) % M;
  C[1] = ((A[0]*B[1]%M) + (A[1]*B[3]%M)) % M;
  C[2] = ((A[2]*B[0]%M) + (A[3]*B[2]%M)) % M;
  C[3] = ((A[2]*B[1]%M) + (A[3]*B[3]%M)) % M;
  return C;
}

unsigned long int * matExp(unsigned long int *A, unsigned long int n){
  if (n==1){
    return A;
  }
  if (n%2==0){
    return matExp(matProd(A,A),n/2);
  }
  return matProd(A,matExp(A,n-1));
}

unsigned long int findFib(unsigned long int n){
  unsigned long int A[4] = {0, 1, 1, 1};
  unsigned long int * C;
  C = malloc(4*sizeof(unsigned long int));
  C = matExp(A,n-2);
  return (C[2]+C[3]);
}

main(){
  unsigned long int n = 300;
  printf("%ld\n",findFib(n));
}

关于正确的编码约定和可以改进的地方，可能存在几个问题。我认为更改为 long int 可能会解决问题，但这似乎并不能解决问题。所以基本上问题是将 M 增加到例如 1000000 不会给我更多的数字，而是给我胡说八道。我犯了什么错误？

P.S。抱歉数学格式不好，我习惯了math.stackexchange。

Answer 1

问题可能是您运行在 long 大小为 32 位的系统上，我认为 Windows 就是这种情况。您可以通过编译和运行 printf("%d\n", sizeof(long)) 来检查它应该输出 4.

由于 M=1000000=10^6，小于 M 的两个数字的乘积可以达到 10^12，自 [=16] 以来计算矩阵条目时会遇到溢出问题=] 最多可以容纳 2^32-1 或大约 4 * 10^9.

要解决此问题，只需使用 unsigned long long 作为您的类型而不是 unsigned long。或者更好的是，uint64_t，保证在所有平台上都是 64 位的（并且需要 #include <stdint.h>）。这应该使您的代码适用于 M 到 sqrt(2^64)~10^9。如果您需要更大的整数，则需要使用大整数库。

Answer 2

如果程序对 M == 10000 有效但对 M == 1000000（甚至 M == 100000）失败，那么这可能意味着您的 C 实现的 unsigned long int 类型是 32 位宽。

如果您的矩阵元素完全取自 Z¹⁰⁰⁰⁰，那么它们最多需要 14 位有效二进制数字。因此，在减少模 M 之前，您在矩阵乘法函数中计算的乘积可能需要多达 28 个二进制数字。但是，如果将 M 增加到 100000，则矩阵元素最多需要 17 个二进制数字，而中间产品最多需要 34 个。减少模数 M 来不及防止溢出 a 32 位整数，因此给你垃圾结果。

您可以考虑将元素类型声明为 uint64_t。如果这是一个溢出问题，那么应该会给你足够的额外数字来处理 M == 1000000.

查找斐波那契数列中第 n 项的最后一位数字的问题

Issue finding last digits of the nth term in the Fibonacci sequence

c

fibonacci

modulus