浅谈几种组合数的求值

rt，讨论一些见过的组合数求值

Case 1：

多次询问$C_n^m\ mod\ p$的值，$n,m \leq 5000$

Solution 1：

由递推式$C_n^m = C_{n-1}^m+C_{n-1}^{m-1}$预处理所有可能出现的询问
$O(n^2)预处理，O(1)$单次查询

Case 2：

多次询问$C_n^m\ mod\ p$的值，$n,m \leq 10^9,p \leq 10^5$且$p$为素数

Solution 2：

直接运用Lucas定理求解，不给出证明，列出一个比较好记忆的形式,
则有$Lucas(n,m)=Lucas(\lfloor \frac np \rfloor,\lfloor \frac mp \rfloor)\ast C_{n\ mod\ p}^{m\ mod\ p}$
对于左式，递归处理，对于右式，根据$C_n^m=\frac {n!}{m!(n-m)!}$，预处理阶乘及逆元即可

Case 3：

多次询问$C_n^m\ mod\ p$的值，$n,m \leq 10^5,1 \leq q \leq 10^6$

Solution 3：

对于模数不是素数的情形，Lucas定理就不再适用了，考虑直接从组合数的公式求解
由$C_n^m=\frac {n!}{m!(n-m)!}$，如果能预处理阶乘及逆元，则同样可以快速求解
但由于模数比较随意，所以逆元可能根本不存在
将模数质因数分解，即令$p = p_1^{a_1}p_2^{a_2}p_3^{a_3}\ldots p_k^{a_k}$
因为$p_k \leq 10^6$，所以本质不同的素因子不超过8个
对于每个阶乘，将它变成$p_1^{b_1}p_2^{b_2}p_3^{b_3}\ldots p_k^{b_k}d$的形式，其中$d$为一个常数
这样，做除法时特殊素因子部分等价于幂的减法操作
对于常数部分，因为不含和$p$相关的素因子，所以可以用扩展欧几里得算法求出逆元
预处理时使用历史版本递推，因此每个自然数只需要分解一次即可，
对一个数只用需要的素数分解，只需$O(8+logn)$，所以总复杂度$O(nlogn)$
单次询问需要求解逆元，复杂度$O(logn)$，当然也可以预处理好降为$O(1)$

代码如下：

#include<iostream>  
#include<cstdio>  
#include<algorithm>  
#include<cmath>  
#include<cstring>  
#include<vector>  
#include<queue>  
#include<set>  
#include<map>  
#include<stack>  
#include<bitset>  
#include<ext/pb_ds/priority_queue.hpp>  
using namespace std;  

const int maxn = 1E5 + 10;  
typedef long long LL;  

int n,mo,tp,Ans,p[10],k[10];  

int Mul(const LL &x,const LL &y) {return x * y % mo;}  
int Add(const int &x,const int &y) {return (x + y) % mo;}  

void Extend_Gcd(int a,LL &x,int b,LL &y)  
{  
    if (!b) {x = 1; y = 0; return;}  
    Extend_Gcd(b,y,a % b,x); y -= a / b * x;  
}  

int GetInv(int a)  
{  
    LL x,y;  
    Extend_Gcd(a,x,mo,y);  
    return (x % mo + mo) % mo;  
}  

struct Fac{  
    int a[10],b; Fac(){memset(a,0,sizeof(a)); b = 0;}  
    Fac operator * (const int &n)  
    {  
        Fac ret; int x = n;  
        for (int i = 1; i <= tp; i++)  
        {  
            ret.a[i] = a[i];  
            while (x % p[i] == 0) ++ret.a[i],x /= p[i];  
        }  
        ret.b = Mul(b,x);  
        return ret;  
    }  
    Fac operator /= (const Fac &B)  
    {  
        b = Mul(b,GetInv(B.b));  
        for (int i = 1; i <= tp; i++) a[i] -= B.a[i];  
    }  
}f[maxn];  

int C(int N,int M)  
{  
    Fac now = f[N];  
    now /= f[M]; now /= f[N - M];  
    int ret = now.b;  
    for (int i = 1; i <= tp; i++)  
        for (int j = 0; j < now.a[i]; j++)  
            ret = Mul(ret,p[i]);  
    return ret;  
}  

void Pre_Work()  
{  
    cin >> n >> mo; int x = mo;  
    if (mo == 1) {puts("0"); return 0;}  
    for (int i = 2; i <= mo; i++)  
    {  
        if (x % i != 0) continue; p[++tp] = i;  
        while (x % i == 0) ++k[tp],x /= i;  
        if (x == 1) break;  
    }  
    f[0].b = 1; for (int i = 1; i < n; i++) f[i] = f[i - 1] * i;  
}  

Case 4：

单次询问$C_n^m\ mod\ p$的值，$n,m\leq 10^{18},p\leq 10^6$

Solution 4：

由于$n,m$太大，已经无法想上面那样预处理阶乘和逆元了
不过还是可以用扩展Lucas定理求解
先将$p$质因数分解成$p=p_1^{a_1}p_2^{a_2}p_3^{a_3}\ldots p_k^{a_k}$
假设对于每个$p_i$都能统计出$C_n^m\ mod\ p_i^{a_i}$的值，
那么就可以使用合并线性同余方程组的方式求出答案了
从$C_n^m = \frac {n!}{m!(n-m)!}$下手，给出一种求解$K!\ mod\ p_i^{a_i}$的方式
举个例子，令$K = 19,p_i = 3,a_i = 2$
$19! = (1\ast 2\ast 4\ast 5\ast \ldots \ast 16\ast 17\ast 19)\ast (3\ast 6\ast \ldots \ast 15\ast 18)$
$\ \ \ \ \ = (1\ast 2\ast 4\ast 5\ast \ldots \ast 16\ast 17\ast 19)\ast 3^6\ast (6!)$
对于左边的部分，
由于$(1\ast 2\ast 4\ast 5\ast 7\ast 8) = (10\ast 11\ast 13\ast 14\ast 16\ast 17)\ (mod\ 3^2)$
可以发现，就是每$p_i^{a_i}$项出现一次循环(是$p$的倍数的项也参与计数但对答案无贡献)
而剩下一小段不在循环的部分，也就是$19$，肯定不超过$p_i^{a_i}$项，暴力统计即可
对于右边的部分，递归处理，对于中间的部分，将所有结果综合
于是最终答案就能构造成$p_i^t\ast d$的形式
这样，常数部分一定不含质因子$p_i$，就可以用扩展欧几里得算法求出逆元了
上述算法在极端情况下($p$本身是个大素数)由于递归单层时需要暴力求解的项比较多
将会达到最坏的复杂度$O(plogn)$

代码如下：

#include<iostream>  
#include<cstdio>  
#include<algorithm>  
#include<cmath>  
#include<cstring>  
#include<vector>  
#include<queue>  
#include<set>  
#include<map>  
#include<stack>  
#include<bitset>  
#include<ext/pb_ds/priority_queue.hpp>  
using namespace std;  

typedef long long LL;  

LL n,m,k,r1,p1;  

int Ksm(int x,int y)  
{  
    int ret = 1;  
    for (; y; y >>= 1)  
    {  
        if (y & 1) ret = ret * x;  
        x = x * x;  
    }  
    return ret;  
}  

LL ksm(LL x,LL y,LL mo)  
{  
    LL ret = 1;  
    for (; y; y >>= 1LL)  
    {  
        if (y & 1LL) ret = ret * x % mo;  
        x = x * x % mo;  
    }  
    return ret;  
}  

void Extend_Gcd(LL a,LL &x,LL b,LL &y)  
{  
    if (!b) {x = 1; y = 0; return;}  
    Extend_Gcd(b,y,a % b,x); y -= a / b * x;  
}  

LL GetInv(LL a,LL b)  
{  
    LL x,y;  
    Extend_Gcd(a,x,b,y);  
    return (x % b + b) % b;  
}  

LL Calc(LL N,LL tmp,LL p,LL x)  
{  
    if (!N) return 1;  
    LL y = N / x,now = ksm(tmp,y,x);  
    for (LL i = x * y + 1; i <= N; i += 1LL)  
        if (i % p != 0) now = i % x * now % x;  
    return now * Calc(N / p,tmp,p,x) % x;  
}  

LL Solve(LL p,LL t,LL x)  
{  
    LL ap = 0;  
    for (LL tmp = n; tmp; tmp /= p) ap += tmp / p;  
    for (LL tmp = k; tmp; tmp /= p) ap -= tmp / p;  
    for (LL tmp = n - k; tmp; tmp /= p) ap -= tmp / p;  
    LL tmp = 1;  
    for (LL i = 1; i <= x; i++)  
        if (i % p != 0) tmp = tmp * i % x;  
    LL a = Calc(n,tmp,p,x);  
    LL b = Calc(k,tmp,p,x);  
    LL c = Calc(n - k,tmp,p,x);  
    b = GetInv(b,x); c = GetInv(c,x);  
    return (a * b % x) * (c * ksm(p,ap,x) % x) % x;  
}  

int main()  
{  
    #ifdef DMC  
        freopen("DMC.txt","r",stdin);  
    #endif  

    cin >> n >> k >> m; int x = m;  
    for (int i = 2; x > 1; i++)  
    {  
        if (x % i != 0) continue;  
        int t = 0; while (x % i == 0) x /= i,++t;  
        LL p2 = Ksm(i,t),r2 = Solve(i,t,p2);  
        if (!p1) {r1 = r2; p1 = p2; continue;}  
        LL k1,Inv = GetInv(p1,p2);  
        k1 = (r2 - r1) * Inv % p2;  
        r1 = k1 * p1 + r1;  
        p1 = p1 * p2; r1 = (r1 % p1 + p1) % p1;  
    }  
    cout << r1 << endl;  
    return 0;  
}  

Case 5：

多次询问$C_n^m\ mod\ 10^9$的值，$n,m \leq 10^{15}$

Solution 5：

这里的模数$10^9$只是随便举得例子。。。
注意到$10^9 = 2^9\ast 5^9$
如果用扩展Lucas定理求解，极端情况下单次询问为$O(5^9log10^9)$
瓶颈在于，每次都需要大力统计$O(5^9)$的余项
如果能把这个暴力的部分去掉，就能快速计算了
注意到$5^9$其实是存得下的。。没错，可以先预处理好
其它部分和扩展Lucas定理的写法就大同小异了
所以说$10^9$是瞎举得例子
总得来说是固定模数且所有质因子分解后最大那个并不很大的都可以这样处理

---

这篇blog可能会持续更新的？？？
如果碰到新的题型的话