话说原理八月初就会了，拖到现在才把代码写出来，是不是颓废之王?
前置知识：多项式乘法（FFT/NTT）（参考阅读：可能是废话最多的 FFT 教程）
一步一步推式子
我们设 $F(x)G(x)\equiv1\pmod{x^n} $，我们设这里的 \(n\) 是 \(2\) 的整数次幂。
然后设 \(H(x)F(x)\equiv1\pmod{x^{\frac{n}{2}}}\)，那么可以确定的是
\[G(x)F(x)\equiv1\pmod{x^{\frac{n}{2}}} \]应该不难理解，因为在模 \(x^n\) 时就没有非常数项余项（\(1\)）了，在模 \(x^{\frac{n}{2}}\) 时也不会有非常数项余项（\(1\)）的
于是我们有：
\[H(x)F(x)-G(x)F(x)\equiv F(x)\left(H(x)-G(x)\right)\equiv0\pmod{x^{\frac{n}{2}}} \]则：
\[H(x)-G(x)\equiv0\pmod{x^{\frac{n}{2}}} \]平方一下（注意模数也要平方）：
\[\left(H(x)-G(x)\right)^2\equiv H^2(x)-2H(x)G(x)+G^2(x)\equiv0\pmod{x^{\frac{n}{2}}} \]那你一定会纳闷了，这一通操作猛如虎，到底是要得到啥？
问的好，我们先考虑一个问题，当 \(n=2^0=1\) 时，设：
\[F(x)=a_0+a_1x+......+a_nx^n \]那么满足 \(F(x)G(x)\equiv1\pmod{x^n=x}\) 的 \(G(x)\) 只用是一个常值多项式即可，那么这个常值，其实就是 \(a_0\) 在模 \(998244353\) 意义下的逆元(因为一次项即以上在模 \(x^n=x\) 意义下不能造成实际贡献，只有常数项是会有贡献的)。
那么如果我们在 \(n=2^1=2\) 时，这个常值多项式就能作为 \(H(x)\) 出现，那么如果我们能得到从 \(H(x)\) 到 \(G(x)\) 的递推关系，就能得到 \(n=2\) 时的 \(G(x)\)，以此类推，就能求出最后的 \(G(x)\)。
这样的目的，也是指示我们选择平方将模数翻倍的动机所在。
那我们争取从前面的式子中解出 \(G(x)\)，但是存在一次方二次方，如果直接开根，那我们平方就没有意义了。
但别忘了在模 \(x^n\) 意义下还有一个关系：
\[G(x)F(x)\equiv1\pmod{x^n} \]为了出现 \(F(x)\) 还有相乘的格式，我们不妨在式子两边同时乘上 \(F(x)\)，即：
\[\begin{aligned}F(x)\left(H^2(x)-2H(x)G(x)+G^2(x)\right)&\equiv F(x)H^2(x)-2F(x)G(x)H(x)+F(x)G^2(x)\\&\equiv F(x)H^2(x)-2H(x)+G(x)\pmod{x^n}\end{aligned} \]于是我们就能解出 \(G(x)\) 关于 \(F(x)\) 与 \(H(x)\) 的递推关系式：
\[G(x)\equiv 2H(x)-F(x)H^2(x)\pmod{x^n} \]实现细节
由上面的结论我们可以得到模 \(x^{\frac{n}{2}}\) 意义下逆多项式到模 \(x^n\) 意义下逆多项式的转换，我们可以通过递推或递归方式实现。
而我是很不会递归的，所以就采用了递推的写法。
在实现时，因为数组不能当做函数的返回值，所以我们将多项式利用一个结构体 poly 存储，便于多项式运算的函数直接返回一个多项式。
另外二进制翻转数组 res[] 在多项式乘法结果最高次次数不同时也会变化，不能直接预处理一劳永逸，应该在每一次乘法之前更新一次。
利用结构体存储多项式易爆栈，可以在 [编译选项] 中设置连接器命令为“-Wl,--stack=128000000”以扩大栈空间，防止爆栈。
其他细节参见代码：
#include"iostream"
#include"cstring"
#include"cstdio"
#include"cmath"
using namespace std;
#define MAXN 100005
#define MOD 998244353
#define ll long long
int N,M,n=1,l=0;
struct poly
{
int a[MAXN*3];
poly(){memset(a,0,sizeof(a));}
}f;
int res[MAXN*3],w[MAXN*3];
int mul(int a,int b){return (ll)a*b%MOD;}//大数相乘取模
int quickpow(int a,int b)//快速幂
{
int ans=1,base=a;
while(b)
{
if(b&1) ans=mul(ans,base);
base=mul(base,base);
b>>=1;
}
return ans;
}
int inv(int a){return quickpow(a,MOD-2);}//数的逆元
void NTT(int *a,int unit,int n,int l)//没什么特别的 NTT
{
w[0]=1,w[1]=quickpow(unit,(MOD-1)/n);
for(int i=2;i<n;i++) w[i]=w[i]=mul(w[1],w[i-1]);
for(int i=0;i<n;i++) if(i>res[i]) swap(a[i],a[res[i]]);
for(int i=2;i<=l+1;i++)
{
int t=n>>(i-1);
for(int j=1;j<=t;j++)
{
int s=n/t;
for(int k=0;k<(s>>1);k++)
{
int op=s*(j-1)+k,G=mul(a[op+(s>>1)],w[k*t]);
a[op+(s>>1)]=(a[op]-G+MOD)%MOD;
a[op]=(a[op]+G)%MOD;
}
}
}
}
//多项式 F,G 的卷积，n 和 l 是长度和 n 对应的幂次
poly poly_mul(poly F,poly G,int n,int l)
{
poly h;
for(int i=0;i<n;i++) res[i]=(res[i>>1]>>1)|((i&1)<<(l-1));
//每次都要重新算一次 res 数组
NTT(F.a,3,n,l),NTT(G.a,3,n,l);
for(int i=0;i<n;i++) F.a[i]=mul(F.a[i],G.a[i]);
NTT(F.a,332748118,n,l);
int m=inv(n);
for(int i=0;i<n;i++) h.a[i]=mul(F.a[i],m);
return h;
}
poly poly_times(poly F,int k,int n)//数乘多项式
{
poly h;
for(int i=0;i<n;i++) h.a[i]=mul(F.a[i],k);
return h;
}
poly poly_min(poly F,poly G,int n)//长度为 n 的多项式 F 和 G 的差多项式
{
poly h;
for(int i=0;i<n;i++) h.a[i]=(F.a[i]-G.a[i]+MOD)%MOD;
return h;
}
poly poly_inv(poly F)//多项式求逆核心过程
{
poly h,g;
for(int i=0;i<n;i++) h.a[i]=0;
h.a[0]=inv(F.a[0]);//初始的 h(x) 就是 a_0 的逆元
for(int i=1;i<l;i++)
{
int s=1<<i;//求 F(x) 在模 x^s 次方意义下的逆元
g=poly_mul(h,h,s,i);
poly q;
for(int j=0;j<s;j++) q.a[j]=F.a[j];
//重新定义一个 q(x) 复制 F(x)，因为在 NTT 中的交换数组操作会破坏 F(x)
g=poly_mul(g,q,s<<1,i+1);
h=poly_min(poly_times(h,2,s>>1),g,s);
}
return h;
}
int main()
{
scanf("%d",&N);
while(n<=(N-1)*2) n<<=1,l++;//考虑一下为什么是这样的一个上界
for(int i=0;i<N;i++) scanf("%d",&f.a[i]);
poly H=poly_inv(f);
for(int i=0;i<N;i++) printf("%d ",H.a[i]);
return puts(""),0;
}
时间复杂度
多项式的时间复杂度主流算法即为主定理，然而我不会，又懒得学。
那我们就通过朴素的方式证明一下，因为多项式次数 \(n-1\) 和不超过 \(2(n-1)\) 的最小 \(2\) 的次幂同阶，于是设 \(n\) 就是一个 \(2\) 的幂次，且大于并最接近多项式次数的两倍，我们不妨讲问题简化为：
\[T(n)=\mathcal O(n\log n)\;\;\;\text{\{即为多项式乘法的复杂度\}} \]\[T(2^m)+T(2^{m-1})+......+T(1)=? \]其中 \(m=\log_2n\)。
先不计常数，即求：
\[M=2^m\log 2^m+2^{m-1}\log 2^{m-1}+......+2^1\log 2^1+2^0\log2^0 \]而：
\[\begin{aligned}M&= \log2\left(m2^m+(m-1)2^{m-1}+......+1\times2^1\right)\\&=\log2\left((m-1)2^{m+1}+2\right)\\&=\log2\left((\log_2n-1)\times2n+2\right)\\&=\mathcal O(n\log n)\end{aligned} \]这说明多项式求逆的时间复杂度为 \(\mathcal O(n\log n)\)。
多项式求导与积分
这一部分看起来最为吓人，但实现起来其实最为简单，我们设：
\[F(x)=a_0+a_1x+a_2x^2+......+a_nx^n \]\[P(x)=F'(x)=b_0+b_1x+b_2x^2+......+b_{n-1}x^{n-1} \]\[Q(x)=\int F(x)=c_0+c_1x+c_2x+......+c_{n+1}x^{n+1} \]注意次数加一减一关系，以及我们一般默认 \(c_0=0\)。
那么由导数公式 \((x^n)'=nx^{n-1}\) 以及导数积分互为逆运算的关系，我们能得到关系：
\[b_i=(i+1)a_{i+1}\;\;\;c_i=\frac{a_{i-1}}{i} \]容易知道，实现复杂度都是线性的。
多项式对数函数（多项式 ln）
借助导数与寄分（bushi）的工具，我们可以知道：
\[\ln F(x)\equiv \int \dfrac{F'(x)}{F(x)}\pmod {x^n} \]借助刚刚学会的多项式求逆就可以轻松解题了，实现上主要还是注意多项式求逆的部分，其余的部分实现难度相对少很多了，实现复杂度容易知道为 \(\mathcal O(n\log n)\)。
code:
#include"iostream"
#include"cstring"
#include"cstdio"
#include"cmath"
using namespace std;
#define MAXN 100005
#define MOD 998244353
#define ll long long
int N,M,n=1,l=0;
struct poly
{
int a[MAXN*3];
poly(){memset(a,0,sizeof(a));}
}f,g;
int res[MAXN*3],w[MAXN*3];
int mul(int a,int b){return (ll)a*b%MOD;}
int quickpow(int a,int b)
{
int ans=1,base=a;
while(b)
{
if(b&1) ans=mul(ans,base);
base=mul(base,base);
b>>=1;
}
return ans;
}
int inv(int a){return quickpow(a,MOD-2);}
void NTT(int *a,int unit,int n,int l)
{
w[0]=1,w[1]=quickpow(unit,(MOD-1)/n);
for(int i=2;i<n;i++) w[i]=w[i]=mul(w[1],w[i-1]);
for(int i=0;i<n;i++) if(i>res[i]) swap(a[i],a[res[i]]);
for(int i=2;i<=l+1;i++)
{
int t=n>>(i-1);
for(int j=1;j<=t;j++)
{
int s=n/t;
for(int k=0;k<(s>>1);k++)
{
int op=s*(j-1)+k,G=mul(a[op+(s>>1)],w[k*t]);
a[op+(s>>1)]=(a[op]-G+MOD)%MOD;
a[op]=(a[op]+G)%MOD;
}
}
}
}
poly poly_mul(poly F,poly G,int n,int l)
{
poly h;
for(int i=0;i<n;i++) res[i]=(res[i>>1]>>1)|((i&1)<<(l-1));
NTT(F.a,3,n,l),NTT(G.a,3,n,l);
for(int i=0;i<n;i++) F.a[i]=mul(F.a[i],G.a[i]);
NTT(F.a,332748118,n,l);
int m=inv(n);
for(int i=0;i<n;i++) h.a[i]=mul(F.a[i],m);
return h;
}
poly poly_times(poly F,int k,int n)
{
poly h;
for(int i=0;i<n;i++) h.a[i]=mul(F.a[i],k);
return h;
}
poly poly_min(poly F,poly G,int n)
{
poly h;
for(int i=0;i<n;i++) h.a[i]=(F.a[i]-G.a[i]+MOD)%MOD;
return h;
}
poly poly_der(poly F,int n)//求导
{
poly h;
for(int i=0;i<n;i++) h.a[i]=mul(i+1,F.a[i+1]);
return h;
}
poly poly_cal(poly F,int n)//积分
{
poly h;
for(int i=1;i<=n+1;i++) h.a[i]=mul(inv(i),F.a[i-1]);
return h;
}
poly poly_inv(poly F)
{
poly h,g;
for(int i=0;i<n;i++) h.a[i]=0;
h.a[0]=inv(F.a[0]);
for(int i=1;i<l;i++)
{
int s=1<<i;
g=poly_mul(h,h,s,i);
poly q;
for(int j=0;j<s;j++) q.a[j]=F.a[j];
g=poly_mul(g,q,s<<1,i+1);
h=poly_min(poly_times(h,2,s>>1),g,s);
}
return h;
}
int main()
{
scanf("%d",&N);
while(n<=(N-1)*2) n<<=1,l++;
for(int i=0;i<N;i++) scanf("%d",&f.a[i]);
g=poly_der(f,N-1);
poly H=poly_inv(f);
g=poly_mul(H,g,n,l);
for(int i=N+1;i<n;i++) g.a[i]=0;
g=poly_cal(g,N),g.a[0]=0;
for(int i=0;i<N;i++) printf("%d ",g.a[i]);
return puts(""),0;
}
呐，就给菜鸡博主一个赞吧，奖励一下自己看完了这么多。