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

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.

Now my questions are:

• Seems like my code expires if the data increases due to overflow of variables. Any ways to handle this?
• Are there any ways to do this without using this theorem?

EDIT: note that as this is an OI or ACM problem, external libs other than original ones are not permitted.

Code below:

#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

This solution assumes that p² fits into an unsigned long long. Since an unsigned long long has at least 64 bits as per standard, this works at least for p up to 4 billion, much more than the question specifies.

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 2

An infinite precision integer seems like the way to go.

If you are in C++, the PicklingTools library has an "infinite precision" integer (similar to Python's LONG type). Someone else suggested Python, that's a reasonable answer if you know Python. if you want to do it in C++, you can use the int_n type:

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

Take a look at the documentation at:

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

and

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

The download for the C++ PicklingTools is here.

Answer 3

You want a bignum (a.k.a. arbitrary precision arithmetic) library.

First, don't write your own bignum (or bigint) library, because efficient algorithms (more efficient than the naive ones you learned at school) are difficult to design and implement.

Then, I would recommend 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. It is packaged in most Linux distributions.

If it is a homework assignment, perhaps your teacher is expecting you to think more on the math, and find, with some proof, a way of solving the problem without any bignums.

Answer 4

Lets suppose that we need to compute a value of (a / b) mod p where p is a prime number. Since p is prime then every number b has an inverse mod p. So (a / b) mod p = (a mod p) * (b mod p)^-1. We can use euclidean algorithm to compute the inverse.

To get (n over k) we need to compute n! mod p, (k!)^-1, ((n - k)!)^-1. Total time complexity is O(n).

UPDATE: Here is the code in c++. I didn't test it extensively though.

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