组合数学与计数 DP （二）

组合数学与计数 DP （二）

其实 FFT 也算组合数学。

Lucas 定理

给定 n,m,pn,m,p ，保证 pp 为质数，求 Cn+mm mod pC_{n+m}^m \bmod p 。

题目： Lucas 定理 。

Lucas 定理：

若 pp 为质数，则有：

Cnm mod p=C⌊np⌋⌊mp⌋Cn mod pm mod p mod pC_n^m \bmod p=C_{\lfloor \frac{n}{p}\rfloor}^{\lfloor \frac{m}{p}\rfloor}C_{n \bmod p}^{m \bmod p} \bmod p

证明：

引理 2 ：

对于质数 pp ，有 (1+x)p≡1+xp(modp)(1+x)^p \equiv 1+x^p \pmod p 。

令 n=ap+b,m=cp+dn=ap+b,m=cp+d 。

则有：

(1+x)n≡(1+x)ap(1+x)b≡(1+xp)a(1+x)b\begin{aligned} (1+x)^n &\equiv(1+x)^{ap}(1+x)^b\\ &\equiv(1+x^p)^a(1+x)^b \end{aligned}

然后对 (1+x)n(1+x)^n 进行二项式展开，观察 xmx^m 的系数 CnmC_n^m ：

Cnmxm≡Cacxcp×Cbdxd≡CacCbdxb\begin{aligned} C_n^mx^m &\equiv C_a^c x^{cp} \times C_b^d x^d &\equiv C_a^cC_b^dx^b \end{aligned}

所以有 Cnm≡CacCbd(modp)C_n^m \equiv C_a^cC_b^d \pmod p 。

所以 Lucas 定理成立。

所以我们就可以在 n,m<pn,m<p 时用阶乘算（因为此时逆元一定存在）， n>pn>p 或 m>pm>p 时递归处理即可。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=2000010;
int n,m,fac[N],inv[N],P,T;
int kpow(int a,int b){
    ll t=1;
    ll la=a;
    while(b){
        if(b&1)t=t*la%P;
        la=la*la%P;
        b>>=1;
    }
    return t;
}
int C(int n,int m,int p){
    if(n<m)return 0;
    if(n<p&&m<p){
        return 1ll*fac[n]*inv[n-m]%P*inv[m]%P;
    }
    return 1ll*C(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main(){
    cin>>T;
    while(T--){
        cin>>n>>m>>P;
        fac[0]=1;
        for(int i=1;i<=P-1;i++){
            fac[i]=1ll*fac[i-1]*i%P;
        }
        inv[P-1]=kpow(fac[P-1],P-2);
        for(int i=P-1;i>=1;i--){
            inv[i-1]=1ll*inv[i]*i%P;
        }
        cout<<C(n+m,m,P)<<"\n";
    }

    return 0;
}

P2480 [SDOI2010] 古代猪文

题目： [SDOI2010] 古代猪文 。

题目描述：

给定 n,gn,g ，求 (g∑d∣nCnd) mod 999911659(g^{\sum\limits_{d\mid n}C_n^d} )\bmod 999911659 ， n,g≤109n,g \le 10^9 。

原问题可以由扩展欧拉定理变为求 ∑d∣nCnd mod 999911658\displaystyle\sum_{d\mid n}C_n^d \bmod 999911658 。

这其中 dd 是可以枚举的，问题变成了求组合数，对一个合数取模。

不难想到把质因数分解 999911658999911658 后 CRT 合并，结果为 2,3,4679,356172,3,4679,35617 。

于是我们可以对于每个质因数用 Lucas 算一遍组合数后 CRT 合并得到答案。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#include <bits/stdc++.h>
using namespace std;
const int N=100100,P=999911659;
const int b[5]={0,2,3,4679,35617};
#define int long long
int n,g,fac[N],inv[N],a[N],m[N],t[N],M=P-1;
int kpow(int a,int b,int P){
    int t=1;
    while(b){
        if(b&1)t=t*a%P;
        a=a*a%P;
        b>>=1;
    }
    return t;
}
int C(int n,int m,int P){
    if(n<m)return 0;
    if(n<P&&m<P)return fac[n]*inv[m]%P*inv[n-m]%P;
    return C(n/P,m/P,P)*C(n%P,m%P,P)%P;
}
signed main(){
    cin>>n>>g;
    for(int i=1;i<=4;i++){
        fac[0]=1;
        for(int j=1;j<b[i];j++)fac[j]=fac[j-1]*j%b[i];
        inv[b[i]-1]=kpow(fac[b[i]-1],b[i]-2,b[i]);
        for(int j=b[i]-2;j>=0;j--)inv[j]=(j+1)*inv[j+1]%b[i];
        for(int j=1;j*j<=n;j++){
            if(n%j==0){
                (a[i]+=C(n,j,b[i]))%=b[i];
                if(j*j!=n){
                    (a[i]+=C(n,n/j,b[i]))%=b[i];
                }
            }
        }
    }
    int res=0;
    for(int i=1;i<=4;i++){
        m[i]=M/b[i];
        t[i]=kpow(m[i],b[i]-2,b[i]);
        (res+=m[i]*t[i]%M*a[i]%M)%=M;
    }
    if(g%P==0){
        cout<<"0\n";
        return 0;
    }
    cout<<kpow(g%P,res,P)<<"\n";
    return 0;
}

P4345 [SHOI2015] 超能粒子炮·改

题目： [SHOI2015] 超能粒子炮·改 。

题目描述：

求 ∑i=0kCni mod 2333\displaystyle\sum_{i=0}^k C_n^i \bmod 2333 。 n,k≤1018n,k \le 10^{18} 。

令 f(n,k)=∑i=0kCni mod 2333,p=2333f(n,k)=\displaystyle\sum_{i=0}^k C_n^i \bmod 2333,p=2333 ，则有：

f(n,k)=∑i=0kCni=∑i=0kC⌊np⌋⌊mp⌋Cn mod pm mod p=∑i=0⌊kp−1⌋C⌊np⌋i∑i=0p−1Cn mod pi+C⌊np⌋⌊kp⌋∑i=0k mod pCn mod pi=f(n mod p,p−1)f(⌊np⌋,⌊kp−1⌋)+C⌊np⌋⌊kp⌋f(n mod p,k mod p)\begin{aligned} f(n,k)&=\displaystyle\sum_{i=0}^k C_n^i \\ &=\displaystyle\sum_{i=0}^k C_{\lfloor \frac{n}{p}\rfloor}^{\lfloor \frac{m}{p}\rfloor}C_{n \bmod p}^{m \bmod p}\\ &=\displaystyle\sum_{i=0}^{\lfloor\frac{k}{p}-1\rfloor}C_{\lfloor\frac{n}{p}\rfloor}^i\displaystyle\sum_{i=0}^{p-1} C_{n\bmod p}^i + C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{k}{p}\rfloor}\displaystyle\sum_{i=0}^{k\bmod p}C_{n\bmod p}^i\\ &=f(n\bmod p,p-1)f(\lfloor\frac{n}{p}\rfloor,\lfloor\frac{k}{p}-1\rfloor)+C_{\lfloor\frac{n}{p}\rfloor}^{\lfloor\frac{k}{p}\rfloor}f(n\bmod p,k\bmod p) \end{aligned}

然后就可以愉快的 Lucas 求解了。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <bits/stdc++.h>
using namespace std;
const int N=10010,P=2333;
int fac[N],inv[N],C[N][N],f[N][N];
long long n,k;
int kpow(int a,int b){
    int t=1;
    while(b){
        if(b&1)t=1ll*t*a%P;
        a=1ll*a*a%P;
        b>>=1;
    }
    return t;
}
int Lucas(long long n,long long m){
    if(n<m)return 0;
    if(n<P&&m<P)return 1ll*fac[n]*inv[m]%P*inv[n-m]%P;
    return 1ll*Lucas(n%P,m%P)*Lucas(n/P,m/P)%P;
}
int F(long long n,long long k){
    if(n<P&&k<P)return f[n][k];
    int t1=F(n%P,P-1);
    int t2=F(n/P,k/P-1);
    int t3=1ll*Lucas(n/P,k/P)*F(n%P,k%P)%P;
    return (1ll*t2*t1%P+t3)%P;
}
signed main(){
    int T;
    C[0][0]=1;
    for(int i=1;i<P;i++){
        C[i][0]=C[i][i]=1;
        for(int j=1;j<i;j++){
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
        }
    }
    for(int i=0;i<P;i++){
        f[i][0]=C[i][0];
        for(int j=1;j<P;j++)
            f[i][j]=(f[i][j-1]+C[i][j])%P;
    }
    fac[0]=1;
    for(int i=1;i<P;i++)fac[i]=1ll*i*fac[i-1]%P;
    inv[P-1]=kpow(fac[P-1],P-2);
    for(int i=P-2;i>=0;i--)inv[i]=1ll*(i+1)*inv[i+1]%P;
    cin>>T;
    while(T--){
        cin>>n>>k;
        cout<<F(n,k)<<"\n";
    }
    return 0;
}

扩展中国剩余定理 exCRT

其实这东西不仅简单好写，而且比 CRT 更实用。

给定 nn 个同余方程：

{x≡a1(modb1)x≡a2(modb2)⋯x≡an(modbn)\begin{cases} x \equiv a_1 \pmod {b_1} \\ x \equiv a_2 \pmod {b_2} \\ \cdots \\ x \equiv a_n \pmod {b_n} \end{cases}

求满足条件的最小非负整数解。

这事实上和求这个等价：

{x≡a(modm1)x≡b(modm2)\begin{cases} x \equiv a \pmod{m_1} \\ x \equiv b \pmod{m_2} \\ \end{cases}

这是因为多个同余方程时我们可以一个一个的按照这个合并。

不妨设 x=m1p+a=m2q+bx=m_1p+a=m_2q+b ，则有 m1p−m2q=b−am_1p-m_2q=b-a 。

这个方程有解的充要条件是 gcd⁡(m1,m2)∣b−a\gcd(m_1,m_2)\mid b-a ，不过这不重要，重要的是这个可以用 exgcd 求解。不妨设满足条件的解为 c,dc,d ，则有 x≡m1c+a(modm1m2gcd⁡(m1,m2))x\equiv m_1c+a \pmod {\frac{m_1m_2}{\gcd(m_1,m_2)}} 。

题目： exCRT 。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define int __int128
int n;
ll tmp=0;
int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
void exgcd(int a,int b,int &x,int &y){
    if(b==0)return (void)(x=1,y=0);
    exgcd(b,a%b,x,y);
    swap(x,y);
    y-=a/b*x;
}
signed main(){
    cin>>tmp;n=tmp;
    int ans,Mod;
    for(int i=1;i<=n;i++){
        int x,y;
        cin>>tmp;y=tmp;
        cin>>tmp;x=tmp;
        if(i==1){
            ans=x;Mod=y;
        }
        else{
            int ansA=0,ansB=0,t=Mod*y/gcd(Mod,y);
            exgcd(Mod,y,ansA,ansB);
            ans=(((x-ans)/gcd(Mod,y)*Mod%t*ansA%t+ans)%t+t)%t;Mod=t;
        }
    }
    tmp=ans;
    cout<<tmp<<"\n";
    return 0;
}

exCRT 其实无论哪方面都比 CRT 好用，不知道 CRT 还有什么用。

exLucas

其实 exLucas 和 Lucas 也没什么关系，不过与 exCRT 和 CRT 关系不同的是， exLucas 极为难写，一般不写封装会几位难看。

给定 n,m,pn,m,p ，求 Cnm mod pC_{n}^m \bmod p 。 p≤106p \le 10^6 ， n,m≤1018n,m \le 10^{18} 。

事实上，由于 exCRT 的存在，我们总能将其等价为模某一质数的幂的意义下的运算。

现在考虑求 Cnm mod paC_n^m \bmod p^a 。

因为阶乘不可写的主要原因其实是模意义下没有逆元，导致 n!m!(n−m)!\frac{n!}{m!(n-m)!} 的公式用不了。

但是，如果设 g(n)g(n) 为 n!n! 中的 pp 的因子个数，那么我们就可以用 n!pg(n)m!pg(m)⋅(n−m)!pg(n−m)×pg(n)−g(m)−g(n−m) mod pa\displaystyle\frac{\frac{n!}{p^{g(n)}}}{\frac{m!}{p^{g(m)}}\cdot\frac{(n-m)!}{p^{g(n-m)}}} \times p^{g(n)-g(m)-g(n-m)} \bmod p^a 求解，因为这样模意义下就一定有逆元了。

考虑如何求解这个。

令 f(n)=n!pg(n)f(n)=\frac{n!}{p^{g(n)}} 。则有：

f(n)=⌊np⌋pg(⌊np⌋)∏i=0,p∤in mod pai×∏i=0,p∤ipa−1if(n)=\frac{\lfloor\frac{n}{p}\rfloor}{p^{g(\lfloor\frac{n}{p}\rfloor)}}\prod_{i=0,p\nmid i}^{n\bmod p^a}i\times \prod_{i=0,p\nmid i}^{p^a-1}i

f(n)=f(⌊np⌋)∏i=0,p∤in mod pai×∏i=0,p∤ipa−1if(n)=f(\lfloor\frac{n}{p}\rfloor)\prod_{i=0,p\nmid i}^{n\bmod p^a}i\times \prod_{i=0,p\nmid i}^{p^a-1}i

定义 ck=∏i=0,p∤ikic_k=\prod_{i=0,p\nmid i}^{k}i ， f(n)f(n) 的计算复杂度为 O(log⁡n)O(\log n) 。

我们还要考虑计算 g(n)g(n) 。

由上式容易看出， g(n)=g(⌊np⌋)+⌊np⌋g(n)=g(\lfloor\frac{n}{p}\rfloor)+\lfloor\frac{n}{p}\rfloor 。

所以我们可以在快速计算 Cnm mod paC_n^m \bmod p^a ，最后 exCRT 合并答案就行。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=1000010;
int n,m,p,IDX,a[N],b[N],c[N];
void exgcd(int a,int b,int &x,int &y){
    if(b==0)return (void)(x=1,y=0);
    exgcd(b,a%b,x,y);
    swap(x,y);
    y-=a/b*x;
}
int kpow(int a,int b,int P){
    int t=1;
    while(b){
        if(b&1)t=1ll*t*a%P;
        a=1ll*a*a%P;
        b>>=1;
    }
    return t;
}
void calc(int P,int PK){
    int num=1;
    c[0]=1;
    for(int i=1;i<PK;i++){
        if(i%P)num=1ll*num*i%PK;
        c[i]=num;
    }
}   
int F(int n,int P,int PK){
    if(n==0)return 1;
    return 1ll*F(n/P,P,PK)*kpow(c[PK-1],n/PK,PK)%PK*c[n%PK]%PK;
}
int G(int n,int P){
    if(n<P)return 0;
    return G(n/P,P)+n/P;
}
int inv(int x,int P){
    int ansA=0,ansB=0;
    exgcd(x,P,ansA,ansB);
    return (ansA%P+P)%P;
}  
int exLucas(int n,int m,int P,int PK){
    if(n<m)return 0;
    int x=F(n,P,PK),y=F(m,P,PK),z=F(n-m,P,PK);
    return 1ll*x*inv(y,PK)%PK*inv(z,PK)%PK*kpow(P,G(n,P)-G(m,P)-G(n-m,P),PK)%PK;
}
signed main(){
    cin>>n>>m>>p;
    int t=p;
    for(int i=2;t>1;i++){
        if(t%i==0){
            int num=i,tnum=1;
            while(t%i==0){
                t/=i;
                tnum*=i;
            }
            calc(num,tnum);
            IDX++;
            a[IDX]=exLucas(n,m,num,tnum);b[IDX]=tnum;
        }
    }
    int res=0;
    for(int i=1;i<=IDX;i++){
        (res+=1ll*a[i]*(p/b[i])%p*inv(p/b[i],b[i])%p)%=p;
    }
    cout<<res<<"\n";
    return 0;
}

P3726 [AHOI2017/HNOI2017] 抛硬币

题目： [AHOI2017/HNOI2017] 抛硬币 。

题目描述：

给定 a,b≤1015,k≤9,0≤a−b≤104a,b \le 10^{15},k\le 9,0\le a-b\le 10^4 ，求 ∑i=0b∑j=i+1aCajCbi mod 10k\displaystyle\sum_{i=0}^b\sum_{j=i+1}^aC_a^jC_b^i\bmod 10^k 。

引理 3 （范德蒙德卷积）：

∑i=0kCniCmk−i=Cn+mk\sum_{i=0}^k C_n^iC_m^{k-i}=C_{n+m}^k

证明可参照 mathstheorem 。

我们有：

∑i=0b∑j=i+1aCajCbi=∑i=0b∑j=0a−i−1CajCbi=∑k=0a−1∑i=0bCak−iCbi=∑k=0a−1Ca+bk\begin{aligned} \sum_{i=0}^b\sum_{j=i+1}^aC_a^jC_b^i &= \sum_{i=0}^b\sum_{j=0}^{a-i-1}C_a^{j}C_b^i\\ &=\sum_{k=0}^{a-1}\sum_{i=0}^{b}C_a^{k-i}C_b^i\\ &=\sum_{k=0}^{a-1} C_{a+b}^k \end{aligned}

于是问题被转化为了一个前缀和问题。

由于 0≤a−b≤1040\le a-b\le 10^4 ，也就是说 a≥a+b2a\ge\frac{a+b}{2} ，那么我们可以直接累加 2a+b−12^{a+b-1} 。注意当 a+ba+b 为偶数时要额外累加 Ca+ba+b22\frac{C_{a+b}^{\frac{a+b}{2}}}{2} 。

然后我们就只要求 ∑i=a+b2+1a−1Ca+bi\sum_{i=\frac{a+b}{2}+1}^{a-1} C_{a+b}^i ，同时又 a−a+b2=a−b2a-\frac{a+b}{2}=\frac{a-b}{2} ，所以枚举范围不超过 50005000 ，于是用 exLucas 求一下即可。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N=2000010;
int a,b,k;
void exgcd(int a,int b,int &x,int &y){
    if(b==0)return (void)(x=1,y=0);
    exgcd(b,a%b,x,y);
    swap(x,y);y-=a/b*x;
}
int kpow(int a,int b,int P){
    int t=1;
    while(b){
        if(b&1)t=1ll*t*a%P;
        a=1ll*a*a%P;
        b>>=1;
    }
    return t;
}
struct ExLucas{
    int P,c[N],PK;
    ExLucas(){P=PK=0;}
    ExLucas(int k){
        P=(k==1?5:2),PK=(k==0?512:(k==1?1953125:1024));
        int num=1;
        for(int i=0;i<PK;i++){
            if(i%P)num=1ll*num*i%PK;
            c[i]=num;
        }
    }
    int kpw(int a,int b){return kpow(a,b,PK);}
    int F(int n){
        if(n==0)return 1;
        return F(n/P)*kpw(c[PK-1],n/PK)%PK*c[n%PK]%PK;
    }
    int G(int n){
        if(n<P)return 0;
        return n/P+G(n/P);
    }
    int inv(int x){
        int ansA=0,ansB=0;
        exgcd(x,PK,ansA,ansB);
        return (ansA%PK+PK)%PK;
    }
    int exLucas(int n,int m){
        if(n<m)return 0;
        int t=kpw(P,G(n)-G(n-m)-G(m));
        if(t==0)return 0;
        int x=F(n),y=F(m),z=F(n-m);
        return x*inv(y)%PK*inv(z)%PK*t%PK;
    }
}A(0),B(1),D(2);
int exCRT(int a,int b,int P1,int P2){
    int M=P1*P2,ansA=0,ansB=0;
    exgcd(P1,P2,ansA,ansB);
    ansA=(ansA%M+M)%M;
    return (ansA*P1%M*(b-a)%M+a)%M;
}
int ans[N];
const int IM=1000000000;
const int IM2=512,IM5=1953125;
int C(int n,int m){return exCRT(A.exLucas(n,m),B.exLucas(n,m),IM2,IM5);}
signed main(){
    while(~scanf("%lld%lld%lld",&a,&b,&k)){
        int t=kpow(2,a+b-1,IM);
        for(int i=(a+b)/2+1;i<a;i++)(t+=C(a+b,i))%=IM;
        if((a+b)&1^1)(t+=exCRT(D.exLucas(a+b,(a+b)/2),B.exLucas(a+b,(a+b)/2),IM2<<1,IM5)/2)%=IM;
        if(a==b)
        t=(kpow(2,a+b-1,IM)-exCRT(D.exLucas(a+b,(a+b)/2),B.exLucas(a+b,(a+b)/2),IM2<<1,IM5)/2+IM)%IM;
        for(int i=0;i<k;i++)ans[i]=t%10,t/=10;
        reverse(ans,ans+k);
        for(int i=0;i<k;i++)printf("%lld",ans[i]);puts("");
    }
    return 0;
}