UOJ #86 mx的组合数 (数位DP+NTT+原根优化)

题目传送门

matthew99神犇的题解讲得非常清楚明白，跪烂Orzzzzzzzzzzzzz

总结一下，本题有很多重要的突破口

1.Lucas定理

看到n,m特别大但模数特别小时，容易想到$lucas$定理

$C_{n}^{m}=C_{n/p}^{m/p}\cdot C_{n\;mod\;p}^{m\;mod\;p}\;(mod\;p)$

但普通的$lucas$显然不适用于多次计算，我们可以把$lucas$定理展开

我们把$n$和$m$都看成两个$p$进制数$a$和$b$

$C_{n}^{m}=\pi C_{a_{i}}^{b_{i}}\;(mod\;p)$

这个展开显然成立

2.数位$DP$

想到了上一条进制转化，很容易就联想到数位$DP$

定义$f[i][j]$表示枚举到第$i$位，余数是$j$的方案数

转移十分经典，设现在要填上的数是$x$，$0$表示没达到上界，$1$达到上界，设$z=C_{x}^{b_{i+1}}$

$f0[i][j\cdot z\;mod\;p]+=f0[i][j]$

$x<a_{i+1}$时，$f0[i][j\cdot z\;mod\;p]+=f1[i][j]$

$x=a_{i+1}$时，$f1[i][j\cdot z\;mod\;p]+=f1[i][j]$

总复杂度$O(p^2log_{p}n)$

3.原根优化与多项式乘法

上面的$p^{2}$转移咋这么无脑呢，是不是有啥优化啊？是的

由于$p$是一个质数，它一定存在一个原根$g$

这就要涉及到离散对数了。其实这是一个映射

$g^{ind(x)}\equiv x\;(mod\;p)$

$ind(x)\;(ind(x)\in[0,p-1))$与$x\;(x\in [1,p))$是一组一一映射

那么$j\cdot z\;mod\;p$

$=g^{ind(j)+ind(z)\;(mod\;p-1)}\;mod\;p$

我们利用映射关系，把底数化成指数

这样转移变成了卷积的形式，用多项式乘法每次$O(plogp)$计算

计算出结果后，再逆映射回来得到实际的答案

而离散对数的映射并不能处理余数等于$0$的情况，我们每次$O(p)$单独讨论即可

总时间$O(plogplog_{p}n)$

 #include <bits/stdc++.h>

 #define N1 (1<<16)+10

 #define dd double

 #define ld long double

 #define ll long long

 #define uint unsigned int

 #define i128 __int128

 using namespace std;

 const int inf=0x3f3f3f3f;

 i128 gi128()

 {

     i128 ret=;int fh=;char c=getchar();

     while(c<''||c>''){if(c=='-')fh=-;c=getchar();}

     while(c>=''&&c<=''){ret=ret*+c-'';c=getchar();}

     return ret*fh;

 }

 ll qpow(ll x,ll y,const int &p)

 {

     ll ans=;

     for(;y;x=x*x%p,y>>=) if(y&) ans=ans*x%p;

     return ans;

 }

 const ll p=;

 namespace NTT{

 ll a[N1],b[N1],c[N1],Wn[N1],_Wn[N1];

 int r[N1];

 ll qpow(ll x,ll y)

 {

     ll ans=;

     for(;y;x=x*x%p,y>>=) if(y&) ans=ans*x%p;

     return ans;

 }

 void NTT(ll *s,int len,int type)

 {

     int i,j,k; ll wn,w,t;

     for(i=;i<len;i++) if(i<r[i]) swap(s[i],s[r[i]]);

     for(k=;k<=len;k<<=)

     {

         wn=(type>)?Wn[k]:_Wn[k];

         for(i=;i<len;i+=k)

         {

             for(j=,w=;j<(k>>);j++,w=w*wn%p)

             {

                 t=w*s[i+j+(k>>)]%p;

                 s[i+j+(k>>)]=(s[i+j]+p-t)%p;

                 s[i+j]=(s[i+j]+t)%p;

             }

         }

     }

 }

 void Pre(int len,int L)

 {

     int i;

     for(i=;i<len;i++) r[i]=(r[i>>]>>)|((i&)<<(L-));

     for(i=;i<=len;i<<=) Wn[i]=qpow(,(p-)/i), _Wn[i]=qpow(Wn[i],p-);

 }

 void Main(int len)

 {

     int i,invl=qpow(len,p-);

     NTT(a,len,); NTT(b,len,);

     for(i=;i<len;i++) c[i]=a[i]*b[i]%p;

     NTT(c,len,-);

     for(i=;i<len;i++) c[i]=c[i]*invl%p;

     memset(a,,sizeof(a)); memset(b,,sizeof(b));

 }

 };

 int T,m,G;

 i128 n,l,r;

 int a[N1],b[N1],tmp[N1];

 ll f0[][N1],f1[][N1],mul[N1],_mul[N1],ans[N1];

 inline ll C(int N,int M)

 {

     if(M>N) return ;

     return mul[N]*_mul[M]%m*_mul[N-M]%m;

 }

 int pr[N1],use[N1],ind[N1],_ind[N1],son[N1];

 void get_ind()

 {

     int i,j,ns=,flag,x,cnt=,sz=;

     for(i=;i<=m-;i++)

     {

         if(!use[i]) pr[++cnt]=i;

         for(j=;j<=cnt&&i*pr[j]<=m-;j++)

         {

             use[i*pr[j]]=;

             if(i%pr[j]==) break;

         }

     }

     for(j=;j<=cnt;j++) if((m-)%pr[j]==) son[++sz]=(m-)/pr[j];

     for(i=;i<=m-;i++)

     {

         flag=;

         for(j=;j<=sz;j++) if(qpow(i,son[j],m)==){ flag=; break;}

         if(flag){ G=i; break; }

     }

     ind[]=; _ind[]=; ind[]=m-;

     for(i=,x=;i<=m-;i++) x=x*G%m, ind[x]=i, _ind[i]=x;

 }

 //using NTT::a; using NTT::b; using NTT::c;

 void solve(i128 w,int type)

 {

     int i,j,k,nw=,nn=,len,L,now=,pst=;

     if(w<n){ ans[]=((w+)*type%p+ans[]+p)%p; return; }

     while(w>) nw++,tmp[nw]=w%m,w/=m;

     for(i=;i<=nw;i++) a[nw-i+]=tmp[i];

     i128 N=n;

     while(N>) nn++,tmp[nn]=N%m,N/=m;

     for(i=;i<=nn;i++) b[nw-i+]=tmp[i];

     for(len=,L=;len<m+m-;len<<=,L++);

     NTT::Pre(len,L);

     memset(f0,,sizeof(f0)); memset(f1,,sizeof(f1));

     for(j=;j<a[];j++) f0[][C(j,b[])]++;

     f1[][C(a[],b[])]++;

     for(i=;i<nw;i++)

     {

         memset(f0[now],,sizeof(f0[now])); memset(f1[now],,sizeof(f1[now])); 

         for(j=;j<m;j++) NTT::a[ind[j]]=f0[pst][j];

         for(j=;j<m;j++) NTT::b[ind[C(j,b[i+])]]++;

         for(j=;j<m;j++) (f0[now][]+=f0[pst][j]*NTT::b[m-])%=p;

         (f0[now][]+=f0[pst][]*m)%=p; NTT::b[m-]=;

         NTT::Main(len);

         for(j=;j<len;j++) (f0[now][_ind[j%(m-)]]+=NTT::c[j])%=p;

         for(j=;j<m;j++) NTT::a[ind[j]]=f1[pst][j];

         for(j=;j<a[i+];j++) NTT::b[ind[C(j,b[i+])]]++;

         for(j=;j<m;j++) (f0[now][]+=f1[pst][j]*NTT::b[m-])%=p;

         (f0[now][]+=f1[pst][]*a[i+])%=p; NTT::b[m-]=;

         NTT::Main(len);

         for(j=;j<len;j++) (f0[now][_ind[j%(m-)]]+=NTT::c[j])%=p;

         for(k=;k<m;k++)

         (f1[now][k*C(a[i+],b[i+])%m]+=f1[pst][k])%=p;

         swap(now,pst);

     }

     for(i=;i<m;i++) (ans[i]+=(f0[pst][i]+f1[pst][i])*type%p+p)%=p;

 }

 int main()

 {

     int i,j,x,cnt=;

     scanf("%d",&m); n=gi128(); l=gi128(); r=gi128(); l--;

     mul[]=mul[]=_mul[]=_mul[]=;

     for(i=;i<m;i++) mul[i]=mul[i-]*i%m, _mul[i]=qpow(mul[i],m-,m);

     get_ind();

     solve(r,);

     solve(l,-);

     for(i=;i<m;i++) printf("%lld\n",ans[i]);

     return ;

 }

巴特西

UOJ #86 mx的组合数 (数位DP+NTT+原根优化)

最新文章

热门文章