Java 中计算大 n 和 k 值模 10^9 + 7 的二项式系数的乘法公式输出错误值

Question

我有一个作业，要求创建一个程序，可以计算给定任意 n、k 的二项式系数，使得 1<=k<=n<=2000. I am able to accomplish this for small n, k, but for some reason its not accurate at larger n and k values.

这是我当前的代码

import java.io.BufferedInputStream;
import java.util.Scanner;

public class Main
{
    public static double bCoeffModded(double nPar, double kPar, double mPar) {
        if ((nPar == kPar) || (kPar == 0)) {
            return 1;
        }
        
        double binomialCoefficient = nPar;
        
        for (int i = 2; i <= kPar; i++) {
            binomialCoefficient *= ((nPar + 1 - i) / i);
            binomialCoefficient %= mPar;
        }
        return binomialCoefficient;
    }
    
    public static void main(String[] args) {
        // Declarations
        double n, k;
        double m = (int) Math.pow(10,9) + 7; // modulus
        int c;
        
        // Input
        Scanner reader = new Scanner(new BufferedInputStream(System.in));
        n = reader.nextInt();
        k = reader.nextInt();
        
        // Calculations
        c = (int) bCoeffModded(n, k, m);
        
        // OutPut
        System.out.println(c);
    }
}

我在每次迭代中应用模数，因为它减少了溢出的机会。我认为无论如何都会存在溢出，因为准确性无论如何都会下降。

二项式系数的公式和概述可以在这里找到。我应该注意我已经尝试过递归方法和阶乘方法。两者也都存在溢出问题。

测试用例1：
6 4
预计15
实际15

测试用例2：
100 50
预计 538992043
实际309695578

Answer 1

嗯，正如我在评论中提到的，你不能将整数除法与模除法一起使用，否则你会得到错误的结果。这意味着您需要找到某种不涉及除法的循环关系或算法，我确实看到了以下选项：

我。明显的循环关系基于帕斯卡三角形：

(n,k) = (n-1,k-1) + (n-1,k)

它不涉及除法，因此我们可以使用以下想法进行计算：

mod((n,k),m) = mod((n-1,k-1),m) + mod((n-1,k),m)

问题是，我们可能达到的最佳时间复杂度是〜

O(n*k)

，而且，我们需要足够的内存来存储中间结果或调用堆栈。然而，我确实相信这种方法（带有一些明显的优化）是从程序员角度来看唯一相关的方法 - 其他方法需要一些研究或数学背景。

二.数学解答

根据定义：

(n,k)

=

n*(n-1)*..*(n-k-1)/(1*2*..*k)

这里我们确实知道如何计算

mod(n*(n-1)*..*(n-k-1),m)

，这似乎没什么大不了的：

mod(mod(mod(n,m)*mod(n-1,m),m)*mod(n-2,m),m)....

真正的挑战是如何处理

1/(1*2*..*k)

部分，主要思想是：是否可以用乘法代替

mod(p/q, m)

中的除法，即找到一些满足以下等式的

r

：

mod(p/q, m) = mod(p*r, m)

？答案是：是的，有可能 - 根据费马定理，它是

pow(q,m-2)

，技术上我们甚至可以用

mod(pow(q,m-2),m)

时间复杂度计算

O(logm)

：

public static int modpow(int x, int n) {
    int modulo = 1000000007;
    if (n == 0) {
        return 1;
    } else if (n == 1) {
        return x % modulo;
    } else if (n == -1) {
        return 1 / x;
    }
    int p = modpow(x, n >> 1);
    long product = ((long) p * p) % modulo;
    return (int) (product * modpow(x, n & 1) % modulo);
}

public static int ncr(int n, int r) {
    int modulo = 1000000007;

    long mult = 1;
    for (int i = n; i > n - r; i--) {
        mult = (mult * i) % modulo;
    }
    
    for (int i = 2; i <= r; i++) {
        mult = (mult * modpow(i, modulo - 2)) % modulo;
    }

    return (int) mult;
}

但是，正如@Dawood ibn Kareem 在评论中实际提到的那样：

BigInteger#modInverse

也返回一些有用的东西，我们可以使用以下形式重写以前的代码：

public static int ncr(int n, int r) {
    int modulo = 1000000007;
    BigInteger bm = BigInteger.valueOf(modulo);

    long mult = 1;
    for (int i = n; i > n - r; i--) {
        mult = (mult * i) % modulo;
    }

    for (int i = 2; i <= r; i++) {
        mult = (mult * BigInteger.valueOf(i).modInverse(bm).intValue()) % modulo;
    }

    return (int) mult;
}

另一方面，这太具体了，我们很可能需要“找到”另一种计算模逆的想法，而扩展欧几里得算法可以拯救：

function inverse(a, n)
    t := 0;     newt := 1
    r := n;     newr := a

    while newr ≠ 0 do
        quotient := r div newr
        (t, newt) := (newt, t − quotient × newt) 
        (r, newr) := (newr, r − quotient × newr)

    if r > 1 then
        return "a is not invertible"
    if t < 0 then
        t := t + n

    return t

在我们的例子中可以简化为以下循环关系（注意，因为

m % x < x

我们可以使用动态规划方法）：

inv(x) = (inv(m % x) * (m - m/x)) % m

现在我们得到：

public static int ncr(int n, int r) {
    int modulo = 1000000007;

    long mult = 1;
    for (int i = n; i > n - r; i--) {
        mult = (mult * i) % modulo;
    }

    int[] invs = new int[r + 1];
    invs[0] = 1;
    if (r > 1) {
        invs[1] = 1;
    }

    for (int i = 2; i <= r; i++) {
        invs[i] = (int) (((long) invs[modulo % i] * (modulo - modulo / i)) % modulo);
        mult = (mult * invs[i]) % modulo;
    }

    return (int) mult;
}