uoj86 mx的组合数 (lucas定理+数位dp+原根与指标+NTT)

uoj86 mx的组合数 (lucas定理+数位dp+原根与指标+NTT)

uoj

题目描述自己看去吧(

题解时间

首先看到 $ p $ 这么小还是质数，第一时间想到 $ lucas $ 定理。

注意 $ lucas $ 定理的另外一种写法是将数转换为 $ p $ 进制后计算$ C_{n}^{m} = Pi C_{a_i}^{b_i} $

所以考虑对于 $ l-1 $ 和 $ r $ 各进行一次数位 $ dp $ 。

$ dp[i][j] $表示从低位起算到 $ i $ 位计算结果取模后为 $ j $ 且保证是在合法范围以内的方案数

$ dg[i][j] $ 表示从低位起算到 $ i $ 位计算结果取模后为 $ j $ 且不保证是在合法范围以内的方案数

转移方法：

对于计算到某一位 $ i $

$ n $ 已经给定，也就是说 $ b_i $ 已经确定

所以枚举 $ x $ 值在这一位对应的 $ a_i $ 设为 $ k $ ，设 $ C_{k}^{b_i}=g $

转移：

$ dg[i][j cdot g mod p]+=dg[i-1][j] $

$ dp[i][j cdot g mod p]+=dg[i-1][j] ( k < a_{i_{max}} ) $

$ dp[i][j cdot g mod p]+=dp[i-1][j] ( k = a_{i_{max}} ) $

时间复杂度$ p^{2}logn $。

这个暴力好像是有50分。

然后考虑优化。

(这么毒瘤咋考虑出来的啊)

上式中的 $ j cdot g $ 可以考虑优化掉。

这时就如毒瘤的数学题一样，我们看到p是质数，考虑直接用指标把它降维就好了。。。（啥？）

还是考虑上面的dp方程。

我们现在枚举到i位，用上面第一个转移式为例。

设 $ f[x] = sum limits_{ k = 0 }^{ p - 1 } [ C_{k}^{b_i} == x] $

那么转移式变成一个乘法卷积 $ dg^{'} [i] = sum limits_{j=0}^{p-1} sum limits_{g=0}^{p-1} dg[j] * f[g] * [j cdot g mod p == i] $

上指标之后$ dg^{'} [i] = sum limits_{j=0}^{p-1} sum limits_{g=0}^{p-1} dg[j] * f[g] * [ (ln[j]+ln[g]) mod phi(p) == i] $

然后上NTT。

注意0没有指标，求出其他答案之后用总数减一下就能求出0的答案了。

#include<bits/stdc++.h>
using namespace std;
typedef long long lint;
typedef __int128 llint;
template<typename TP>inline void read(TP &tar)
{
	TP ret=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){ret=(TP)ret*10+ch-'0';ch=getchar();}
	tar=ret*f;
}
namespace LarjaIX
{
const int N=70011,maxn=65536,P=30011,B=150;
const int mo=998244353,G=3;
lint fpow(lint a1,lint p1,lint m1);
void ntt(lint *f1,int tp);
int p,phi,len=1,g;llint n,l,r;
int rev[N];
int bitn[B],bitm[B],maxbit;
int fac[P],inv[P],facinv[P];
int ln[P];
int c[B][P];
lint ans[P];
lint f1[N],f2[N],dp[N],dg[N],dt[N];
lint wg[N],iwg[N];
void work(llint lim)
{
	memset(dp,0,sizeof(dp));
	memset(dg,0,sizeof(dg));
	memset(bitn,0,sizeof(bitn));
	for(int i=1;lim;i++) bitn[i]=lim%p,lim/=p,maxbit=max(maxbit,i);
	dp[1]=dg[1]=1;
	for(int b=1;b<=maxbit;b++)
	{
		memset(f1,0,sizeof(f1)),memset(f2,0,sizeof(f2));
		for(int i=1;i<p;i++) f1[ln[i]]=dg[i];
		for(int i=bitm[b];i<p;i++)if(c[b][i]) f2[ln[c[b][i]]]++;
		ntt(f1,1),ntt(f2,1);
		for(int i=0;i<len;i++) f2[i]=f1[i]*f2[i]%mo;
		ntt(f2,-1);
		memset(dg,0,sizeof(dg));
		for(int i=0;i<len;i++) (dg[fpow(g,i%phi,p)]+=f2[i])%=mo;
		memset(f2,0,sizeof(f2));
		for(int i=bitm[b];i<bitn[b];i++)if(c[b][i]) f2[ln[c[b][i]]]++;
		ntt(f2,1);
		for(int i=0;i<len;i++) f2[i]=f1[i]*f2[i]%mo;
		ntt(f2,-1);
		memset(dt,0,sizeof(dt));
		for(int i=0;i<len;i++) (dt[fpow(g,i%phi,p)]+=f2[i])%=mo;
		if(c[b][bitn[b]])for(int i=1;i<p;i++) (dt[c[b][bitn[b]]*i%p]+=dp[i])%=mo;
		memcpy(dp,dt,sizeof(dp));
	}
}
int pri[P];
bool gck(int i){for(int j=1;j<=pri[0];j++) if(fpow(i,phi/pri[j],p)==1) return 0;return 1;}
int maid()
{
//	freopen("sample.in","r",stdin);
//	freopen("u.out","w",stdout);
	read(p),read(n),read(l),read(r),phi=p-1;
	//-----------------------------------------------------------------------------------------------------
	fac[1]=fac[0]=inv[1]=facinv[1]=facinv[0]=1;
	for(int i=2;i<p;i++) fac[i]=fac[i-1]*i%p,inv[i]=inv[p%i]*(p-p/i)%p,facinv[i]=facinv[i-1]*inv[i]%p;
	//-----------------------------------------------------------------------------------------------------
	while(len<=p*2) len<<=1;
	for(int i=1;i<=len;i<<=1) wg[i]=fpow(G,(mo-1)/(i<<1),mo),iwg[i]=fpow(wg[i],mo-2,mo);
	for(int i=1;i<len;i++) rev[i]=(rev[i>>1]>>1)|((len>>1)*(i&1));
	//-----------------------------------------------------------------------------------------------------
	//just get the g and ln of p
	{
		int tmp=phi;
		for(int i=2;i*i<=tmp;i++)if(tmp%i==0)
		{
			pri[++pri[0]]=i;
			while(tmp%i==0) tmp/=i;
		}
		if(tmp!=1) pri[++pri[0]]=tmp;
		for(int i=1;i<p;i++) if(gck(i)){g=i;break;}
		tmp=1;
		for(int i=0;i<phi;i++) ln[tmp]=i,(tmp*=g)%=p;
	}
	//-----------------------------------------------------------------------------------------------------
	{
		llint tmp=n;
		for(int bi=1;tmp;bi++) bitm[bi]=tmp%p,tmp/=p,maxbit=bi;//every bit of n
	}
	//-----------------------------------------------------------------------------------------------------
	for(int i=1;i<128;i++)for(int j=bitm[i];j<p;j++) c[i][j]=fac[j]*facinv[j-bitm[i]]%p*facinv[bitm[i]]%p;
	//-----------------------------------------------------------------------------------------------------
	work(r);
	for(int i=1;i<p;i++) ans[i]=dp[i];
	work(l-1);
	for(int i=1;i<p;i++) (ans[i]+=mo-dp[i])%=mo;
	ans[0]=(r-l+1)%mo;
	for(int i=1;i<p;i++) (ans[0]+=mo-ans[i])%=mo;
	for(int i=0;i<p;i++) printf("%lld
",ans[i]);
	return 0;
}
lint fpow(lint a1,lint p1,lint m1)
{
	lint ret=1;
	while(p1)
	{
		if(p1&1ll) (ret*=a1)%=m1;
		(a1*=a1)%=m1,p1>>=1;
	}
	return ret;
}
void ntt(lint *f1,int tp)
{
	for(int i=0;i<len;i++) if(i<rev[i]) swap(f1[i],f1[rev[i]]);
	lint ilen=fpow(len,mo-2,mo);
	for(int i=1;i<len;i<<=1)
	{
		lint w0=~tp?wg[i]:iwg[i];
		for(int j=0;j<len;j+=(i<<1))
		{
			lint w=1;
			for(int k=0;k<i;k++,(w*=w0)%=mo)
			{
				lint w1=f1[j+k],w2=w*f1[j+k+i]%mo;
				f1[j+k]=(w1+w2)%mo,f1[j+k+i]=(w1-w2+mo)%mo;
			}
		}
	}
	if(tp==-1) for(int i=0;i<len;i++) (f1[i]*=ilen)%=mo;
}
}
int main(){return LarjaIX::maid();}