uoj86 mx的组合数 (lucas定理+数位dp+原根与指标+NTT)
题目描述自己看去吧(
题解时间
首先看到 $ 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();}