LibreOJ #6268. 分拆数 生成函数+多项式exp

题意

求1到n的分拆数。
$n=10^5$

分析

设分拆数的生成函数为$$F(x)=\sum _{i=1}^{n}f(i)x^i$$
显然有
$$F(x)=\prod _{k=1}^n\frac{1}{1-x^k}$$
两边取对数得
$$\ln F(x)=\sum _{i=1}^n\ln \frac{1}{1-x^k}$$
将右边的式子泰勒展开一下可得
$$\ln F(x)=\sum _{k=1}^n\sum _{i\ge 1}\frac{x^{ki}}{i}$$
先$O(nlogn)$预处理出右边然后再exp回去即可。

如何求$e^{A(x)}$：
设$$\ln F(x)\equiv A(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
$$\ln H(x)\equiv A(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}})$$
若已经求出了$H(x)$，由牛顿迭代可得
$$F(x)=H(x)-\frac{\ln H(x)-A(x)}{\frac{1}{H(x)}}=H(x)(1-\ln H(x)+A(x))$$

如何求$\ln A(x)$：
设$$B(x)=\ln A(x)$$
两边求导得
$$B^{\prime }(x)=\frac{A^{\prime }(x)}{A(x)}$$
先把$A(x)$分别求导，求逆后相乘，然后积分回去就好了。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>

typedef long long LL;

const int N=400005;
const int MOD=998244353;

int n,a[N],rev[N],G[N],L,lg,b[N],H[N],F[N],ny[N];

int ksm(int x,int y)
{
	int ans=1;
	while (y)
	{
		if (y&1) ans=(LL)ans*x%MOD;
		x=(LL)x*x%MOD;y>>=1;
	}
	return ans;
}

void NTT(int *a,int f)
{
	for (int i=0;i<L;i++) if (i<rev[i]) std::swap(a[i],a[rev[i]]);
	for (int i=1;i<L;i<<=1)
	{
		int wn=ksm(3,f==1?(MOD-1)/i/2:MOD-1-(MOD-1)/i/2);
		for (int j=0;j<L;j+=(i<<1))
		{
			int w=1;
			for (int k=0;k<i;k++)
			{
				int u=a[j+k],v=(LL)w*a[j+k+i]%MOD;
				a[j+k]=(u+v)%MOD;a[j+k+i]=(u+MOD-v)%MOD;
				w=(LL)w*wn%MOD;
			}
		}
	}
	if (f==-1)
	{
		int ny=ksm(L,MOD-2);
		for (int i=0;i<L;i++) a[i]=(LL)a[i]*ny%MOD;
	}
}

void get_inv(int *a,int n)
{
	if (n==1) {G[0]=ksm(a[0],MOD-2);return;}
	get_inv(a,n/2);
	for (L=1,lg=0;L<=n;L<<=1,lg++);
	for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1)),b[i]=0;
	for (int i=0;i<n;i++) b[i]=a[i];
	NTT(G,1);NTT(b,1);
	for (int i=0;i<L;i++) G[i]=(G[i]*2%MOD+MOD-(LL)G[i]*G[i]%MOD*b[i]%MOD)%MOD;
	NTT(G,-1);
	for (int i=n;i<L;i++) G[i]=0;
}

void get_ln(int *a,int n)
{
	for (int i=0;i<=n*2;i++) G[i]=0;
	get_inv(a,n);
	for (int i=0;i<n;i++) H[i]=(LL)a[i+1]*(i+1)%MOD;
	for (L=1,lg=0;L<=n;L<<=1,lg++);
	for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	NTT(H,1);NTT(G,1);
	for (int i=0;i<L;i++) H[i]=(LL)H[i]*G[i]%MOD;
	NTT(H,-1);
	for (int i=n-1;i>=1;i--) H[i]=(LL)H[i-1]*ksm(i,MOD-2)%MOD;
	H[0]=0;
	for (int i=n;i<L;i++) H[i]=0;
}

void get_exp(int *a,int n)
{
	if (n==1) {F[0]=1;return;}
	get_exp(a,n/2);
	for (int i=0;i<=n*2;i++) H[i]=0;
	get_ln(F,n);
	for (L=1,lg=0;L<=n;L<<=1,lg++);
	for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1)),b[i]=0;
	b[0]=1;
	for (int i=1;i<n;i++) b[i]=(a[i]+MOD-H[i])%MOD;
	NTT(F,1);NTT(b,1);
	for (int i=0;i<L;i++) F[i]=(LL)F[i]*b[i]%MOD;
	NTT(F,-1);
	for (int i=n;i<L;i++) F[i]=0;
}

int main()
{
	scanf("%d",&n);
	ny[0]=ny[1]=1;
	for (int i=2;i<=n;i++) ny[i]=(LL)(MOD-MOD/i)*ny[MOD%i]%MOD;
	for (int k=1;k<=n;k++)
		for (int i=1;i*k<=n;i++)
			(a[k*i]+=ny[i])%=MOD; 
	int k=1;
	for (;k<=n;k<<=1);
	get_exp(a,k);
	for (int i=1;i<=n;i++) printf("%d\n",F[i]);
	return 0;
}