JOJ1595：Combinations（求解大组合数，Lucas定理）

题目传送门：http://acm.jlu.edu.cn/joj/showproblem.php?pid=1595

大概意思很简单，就是求100以内组合数的解，这里顺便说说大组合数及其求余吧。

关于Lucas定理及其推广可以去看陆师兄的BLOG -> 我是传送门

这里就只放代码了

#include <cstdio>
#include <cmath>
#include <memory.h>

typedef long long typec;
///Lib functions
typec GCD(typec a, typec b)
{
    return b ? GCD(b, a % b) : a;
}
typec extendGCD(typec a, typec b, typec& x, typec& y)
{
    if(!b) return x = 1, y = 0, a;
    typec res = extendGCD(b, a % b, x, y), tmp = x;
    x = y, y = tmp - (a / b) * y;
    return res;
}
///for x^k
typec power(typec x, typec k)
{
    typec res = 1;
    while(k)
    {
        if(k&1) res *= x;
        x *= x, k >>= 1;
    }
    return res;
}
///for x^k mod m
typec powerMod(typec x, typec k, typec m)
{
    typec res = 1;
    while(x %= m, k)
    {
        if(k&1) res *= x, res %= m;
        x *= x, k >>=1;
    }
    return res;
}
/***************************************
Inverse in mod p^t system
***************************************/
typec inverse(typec a, typec p, typec t = 1)
{
    typec pt = power(p, t);
    typec x, y;
    y = extendGCD(a, pt, x, y);
    return x < 0 ? x += pt : x;
}
/***************************************
Linear congruence theorem
x = a (mod p)
x = b (mod q)
for gcd(p, q) = 1, 0 <= x < pq
***************************************/
typec linearCongruence(typec a, typec b, typec p, typec q)
{
    typec x, y;
    y = extendGCD(p, q, x, y);
    while(b < a) b += q / y;
    x *= b - a, x = p * x + a, x %= p * q;
    if(x < 0) x += p * q;
    return x;
}
/***************************************
prime table
O(n)
***************************************/
const int PRIMERANGE = 1000000;
int prime[PRIMERANGE + 1];
int getPrime()
{
    memset (prime, 0, sizeof (int) * (PRIMERANGE + 1));
    for (int i = 2; i <= PRIMERANGE; i++)
    {
        if (!prime[i]) prime[++prime[0]] = i;
        for (int j = 1; j <= prime[0] && prime[j] <= PRIMERANGE / i; j++)
        {
            prime[prime[j]*i] = 1;
            if (i % prime[j] == 0) break;
        }
    }
    return prime[0];
}
/***************************************
get factor of n
O(sqrt(n))
factor[][0] is prime factor
factor[][1] is factor generated by this prime
factor[][2] is factor counter

need: Prime Table
***************************************/
///you should init the prime table before
int factor[100][3], facCnt;
int getFactors(int x)
{
    facCnt = 0;
    int tmp = x;
    for(int i = 1; prime[i] <= tmp / prime[i]; i++)
    {
        factor[facCnt][1] = 1, factor[facCnt][2] = 0;
        if(tmp % prime[i] == 0)
            factor[facCnt][0] = prime[i];
        while(tmp % prime[i] == 0)
            factor[facCnt][2]++, factor[facCnt][1] *= prime[i], tmp /= prime[i];
        if(factor[facCnt][1] > 1) facCnt++;
    }
    if(tmp != 1)
        factor[facCnt][0] = tmp, factor[facCnt][1] = tmp, factor[facCnt++][2] = 1;
    return facCnt;
}
/***************************************
Easy Combination for
C(n, k) that not exceeds limit of typec
****************************************/
typec combination(int n, int k)
{
    typec res = 1, g;
    if(k > n) return 0;
    if(n - k < k) k = n - k;
    for(int i = 0; i < k; i++)
    {
        if(res % (i + 1) == 0)
            res /= i + 1, res *= n - i;
        else if( (n - i) % (i + 1) == 0)
            res *= (n - i) / (i + 1);
        else
        {
            g = GCD(res, i + 1), res /= g;
            res *= (n - i) / ((i + 1) / g);
        }
    }
    return res;
}

/***************************************
C(n, k) mod m for k < 50,
O(k*k*lgk) m * m < typec_MAX
***************************************/
typec combinationModS(typec n, typec k, typec m)
{
    ///larger gate for more optimization
    ///too large gate may overflow somewhere
    const typec gate = 1LL << 50;
    if(k > n || m == 1) return 0;
    if(n - k < k) k = n - k;
    typec d[k],  tmp = 1, i = 0, j , h, g;
    for(i = 0, j = n - k + 1; n - i >= j; i++)
    {
        d[i] = n - i;
        while(gate / d[i] >= j && n - i != j)
            d[i] *= j, j++;
    }
    for( j = 2, h = k; j <= h; h--)
    {
        tmp *= h;
        while(gate /tmp >= j && h != j)
            tmp *= j, j++;
        for(int s = 0; tmp != 1; s++)
        {
            g = GCD(tmp, d[s]);
            d[s] /= g, tmp /= g;
        }
    }
    int res = 1;
    for(j = 0; j < i; j++)
        d[j] %= m, res *= d[j], res %= m;
    return res;
}
/*****************************************
C(n, k) mod m using prime table for sieving
O(n)  m * m < typecMAX
limited by primetable, n could not be too large
need: prime table
*****************************************/
typec combinationModPri (typec n, typec k, typec m)
{
    if(k > n || m == 1) return 0;
    typec result = 1, cnt = 0, temp;
    for(int i = 1; i < prime[0] && prime[i] <= n; i++)
    {
        temp = n, cnt = 0;
        while(temp)
            temp /= prime[i], cnt += temp;
        temp = n - k;
        while(temp)
            temp /= prime[i], cnt -= temp;
        temp = k;
        while(temp)
            temp /= prime[i], cnt -= temp;
        temp = prime[i];
        while(cnt)
        {
            if(cnt & 1)
                result *= temp, result %= m;
            temp *= temp, cnt >>= 1, temp %= m;
        }
        if(result == 0) return 0;
    }
    return result;
}
/*************************************************
C(n, k) mod m
O(k*lgm) m * m < typec_MAX
*************************************************/
typec combinationModN(typec n, typec k, typec m)
{
    if(k > n || m == 1) return 0;
    k = (n - k < k)? n - k : k;
    int pcnt = 0;
    typec a = 1, b = 1, x, y, g;
    typec pa = 1, pb = 1; ///may over flow
    for(int i = 1; i <= k; i++)
    {
        a *= n - i + 1, b *= k - i + 1;
        while( (g = GCD(a, m)) > 1) pa *= g, a /= g;
        while( (g = GCD(b, m)) > 1) pb *= g, b /= g;
        g = GCD(pa, pb), pa /= g, pb /= g;
        while(pa % m == 0) pa /= m, pcnt++;
        while(pb % m == 0) pb /= m, pcnt--;
        b %= m, a %= m;
    }
    a *= pa / pb, a %= m;
    while(pcnt) return 0;
    extendGCD(b, m, x, y);
    if(x < m) x += m;
    x *= a, x %= m;
    return x;
}

/*******************************************
C(n, k) mod p
O(k) p*p <= typecMAX
*******************************************/
typec combinationModP(typec n, typec k, typec p)
{
    if(k > n) return 0;
    if(n - k < k) k = n - k;
    typec a = 1, b = 1, x, y;
    int pcnt = 0;
    for(int i = 1; i <= k; i++)
    {
        x = n - i + 1, y = i;
        while(x % p == 0) x /= p, pcnt++;
        while(y % p == 0) y /= p, pcnt--;
        x %= p, y %= p, a *= x, b *= y;
        b %= p, a %= p;
    }
    if(pcnt) return 0;
    extendGCD(b, p, x, y);
    if(x < 0) x += p;
    a *= x, a %= p;
    return a;
}
/*******************************************
C(n, k) mod p^t
O(k*lgn/lgp) p^2t < typecMAX
*******************************************/
typec combinationModPt(typec n, typec k, typec p, typec t)
{
    if(k > n) return 0;
    if(n - k < k) k = n - k;
    typec pt = power(p, t);
    typec a = 1, b = 1, x, y;
    int pcnt = 0;
    for(int i = 1; i <= k; i++)
    {
        x = n - i + 1, y = i;
        while(x % p == 0) pcnt++, x /= p;
        while(y % p == 0) pcnt--, y /= p;
        x %= pt, y %= pt, a *= x, b *= y;
        a %= pt, b %= pt;
    }
    if(pcnt >= t) return 0;
    extendGCD(b, pt, x, y);
    if(x < 0) x += pt;
    a *= x, a %= pt;
    return a * power(p, pcnt) % pt;
}

/*******************************************
C(n, k) mod m
O(k*lgn/lgp) m * m < typecMAX
    p is fractor of m (depends on the smallest one)
need:
    prime table
    factor table
    combinationModPt()
    linearCongruence
*******************************************/
///you need to init the prime table
typec combinationModLi(typec n, typec k, typec m)
{
    if(k > n || m == 1) return 0;
    getFactors(m);
    typec a, b, p, q;
    for(int i = 0; i < facCnt; i++)
    {
        if(!i) a = combinationModPt(n, k, factor[i][0], factor[i][2]), p = factor[i][1];
        else b = combinationModPt(n, k, factor[i][0], factor[i][2]), q = factor[i][1];
        if(!i) continue;
        a = linearCongruence(a, b, p, q), p *= q;
    }
    return a;
}
/*******************************************
C(n, k) mod p
Lucas's theorem for combination mod p
O(p * lgn/lgp)
*******************************************/
typec lucas(typec n, typec k, typec p)
{
    typec res = 1;
    while (n && k && res)
    {
        res *= combinationModP(n % p, k % p, p);
        res %= p, n /= p, k /= p;
    }
    return res;
}
/*******************************************
a = n * (n - 1) * ... * (n - k + 1)
b = m * (m - 1) * ... * (m - k + 1)
c = ? from input
if a * c / b is an integer
this function will calculate this value module p^t
and the c input is moduled by p^t, be sure that gcd(c, p^t) = 1
O(len * lgn/lgp) , p^2t < typecMAX
*******************************************/
///the parameter &pcnt caches the factors consists of p
typec productQuotient(typec n, typec m, typec len, typec p, typec pt, typec &c, typec &pcnt)
{
    if(!c || n < len) return c = 0;
    typec &a = c, b = 1, x, y;
    for(int i = 1; i <= len; i++)
    {
        x = n - i + 1, y = m - i + 1;
        while(x % p == 0) x /= p, pcnt++;
        while(y % p == 0) y /= p, pcnt--;
        x %= pt, y %= pt, a *= x, b *= y;
        a %= pt, b %= pt;
    }
    extendGCD(b, pt, x, y);
    if(x < 0) x += pt;
    a *= x, a %= pt;
    return a;
}

/*******************************************
C(n, k) mod p^t
generalized Lucas's theorem for combination mod p
O(p^t * lgn/lgp), p^2t < typecMAX
*******************************************/
typec generalizedLucas(typec n, typec k, typec p, typec t)
{
    if(k > n) return 0;
    if(n - k < k) k = n - k;
    if(t == 1) return lucas(n, k, p);
    typec pt = power(p, t);
    typec c = 1, pcnt = 0, ktable[100], ntable[100], ltable[100];
    int cnt = 0;
    for(; n || k; cnt++)
    {
        ktable[cnt] = k, ntable[cnt] = n, ltable[cnt] = k % pt;
        n -= k % pt, k -= k % pt, n /= p, k /= p;
    }
    for(--cnt; c && cnt >= 0; cnt--)
        productQuotient(ntable[cnt], ktable[cnt], ltable[cnt], p, pt, c, pcnt);
    if(!c || pcnt >= t) return 0;
    return c * power(p, pcnt) % pt;
}
/*******************************************
C(n, k) mod m
O(min(k, p^t * lgn/lgp)) m * m < typecMAX
    p^t is fractor of m
need:
    prime table
    factor table
    generalizedLucas
    linearCongruence
*******************************************/
///you need to init the prime table
typec combinationModLucas(typec n, typec k, typec m)
{
    if(m == 1 || k > n) return 0;
    if(n - k < k) k = n - k;
    getFactors(m);
    typec a, b, p, q;
    for(int i = 0; i < facCnt; i++)
    {
        if(!i) a = generalizedLucas(n, k, factor[i][0], factor[i][2]), p = factor[i][1];
        else b = generalizedLucas(n, k, factor[i][0], factor[i][2]), q = factor[i][1];
        if(!i) continue;
        a = linearCongruence(a, b, p, q), p *= q;
    }
    return a;
}

/********************************************
********************************************/
const typec PTMAX = 10000;
typec facmod[PTMAX];
void initFacMod(typec p, typec t = 1)
{
    typec pt = power(p, t);
    facmod[0] = 1 % pt;
    for(int i = 1; i < pt; i++)
    {
        if(i % p) facmod[i] = facmod[i - 1] * i % pt;
        else facmod[i] = facmod[i - 1];
    }
}
///you should init the facmod[] before
typec factorialMod(typec n, typec &pcnt, typec p, typec t = 1)
{
    typec pt = power(p, t), res = 1;
    typec stepCnt = 0;
    while(n)
    {
        res *= facmod[n % pt], res %= pt;
        stepCnt += n / pt, n /= p, pcnt += n;
    }
    res *= powerMod(facmod[pt - 1], stepCnt, pt);
    return res %= pt;
}
typec combinationModPtFac(typec n, typec k, typec p, typec t = 1)
{
    if(k > n || p == 1) return 0;
    if(n - k < k) k = n - k;
    typec pt = power(p, t), pcnt = 0, pmcnt = 0;
    if(k < pt) return combinationModPt(n, k, p, t);
    initFacMod(p, t);
    typec a = factorialMod(n, pcnt, p, t);
    typec b = factorialMod(k, pmcnt, p, t);
    b *= factorialMod(n - k, pmcnt, p, t), b %= pt;
    pcnt -= pmcnt;
    if(pcnt >= t) return 0;
    a *= inverse(b, p, t), a %= pt;
    return a * power(p, pcnt) % pt;
}
typec combinationModFac(typec n, typec k, typec m)
{
    getFactors(m);
    typec a, b, p, q;
    for(int i = 0; i < facCnt; i++)
    {
        if(!i) a = combinationModPtFac(n, k, factor[i][0], factor[i][2]), p = factor[i][1];
        else b = combinationModPtFac(n, k, factor[i][0], factor[i][2]), q = factor[i][1];
        if(!i) continue;
        a = linearCongruence(a, b, p, q), p *= q;
    }
    return a;
}
/*********************************************************
C(n, k) mod m with genelizedLucas and combinationModFac
O(min(k, p^t))
*********************************************************/
typec combinationModPtAuto(typec n, typec k, typec p, typec t)
{
    if(t > 6) return combinationModPtFac(n, k, p, t);
    return generalizedLucas(n, k, p, t);
}
typec combinationModOdds(typec n, typec k, typec m)
{
    if(k > n || m == 1) return 0;
    getFactors(m);
    typec a, b, p, q;
    for(int i = 0; i < facCnt; i++)
    {
        if(!i) a = combinationModPtAuto(n, k, factor[i][0], factor[i][2]), p = factor[i][1];
        else b = combinationModPtAuto(n, k, factor[i][0], factor[i][2]), q = factor[i][1];
        if(!i) continue;
        a = linearCongruence(a, b, p, q), p *= q;
    }
    return a;
}