洛谷P4491组合数学+NTT_洛谷组合数学-CSDN博客

本文链接：https://blog.csdn.net/animalcoder/article/details/81267633

第一道组合数学+NTT

例题1.洛谷P4491 ，给NMS<=1e7

题意：N个位置每个位置可以染为M种染色，若位置恰有K种颜色出现S次，则获得愉悦度WK

问所有愉悦度之和

K的范围：0<=E<=min(M,N/S),设F(i)为位置恰有i种颜色出现S次,由加法原理

$ans=\sum_{i=0}^{E}wi*F(i)$

套路求F(i):从N个位置选iS个位置从M种颜色选i个颜色选出来的iS个位置进行可重集排列

$F(i)=\binom{n}{iS}\binom{m}{i}\frac{(iS)!}{(S!)^{i}} *G(i)$ ，G（i）为剩下的位置没有颜色出现S次

剩下没有颜色出现S次可以容斥=全体-有1种颜色出现S次+有2种颜色出现S次-...

$|\bigcap_{j=1}^{??}{\bar{A_{j}}}|=\sum_{j=0}^{E-i}(-1)^{j}\binom{n-iS}{jS}\binom{m-i}{j}\frac{(jS)!}{(S!)^{j}}*(m-i-j)^{n-iS-jS}$

$\small ans=\sum_{i=0}^{E}wi\binom{n}{iS}\binom{m}{i}\frac{(iS)!}{(S!)^{i}}\sum_{j=0}^{E-i}(-1)^{j}\binom{n-iS}{jS}\binom{m-i}{j}\frac{(jS)!}{(S!)^{j}}(m-i-j)^{n-iS-jS}$

$\small ans=$ 懒得写了= =截了个图-》

化简到最后，换个求和次序就是NTT卷积形式

就是说： $\small ans=\sum_{k}f(k)\sum_{i=0}^{k}a_{i}b_{k-i}=\sum_{k}f_{k}c_{k}$ ,ck为ai bi的卷积后第k项的系数


//如果 r* 2^k +1 是个素数，那么在modr*2^k+1意义下，可以处理 2^k 以内规模的数据，

//2281701377=17*2^27+1 是一个挺好的数，平方刚好不会爆 long long

//1004535809=479*2^21+1 加起来刚好不会爆 int 也不错

//还有就是 998244353=119*2^23+1

//2^21=2e6 
#include<bits/stdc++.h> 
#define ll long long
using namespace std;

inline int gi(){//快读 
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
    return f?x:-x;
}
const int _ = 1e7+5;
const int mod = 1004535809;//模数 
int n,m,s,E,lim,ans;int jc[_];
int inv[_],a[_],b[_],rev[_],l,og[_],len;//NTT本身用到的数组 
int qmod(int a,int b){//快速幂 
    int res=1;
    while (b) {if (b&1) res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}
    return res;
}
void ntt(int *P,int opt){//基本上不用改 
    for (int i=0;i<len;++i) if (i<rev[i]) swap(P[i],P[rev[i]]);
    for (int i=1;i<len;i<<=1){
        int W=qmod(3,(mod-1)/(i<<1));//3 根据模数的原根  可能要改 
        if (opt==-1) W=qmod(W,mod-2);
        og[0]=1;
        for (int j=1;j<i;++j) og[j]=1ll*og[j-1]*W%mod;
        for (int p=i<<1,j=0;j<len;j+=p)
            for (int k=0;k<i;++k){
                int x=P[j+k],y=1ll*og[k]*P[j+k+i]%mod;
                P[j+k]=(x+y)%mod,P[j+k+i]=(x-y+mod)%mod;
            }
    }
    if (opt==-1) for (int i=0,Inv=qmod(len,mod-2);i<len;++i) P[i]=1ll*P[i]*Inv%mod;
}
int main(){
    n=gi();m=gi();s=gi();//快读
	 
	E=min(m,n/s);        //求和表达式上限 (3.会用到) 
	lim=max(n,max(m,s)); //预处理上限(非NTT)
	
	//预处理阶乘 跟逆元阶乘 (非NTT)
    jc[0]=1; for (int i=1;i<=lim;++i) jc[i]=1ll*jc[i-1]*i%mod;
    inv[lim]=qmod(jc[lim],mod-2);
    for (int i=lim;i;--i) inv[i-1]=1ll*inv[i]*i%mod;
    
    //1.NTT   固定操作，处理数据的两倍 对2^k向上取整=len   len卷积上限 
    for (len=1;len<=(E<<1);len<<=1) ++l;--l;
    
    //2.NTT	  固定操作，预处理rev 蝶形变换 
    for (int i=0;i<len;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
    
    //3.NTT   灵活操作  ai bi的表达式 
    for (int i=0;i<=E;++i) a[i]=1ll*gi()*inv[i]%mod;
    for (int i=0;i<=E;++i) b[i]=i&1?mod-inv[i]:inv[i];
    
    //4.NTT   固定操作  ai卷积bi  也就是(a0+a1x+a2x^2...)*(b0+b1x+b2x^2....)=sum_{i=0...n}sum_{j=0...i}ai*b(j-i)
    ntt(a,1);ntt(b,1); 
    for (int i=0;i<len;++i) a[i]=1ll*a[i]*b[i]%mod;
    ntt(a,-1);
    
    //5.NTT   非固定操作  计算求和表达式的答案ans  
    for (int i=0;i<=E;++i) (ans+=1ll*qmod(m-i,n-i*s)*inv[m-i]%mod*qmod(inv[s],i)%mod*inv[n-i*s]%mod*a[i]%mod)%=mod;
    ans=1ll*jc[n]*jc[m]%mod*ans%mod;
    printf("%d\n",ans);
    return 0;
}