如何在模素数情况下计算极大的二项式系数？

Question

这个问题的答案是使用Lucas定理计算模质数的大二项式系数。以下是使用该技术解决该问题的解决方案：这里。

现在我的问题是：

• 似乎我的代码会因变量溢出而在数据增加时过期。有没有处理这种情况的方法？
• 是否有其他方法可以不使用定理来解决此问题？

编辑：请注���，由于这是一个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;
}

- Johnson Steward

---

我们在谈论一个“大”数字时需要多少位（链接固然好，但最好将所有相关信息都放在问题中）。如果任何一位都合理的话，您可以直接使用 http://www.boost.org/doc/libs/1_59_0/libs/multiprecision/doc/html/boost_multiprecision/tut/ints/cpp_int.html 并完成它。 - IdeaHat

你可以在中间结果上应用模算术，以避免溢出。你发布的解决方案使用mod_pow函数实现了这一点，并成功地使用了long long。 - M Oehm

如果不能使用外部库，请在问题中加入此说明，而不仅是发布评论。另外，binom.coeff.的n、k是否有限制？由于速度很重要，选择解决方案时可能需要了解这一点。 - deviantfan

@deviantfan 引用问题描述：1 <= n，m <= 1000000000，1 < p < 100000，并且保证p是一个质数。 - Johnson Steward

Lucas定理绝对是计算质数模下二项式系数的最佳方法。Lucas定理就是你要寻找的解决方案。 - fuz

Answer 1

这个解决方案假设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;
}

- fuz

---

太好了。这似乎消除了我所有的担忧。 - Johnson Steward

Answer 2

无限精度整数似乎是正确的选择。

如果您使用C++， PicklingTools库具有“无限精度”整数（类似于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 的下载链接在这里。

- rts1

---

请注意，您无法使用int_n表示大常量，因此必须使用字符串来表示常量。 - rts1

我希望无限精度的整数是C++的一部分：它们应该是的。 ;) - rts1

我也是。这样我们就可以节省很多编写（或背诵）自己的 MP 实现的时间，而将更多时间花在思考算法上。 - Johnson Steward

你可以使用整数向量（类似于十进制数字）来实现一个简单的无限精度。然后，你可以像高中时一样实现长乘法和加法（包括进位等）。 - rts1

另一个想法是，您可以因式分解需要相乘的所有小数字（并且您知道n choose k始终是int），然后划掉因子。剩下的任何东西可能适合于long long。 - rts1

Answer 3

假设我们需要计算 (a / b) mod p 的值，其中p是一个质数。由于p是质数，则每个数字b都具有模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);
}

- piotrekg2

Answer 4

你需要一个bignum（又称为任意精度算术）库。

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

然后，我建议使用GMPlib。它是免费软件，文档齐全，经常使用，相当高效，设计良好（可能存在一些缺陷，特别是不能插件自己的内存分配器来替换系统的malloc；但你可能不关心，除非你想捕捉罕见的内存不足情况...）。它有一个简单的C++接口。它已打包在大多数Linux发行版中。

如果这是一个家庭作业，也许你的老师希望你更多地考虑数学，并找到一种方法来解决问题没有使用bignums。

- Basile Starynkevitch

---

在普通的Linux中，我可以简单地包含gmpxx.h吗？ - Johnson Steward

@JohnsonSteward 是的。然后链接使用-lgmp。您可能需要先安装libgmp库。 - fuz

@FUZxxl的问题是我无法更改评估器使用的标志。 - Johnson Steward

@JohnsonSteward 好吧，那么你也不能使用libgmp。 - fuz