如何以素数为模计算非常大的二项式系数？

Question

This problem's answer turns out to be calculating large binomial coefficients modulo prime number using Lucas' theorem. Here's the solution to that problem using this technique: here.

现在我的问题是：

如果数据因变量溢出而增加，我的代码似乎会过期。有什么办法可以解决这个问题？
有没有什么方法可以不使用这个 theorem？

编辑： 请注意，由于这是一个 OI 或 ACM 问题，因此不允许使用原始库以外的外部库。

代码如下：

#include <iostream>
#include <string.h>
#include <stdio.h>
using namespace std;

#define N 100010

long long mod_pow(int a,int n,int p)
{
    long long ret=1;
    long long A=a;
    while(n)
    {
        if (n & 1)
            ret=(ret*A)%p;
        A=(A*A)%p;
        n>>=1;
    }
    return ret;
}

long long factorial[N];

void init(long long p)
{
    factorial[0] = 1;
    for(int i = 1;i <= p;i++)
        factorial[i] = factorial[i-1]*i%p;
    //for(int i = 0;i < p;i++)
        //ni[i] = mod_pow(factorial[i],p-2,p);
}

long long Lucas(long long a,long long k,long long p) 
{
    long long re = 1;
    while(a && k)
    {
        long long aa = a%p;long long bb = k%p;
        if(aa < bb) return 0; 
        re = re*factorial[aa]*mod_pow(factorial[bb]*factorial[aa-bb]%p,p-2,p)%p;
        a /= p;
        k /= p;
    }
    return re;
}

int main()
{
    int t;
    cin >> t;
    while(t--)
    {
        long long n,m,p;
        cin >> n >> m >> p;
        init(p);
        cout << Lucas(n+m,m,p) << "\n";
    }
    return 0;
}

Answer 1

无限精度整数似乎是可行的方法。

如果您使用 C++， PicklingTools 库有一个 "infinite precision" 整数（类似于 Python 的 LONG 类型）。别人推荐了Python，合理知道就回答Python。如果你想用 C++ 来做，你可以使用 int_n 类型：

#include "ocval.h"
int_n n="012345678910227836478627843";
n = n + 1;   // Can combine with other plain ints as well

查看文档：

http://www.picklingtools.com/html/usersguide.html#c-int-n-and-the-python-arbitrary-size-ints-long

和

http://www.picklingtools.com/html/faq.html#c-and-otab-tup-int-un-int-n-new-in-picklingtools-1-2-0

C++ PicklingTools 的下载地址是 here。

Answer 2

你想要一个bignum（a.k.a。任意精度算术）库。

首先，不要编写自己的 bignum（或 bigint）库，因为高效算法（比你在学校学到的朴素算法更高效）很难设计和实现。

那么，我会推荐GMPlib. It is free software, well documented, often used, quite efficient, and well designed (with perhaps some imperfections, in particular the inability to plugin your own memory allocator in replacement of the system malloc; but you probably don't care, unless you want to catch the rare out-of-memory condition ...). It has an easy C++ interface。它被打包在大多数 Linux 发行版中。

如果是家庭作业，也许你的老师希望你多思考数学，并通过一些证据找到解决问题的方法没有 bignums.

Answer 3

此解决方案假设 p² 适合 unsigned long long。由于 unsigned long long 按照标准至少有 64 位，这至少适用于 p 高达 40 亿，远远超过问题指定的数量。

typedef unsigned long long num;

/* x such that a*x = 1 mod p */
num modinv(num a, num p)
{
    /* implement this one on your own */
    /* you can use the extended Euclidean algorithm */
}

/* n chose m mod p */
/* computed with the theorem of Lucas */
num modbinom(num n, num m, num p)
{
    num i, result, divisor, n_, m_;

    if (m == 0)
        return 1;

    /* check for the likely case that the result is zero */
    if (n < m)
        return 0;

    for (n_ = n, m_ = m; m_ > 0; n_ /= p, m_ /= p)
        if (n_ % p < m_ % p)
            return 0;

    for (result = 1; n >= p || m >= p; n /= p, m /= p) {
        result *= modbinom(n % p, m % p, p);
        result %= p;
    }

    /* avoid unnecessary computations */
    if (m > n - m)
         m = n - m;

    divisor = 1;
    for (i = 0; i < m; i++) {
         result *= n - i;
         result %= p;

         divisor *= i + 1;
         divisor %= p;
    }

    result *= modinv(divisor, p);
    result %= p;

    return result;
}

Answer 4

假设我们需要计算 (a / b) mod p 的值，其中 p 是质数。由于 p 是素数，因此每个数字 b 都有一个倒数 mod p。所以(a / b) mod p = (a mod p) * (b mod p)^-1。我们可以使用欧几里得算法来计算逆。

要获得 (n over k)，我们需要计算 n! mod p、(k!)^-1、((n - k)!)^-1。总时间复杂度为 O(n).

更新： 这是 C++ 中的代码。不过我没有对其进行广泛的测试。

int64_t fastPow(int64_t a, int64_t exp, int64_t mod)                               
{                                                                                  
    int64_t res = 1;                                                               
    while (exp)                                                                    
    {                                                                              
        if (exp % 2 == 1)                                                          
        {                                                                          
            res *= a;                                                              
            res %= mod;                                                            
        }                                                                          

        a *= a;                                                                    
        a %= mod;                                                                  
        exp >>= 1;                                                                 
    }                                                                              
    return res;                                                                    
}                                                                                  

// This inverse works only for primes p, it uses Fermat's little theorem                                                                                   
int64_t inverse(int64_t a, int64_t p)                                              
{                                                                                  
    assert(p >= 2);                                                                
    return fastPow(a, p - 2, p);                                                   
}                                                                                  

int64_t binomial(int64_t n, int64_t k, int64_t p)                                  
{                                                                                  
    std::vector<int64_t> fact(n + 1);                                              
    fact[0] = 1;                                                                   
    for (auto i = 1; i <= n; ++i)                                                  
        fact[i] = (fact[i - 1] * i) % p;                                           

    return ((((fact[n] * inverse(fact[k], p)) % p) * inverse(fact[n - k], p)) % p);
}

如何以素数为模计算非常大的二项式系数？

How to calculate EXTREMELY big binomial coefficients modulo by prime number?

c++

algorithm

binomial-coefficients