自动搬运

查看原文

来自洛谷，原作者为

	command_block 众水皆昂首，饮月唯我一。

搬运于2025-08-24 22:10:21，当前版本为作者最后更新于2019-07-11 18:19:06，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

本题小常数倍增NTT$O(nlogn)$解法!

虽然本蒟蒻人傻常数大跑了800+ms

但是作为一个求逆模板没有进入最优解前1000的人,这已经很快乐了

安利斯特林数小结,欢迎各位大佬前来吊打。

没有什么优美的公式帮忙了……我们来看看递推公式吧:

$\begin{Bmatrix}n \\ m\end{Bmatrix}=\begin{Bmatrix}n-1 \\ m-1\end{Bmatrix}+m \begin{Bmatrix}n-1 \\ m\end{Bmatrix}$

我们还是按套路,把答案设成多项式$H_k=\sum\limits_{i=0}\begin{Bmatrix} i \\k \end{Bmatrix}x^i$。

然后,我们发现$H_k=\sum\limits_{i=0}\left(\begin{Bmatrix}i-1 \\ k-1\end{Bmatrix}+k \begin{Bmatrix}i-1 \\ k\end{Bmatrix}\right)x^i$

$=\sum\limits_{i=0}\begin{Bmatrix}i-1 \\ k-1\end{Bmatrix}x^i+k\sum\limits_{i=0}\begin{Bmatrix}i-1 \\ k\end{Bmatrix}x^i$

$=x\sum\limits_{i=0}\begin{Bmatrix}i \\ k-1\end{Bmatrix}x^i+kx\sum\limits_{i=0}\begin{Bmatrix}i \\ k\end{Bmatrix}x^i=xH_{k-1}+kxH_k$

总结一下就是$H_k=xH_k+kxH_{k-1}$,那么$H_k=\dfrac{x}{1-kx}H_{k-1}$

根据递推公式能得到边界就是$H_0=1$,这里不再赘述。

那么,$H_k=H_0\prod\limits_{i=1}^k\dfrac{x}{1-ix}=x^k\left(\prod\limits_{i=1}^k(1-ix)\right)^{-1}$

重点就是怎么算$\prod\limits_{i=1}^k(1-ix)$。

$\prod\limits_{i=1}^k(1-ix)=\prod\limits_{i=1}^kx(\dfrac{1}{x}-i)=x^k\prod\limits_{i=1}^k(x^{-1}-i)$

使用一些小技巧,$\prod\limits_{i=1}^k(x-i)$翻转过来之后就是$x^k\prod\limits_{i=1}^k(x^{-1}-i)$

那么注意到$\prod\limits_{i=1}^k(x-i)=\dfrac{x^{\underline{k+1}}}{x}$我们只要求出$x^{\underline{k+1}}$就好。

咋整呢?我们还是考虑像这样来倍增(Orz rqy!):

我们知道$x^{\underline{2n}}=x^{\underline n}*(x-n)^{\underline n}$

我们知道一个多项式$F(x)$可以快速求出$F(x-c)$,方法:

首先$c=mod-c$,那么问题变成$F(x+c)$

设$F(x)$的系数为$f[0...n]$

$F(x+c)=\sum\limits_{i=0}^nf[i](x+c)^i$,显而易见吧!

$=\sum\limits_{i=0}^nf[i]\sum\limits_{j=0}^i\dbinom{i}{j}x^jc^{i-j}$(二项式定理)

$=\sum\limits_{i=0}^nf[i]\sum\limits_{j=0}^i\dfrac{i!}{j!(i-j)!}x^jc^{i-j}$

$=\sum\limits_{i=0}^nf[i]i!\sum\limits_{j=0}^i\dfrac{x^j}{j!}\dfrac{c^{i-j}}{(i-j)!}$

$=\sum\limits_{j=0}^n\dfrac{x^j}{j!}\sum\limits_{i=j}^nf_ii!\dfrac{c^{i-j}}{(i-j)!}$(交换求和号)

我们要求的就是$g[j]*j!=\sum\limits_{i=j}^nf[i]i!\dfrac{c^{i-j}}{(i-j)!}$

我们设$p[n]=f[n]n!;\ h[n]=\dfrac{c^{n}}{n!}$,则$g[j]j!=\sum\limits_{i=j}^np[i]h[i-j]$

但是我们发现$i+i-j≠$定值 , 我们把$p$翻转,得到$g[j]j!=\sum\limits_{i=j}^np'[n-i]h[i-j]$

那么$n-i+i-j=n-j$这就和$i$无关了。

卷起来之后翻转即可,记得最后乘上$j!$。

复杂度是$T(n)=T(n/2)+O(nlogn)=O(nlogn)$。

这东西跑得比EXP要快的多了,但是思考过程肯定没有EGF+EXP那样简洁。

于是乎现在可以考虑卡O(nlog^2n)的错解了(雾)

Code:

#include<algorithm>
#include<cstdio>
#define mod 167772161
#define G 3
#define Maxn 270000
using namespace std;
int n,k,r[Maxn<<2];
long long invn,invG;
long long fac[Maxn],inv[Maxn];
long long powM(long long a,long long t=mod-2)
{
  long long ans=1,buf=a;
  while(t){
  	if(t&1)ans=(ans*buf)%mod;
  	buf=(buf*buf)%mod;
  	t>>=1;
  }return ans;
}
void NTT(long long *f,bool op,int n)
{
  for (int i=0;i<n;i++)
    if (r[i]<i)swap(f[r[i]],f[i]);
  for (int len=1;len<n;len<<=1){
    int w=powM(op==1 ? G:invG,(mod-1)/len/2);
  	for (int p=0;p<n;p+=len+len){
  	  long long buf=1;
  	  for (int i=p;i<p+len;i++){
  	  	int sav=f[i+len]*buf%mod;
  	  	f[i+len]=f[i]-sav;
  	  	if (f[i+len]<0)f[i+len]+=mod;
  	  	f[i]=f[i]+sav;
  	  	if (f[i]>=mod)f[i]-=mod;
  	  	buf=buf*w%mod;
	    }//F(x)=FL(x^2)+x*FR(x^2)
	     //F(W^k)=FL(w^k)+W^k*FR(w^k)
	     //F(W^{k+n/2})=FL(w^k)-W^k*FR(w^k)
    }
  }
}
long long g[Maxn<<2];
void rev(long long *f,int len)
{
  for (int i=0;i<len;i++)g[i]=f[i];
  for (int i=0;i<len;i++)f[len-i-1]=g[i];
}
//令f=f*g (mod x^lim) 
void times(long long *f,long long *gg,int len,int lim)
{
  int m=len+len,n;
  for(int i=0;i<len;i++)g[i]=gg[i];
  for(n=1;n<m;n<<=1);invn=powM(n);
  for(int i=len;i<n;i++)g[i]=0;
  for(int i=0;i<n;i++)
    r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
  NTT(f,1,n);NTT(g,1,n);
  for(int i=0;i<n;++i)f[i]=(f[i]*g[i])%mod;
  NTT(f,0,n);
  for(int i=0;i<lim;++i)f[i]=(f[i]*invn)%mod;
  for(int i=lim;i<n;++i)f[i]=0;
}
void Init(int lim)
{
  inv[1]=inv[0]=fac[0]=1;
  for (int i=1;i<=lim;i++)fac[i]=fac[i-1]*i%mod;
  for (int i=2;i<=lim;i++)
  	inv[i]=inv[mod%i]*(mod-mod/i)%mod;
  for (int i=2;i<=lim;i++)inv[i]=inv[i-1]*inv[i]%mod;
  for (int i=1;i<=lim;i++)inv[i]=powM(fac[i]);
}
long long p[Maxn<<2];
//求出F(x-c) 
void fminus(long long *s,long long *f,int len,int c)
{
  c=mod-c;
  for (int i=0;i<len;i++)
    p[len-i-1]=f[i]*fac[i]%mod;
  long long buf=1;
  for (int i=0;i<len;i++,buf=buf*c%mod)
    s[i]=buf*inv[i]%mod;
  times(p,s,len,len);
  for (int i=0;i<len;i++)s[len-i-1]=p[i]*inv[len-i-1]%mod;
  for (int i=len;i<len+len;i++)s[i]=0;
}
long long f[Maxn<<2],s[Maxn<<2];
void solve(long long *f,int n)
{
  if (n==1){f[0]=0;f[1]=1;}
  else if (n&1){
    solve(f,n-1);f[n]=0;
    //再乘上(x-n+1)就好了
    for (int i=n;i>0;i--)
      f[i]=(f[i-1]+(mod-n+1)*f[i])%mod;
    f[0]=f[0]*(mod-n+1)%mod;
  }else {
    solve(f,n/2);
    //S(x)=F(x+n/2)
    fminus(s,f,n/2+1,n/2);
    times(f,s,n/2+1,n+1);
  }
}
void invp(long long *f,int len)
{
  for (int i=0;i<k+1;i++)s[i]=p[i]=0;
  //注意清空 
  long long *r=s,*rr=p;
  int n=1;for(;n<len;n<<=1);
  rr[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
  	for (int i=0;i<len;i++)
  	  r[i]=rr[i]*2%mod;
  	times(rr,rr,len/2,len);
  	times(rr,f,len,len);
  	for (int i=0;i<len;i++)
  	  rr[i]=(r[i]-rr[i]+mod)%mod;
  }for (int i=0;i<len;i++)
    f[i]=rr[i];
}
int main()
{
  scanf("%d%d",&n,&k);
  if (k>n){
    for (int i=0;i<=n;i++)printf("0 ");
    return 0;
  }invG=powM(G);
  Init(k);solve(f,k+1);
  for (int i=0;i<k+1;i++)f[i]=f[i+1];
  rev(f,k+1);
  for (int i=n-k+1;i<k+1;i++)f[i]=0;
  for (int i=k+1;i<n-k+1;i++)f[i]=0;
  invp(f,n-k+1);
  for (int i=0;i<k;i++)printf("0 ");
  for (int i=0;i<n-k+1;i++)printf("%lld ",f[i]);
  return 0;
}