这里都是一些论文级别的玩意，基本不是给正常人类看的
注意：这里仅对模数是 (998244353) 的部分进行介绍。
零.让我们开始
这里是一些基础的东西，不怎么需要想，这里就略过了。
Code
#include<bits/stdc++.h>
#define endl '
'
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
using namespace std;
const int P=998244353,G=3,LIMIT=50;
typedef vector<int> vec;
struct IO_Tp {
const static int _I_Buffer_Size = 2 << 22;
char _I_Buffer[_I_Buffer_Size], *_I_pos = _I_Buffer;
const static int _O_Buffer_Size = 2 << 22;
char _O_Buffer[_O_Buffer_Size], *_O_pos = _O_Buffer;
IO_Tp() { fread(_I_Buffer, 1, _I_Buffer_Size, stdin); }
~IO_Tp() { fwrite(_O_Buffer, 1, _O_pos - _O_Buffer, stdout); }
IO_Tp &operator>>(int &res) {
int f=1;
while (!isdigit(*_I_pos)&&(*_I_pos)!='-') ++_I_pos;
if(*_I_pos=='-')f=-1,++_I_pos;
res = *_I_pos++ - '0';
while (isdigit(*_I_pos)) res = res * 10 + (*_I_pos++ - '0');
res*=f;
return *this;
}
IO_Tp &operator<<(int n) {
if(n<0)*_O_pos++='-',n=-n;
static char _buf[10];
char *_pos = _buf;
do
*_pos++ = '0' + n % 10;
while (n /= 10);
while (_pos != _buf) *_O_pos++ = *--_pos;
return *this;
}
IO_Tp &operator<<(char ch) {
*_O_pos++ = ch;
return *this;
}
} IO;//快读
void chkmax(int &x,int y){if(x<y)x=y;}
void chkmin(int &x,int y){if(x>y)x=y;}
int qpow(int a,int k,int p=P){//快速幂
int ret=1;
while(k){
if(k&1)ret=1ll*ret*a%p;
a=1ll*a*a%p;
k>>=1;
}
return ret;
}
int norm(int x){return x>=P?x-P:x;}
int reduce(int x){return x<0?x+P:x;}
void add(int&x,int y){if((x+=y)>=P)x-=P;}//取模
struct Maths{
int n;
vec fac,invfac,inv;
void build(int n){
this->n=n;
fac.resize(n+1);
invfac.resize(n+1);
inv.resize(n+1);
fac[0]=1;
rep(k,1,n)fac[k]=1ll*fac[k-1]*k%P;
inv[1]=inv[0]=1;
rep(k,2,n)inv[k]=P-1ll*(P/k)*inv[P%k]%P;
invfac[0]=1;
rep(k,1,n)invfac[k]=1ll*invfac[k-1]*inv[k]%P;
}
Maths(){build(1);}
void chk(int k){
int lmt=n;
if(k>lmt){while(k>lmt)lmt<<=1;build(lmt);}
}
int cfac(int k){return chk(k),fac[k];}
int cifac(int k){return chk(k),invfac[k];}
int cinv(int k){return chk(k),inv[k];}
int binom(int n,int m){
if(m<0||m>n)return 0;
return 1ll*cfac(n)*cifac(m)%P*cifac(n-m)%P;
}
}math;//普通数论部分
struct poly{
vec a;
poly(int v=0):a(1){
if((v%=P)<0)v+=P;
a[0]=v;
}
poly(const vec&a):a(a){}
poly(initializer_list<int>init):a(init){}
int operator[](int k)const{return k<a.size()?a[k]:0;}
int&operator[](int k){
if(k>=a.size())a.resize(k+1);
return a[k];
}
int deg()const{return a.size()-1;}
void redeg(int d){a.resize(d+1);}
poly slice(int d)const{
if(d<a.size())return vec(a.begin(),a.begin()+d+1);
vec res(a);
res.resize(d+1);
return res;
}
int*base(){return a.data();}
const int*base()const{return a.data();}
poly println(FILE* fp)const{
fprintf(fp,"%d",a[0]);
rep(i,1,a.size()-1)fprintf(fp," %d",a[i]);
fputc('
',fp);
return *this;
}
poly operator+(const poly&rhs)const{
vec res(max(a.size(),rhs.a.size()));
rep(i,0,res.size()-1)if((res[i]=operator[](i)+rhs[i])>=P)res[i]-=P;
return res;
}
poly operator-()const{
poly ret(a);
rep(i,0,a.size()-1)if(ret[i])ret[i]=P-ret[i];
return ret;
}
poly operator-(const poly&rhs)const{return operator+(-rhs);}
/*
这里应该有一堆屎山声明，可是我懒得罗列了所以就没写
*/
poly shift(int k)const;
};//声明+部分简单函数
poly zeroes(int deg){return vec(deg+1);}//0函数
poly operator "" _z(unsigned long long a){return {0,(int)a};}
poly operator+(int v,const poly&rhs){return poly(v)+rhs;}//多项式加整数
poly operator*(int v,const poly&rhs){//多项式乘整数
poly ret=zeroes(rhs.deg());
rep(i,0,rhs.deg())ret[i]=1ll*rhs[i]*v%P;
return ret;
}
poly operator*(const poly&lhs,int v){return v*lhs;}
poly poly::shift(int k)const{//多项式乘 x^k
poly g=zeroes(deg()+k);
rep(i,0,k-1)g[i]=0;
rep(i,min(0,-k),deg()-1)g[i+k]=a[i];
return g;
}
template<class T>
IO_Tp& operator>>(IO_Tp& IO,vector<T>&v){//输入 vector
for(T&x:v)IO>>x;
return IO;
}
template<class T>
IO_Tp& operator<<(IO_Tp& IO,vector<T>&v){//输出 vector
for(T&x:v)IO<<x;
return IO;
}
一.多项式乘法
原理和普通的多项式乘法一致，没啥好说的，我们看看哪些地方可以优化。
我们注意到可以在初始化的时候做一些预处理，这样大概可以减少一定的常数，在 P3803 这道题上面总时间 (1.6s o 1s)，快了 0.6s。
补充：EI 认为 “经过测试，某些形式特殊的数组的 NTT 改良版本，看似省略了部分计算，实则缓存不友好，还不如直接做”。（虽然经过测试，某一版本的 NTT 比这一版块 0.2s，但是码量差不多翻了一倍（在 NTT 部分），因此不予使用。）
Code
struct NTT{
int L,brev[1<<11];
vec root;
NTT():L(-1){
rep(i,1,(1<<11)-1)brev[i]=brev[i>>1]>>1|((i&1)<<10);
}
void preproot(int l){
L=l;
root.resize(2<<L);
rep(i,0,L){
int *w=root.data()+(1<<i);
w[0]=1;
int omega=qpow(G,(P-1)>>i);
rep(j,1,(1<<i)-1)w[j]=1ll*w[j-1]*omega%P;
}
}
void dft(int*a,int lgn,int d=1){
if(L<lgn)preproot(lgn);
int n=1<<lgn;
rep(i,0,n-1){
int rev=(brev[i>>11]|(brev[i&((1<<11)-1)]<<11))>>((11<<1)-lgn);
if(i<rev)swap(a[i],a[rev]);
}
for(int i=1;i<n;i<<=1){
int *w=root.data()+(i<<1);
for(int j=0;j<n;j+=i<<1)rep(k,0,i-1){
int aa=1ll*w[k]*a[i+j+k]%P;
a[i+j+k]=norm(a[j+k]+P-aa);
add(a[j+k],aa);
}
}
if(d==-1){
reverse(a+1,a+n);
int inv=nt.inv(n);
rep(i,0,n-1)a[i]=1ll*a[i]*inv%P;
}
}
}ntt;
二.多项式乘法逆
首先列一下牛迭的式子：假设 (fin mathbb R[[x]],Ain mathbb R[[x,y]])，满足 (A(x,f)=0)，令 (f_0=fmod x^n)，则
如果定义 (operatorname{ord}(f)=min{n|[x^n]f eq 0})，那么我们可以观察到 (operatorname{ord}(A(f_0))ge n)，因此计算 (A'(f_0)) 的精度只需达到 (n) 即可。
下面讨论操作的优化。在以下讨论中，内容分为三部分：
- 直接按式子计算的时间。下面定义 (E(n)) 为一次长度为 (n) 的 DFT 所需的时间，(M(n)) 为一次两个精度为 (n) 的形式幂级数的乘法所需要的时间，因此 (M(n)=(3+o(1))E(2n)=(6+o(1))E(n))。在下文中，一切 (o(1)) 会被省略。
- 利用循环卷积优化。注意到在大部分情况下，我们已经得到了结果中的一部分系数，而长度为 (n) 的 DFT 解决了循环卷积问题 (fgmod (x^n-1))，仅用于计算卷积很浪费，所以可以考虑先计算循环卷积，必要时进行一些处理，最后得到所需的系数。
- 减少 DFT 次数。注意到很多时候计算了相同的 DFT，或者可以用线性变换的性质合并几次 IDFT，考虑减少这些额外的开销。
在这一部分，我们研究的问题是倒数。
1
(fin mathbb R[[x]])，令 (g=1/f)，求 (g)。
取 (A(g)=fg-1)，代入牛迭式子可得：
这就是我们一般使用的方法。
这里一共使用了 (1) 次长度 (2n) 的乘法，(1) 次长度 (4n) 的乘法，用时 (M(2n)+M(4n)=18E(n))。（其实我觉得做三次长度为 (2n) 的 DFT 就行了，这样是 (12E(n)) 的，不知道为什么没有这么写）
2
(fin mathbb R[[x]])，令 (g=1/f)，求 (g)。
考虑 (gmod x^{2n}=g_0-(fg_0-1)g_0mod x^{2n})，显然 (deg((fmod x^2n)g_0-1)<3n,operatorname{ord}((fmod x^2n)g_0-1)ge n)，因此只需计算 ((fmod x^{2n})g_0mod (x^{2n}-1)) 即可，同理 ((fg_0-1)g_0mod x^{2n}) 也只需要长为 (2n) 的循环卷积，用时 (12E(n))。所以计算 (gmod x^n) 的时间是 (12E(n))。
3
观察上述过程，有两次和 (g_0) 有关的长为 (2n) 的循环卷积，可以记录下来而不是重新算，用时 (10E(n))。
EI 的代码应该是按照这个实现的，不过把普通的递归换成了迭代，因此总时间 (0.6s o 0.2s)，优化了 (0.4s)。
Code
struct Newton{
void inv(const poly&f,const poly&nttf,poly&g,const poly&nttg,int t){//given f,g,nttf,nttg
int n=1<<t;
poly prod=nttf;
rep(i,0,(n<<1)-1)prod[i]=1ll*prod[i]*nttg[i]%P;
ntt.dft(prod.base(),t+1,-1);//calculate fg-1
rep(i,0,n-1)prod[i]=0;//prod=
ntt.dft(prod.base(),t+1,1);
rep(i,0,(n<<1)-1)prod[i]=1ll*prod[i]*nttg[i]%P;
ntt.dft(prod.base(),t+1,-1);//calculate (fg-1)g
rep(i,0,n-1)prod[i]=0;
g=g-prod;//calculate g-(fg-1)g
}
void inv(const poly&f,const poly&nttf,poly&g,int t){//given f,nttf,g
poly nttg=g;
nttg.redeg((2<<t)-1);
ntt.dft(nttg.base(),t+1,1);//calc nttg
inv(f,nttf,g,nttg,t);
}
void inv(const poly&f,poly&g,int t){//given f,g
poly nttg=g;
nttg.redeg((2<<t)-1);
ntt.dft(nttg.base(),t+1,1);//calc nttg
poly nttf=f;
nttf.redeg((2<<t)-1);
ntt.dft(nttf.base(),t+1,1);//calc nttf
inv(f,nttf,g,nttg,t);
}
}nit;
poly poly::inv()const{
poly g=nt.inv(a[0]);
for(int t=0;(1<<t)<=deg();++t)nit.inv(slice((2<<t)-1),g,t);
g.redeg(deg());
return g;
}
还可以继续优化吗？
如果允许长度为 (3n) 的 DFT，那么考虑 (gmod x^{2n}=g_0-(fg_0^2-g_0)mod x^{2n})，用长度为 (3n) 的循环卷积计算 (fg_0^2) 即可达到 (9E(n)) 的用时，可惜大部分时候（比如 (998244353)）都不能做。
注意到我们不是必须要循环卷积才能解决问题，对 (ain mathbb R,a^{2n} eq 1)，考虑在 (mathbb R[x]/(x^{2n}-1)(x^n-a^n)) 中计算卷积，即在 (1,zeta_{2n},zeta_{2n}^2,dots,zeta_{2n}^{2n-1},a,azeta_n,azeta_n^2,dots,azeta_n^{n-1}) 上多点求值和插值。对 (fin mathbb R[x]/(x^{2n}-1)(x^n-a^n)) 进行多点求值只需用 FFT 计算 (mathcal F_{2n}(f)) 和 (mathcal F_{2n}(f(ax)))，而插值只需分别还原 CRT 合并。
容易发现，如果在 (mathbb R[x]/(x^{2n}-1)(x^n-a^n)) 中进行卷积，仍然可以处理超出长度部分的影响，且不需要长度为 (3n) 的 DFT，同时也计算了 (mathcal F_{2n}(fmod x^{2n})) 或 (mathcal F_{2n}(g_0))，所需时间仍是改进前的 (9E(n)=frac{3}{2}M(n))，所以可以几乎完全代替前一种做法。
简单描述一下思路：为了算出所需结果 (f)，先算出 (fmod (x^{2n}−1)(x^n−a^n))，考虑超出部分对前 (n) 项（本应全是 (0)）的贡献，利用这些信息还原出这一部分，然后即可把这一部分对所需部分的影响消除掉。
算出结果需在 (1, zeta_{2n}, zeta_{2n}^2, dots, zeta_{2n}^{2n-1}, a, azeta_n, azeta_n^2, dots, azeta_n^{n-1}) 上多点求值和插值，多点求值即计算 (mathcal F_{2n}(f)) 和 (mathcal F_n(f(ax)))，插值即分别还原并 CRT 合并。
实际实现并不需要考虑这些，最后的推荐实现也可以这样描述。考虑将所需结果 (f) 表示为 (ax^n+bx^{2n}+cx^{3n})，其中 (a,b,cin mathbb R[x])，(deg(a),deg(b),deg(c)<n)，那么可以用循环卷积计算出 (fmod (x^{2n}-1)) 和 (fmod (x^n-i)=f(zeta_{4n}x)mod (x^n-1))，也就相当于算出了 (b,a+c,ia-b-ic)，还原出 (a) 即可。
三.多项式 ln
哈哈这个我熟，直接两边求导再积分，就可以得到 (ln f=int f'/f)，然后利用多项式求逆就可以做到 (9E(n)+6E(n)=15E(n)) 了！
一翻 EI 的代码：wtf 怎么这么长，怎么还有个 quo 函数？？？
EI：You are too naive.
翻了翻博客，发现求商数居然还有科技，学废了。
1
商数：对于 (f,hin mathbb R[[x]])，令 (g=1/f,q=hg=h/f)，求 (q)。
显然先求 (g=1/f)，再求 (q=hg) 即可，总时间是 (18E(n)+6E(n)=24E(n))。
对数：略。
2
对于 (f,hin mathbb R[[x]])，令 (g=1/f,q=hg=h/f)，求 (q)。
直接求 (g) 然后卷积的总时间是 (18E(n))。
如果不求 (g)，注意到 (A(q)=fq-h=0)，令 (g_0=gmod x^n,h_0=hmod x^n,q_0=qmod x^n=g_0h_0mod x^n)，有
计算 (g_0) 是 (12E(n))，计算 (q_0) 是 (6E(n))，计算 ((fq_0-h)) 和倒数类似需要 (12E(n)) 的时间，因此计算 (qmod x^{2n}) 是 (30E(n))，计算 (qmod x^n) 就是 (15E(n))。
3
对于 (f,hin mathbb R[[x]])，令 (g=1/f,q=hg=h/f)，求 (q)。
令 (g_0=gmod x^n,g_1=(gmod x^{2n}-g_0)/x^n,h_0=hmod x^n,h_1=(hmod x^{2n}-h_0)/x^n)。
如果需要求 (g)，考虑计算
- 计算 (mathcal F_{2n}(g_0),g_0,g_1) 需要 (18E(n)) 时间。
- 计算 (mathcal F_{2n}(g_1),mathcal F_{2n}(h_0),mathcal F_{2n}(h_1)) 需要 (6E(n)) 时间。
- 计算 (g_0h_0,g_0h_1+g_1h_0) 需要 (4E(n)) 时间。
总共需要 (28E(n)) 时间计算 (qmod x^{2n})，因此计算 (qmod x^n) 就需要 (14E(n)) 时间。
这里用到的技巧可以表述为：对于 (f,gin mathbb R[[x]])，已知 (fmod x^n,gmod x^n,mathcal F_n(fmod x^{n/2}))，则需要 (5E(n)) 时间计算出 (fgmod x^n)。
如果不求 (g)，仍然考虑
与第二部分的做法相比，可以用更快计算倒数的方法，计算 (q_0) 时用的 (mathcal F_{2n}(g_0)) 可以用于计算 ((fq_0-1)g_0)，所以计算 (qmod x^{2n}) 所用时间为 (24E(n))，计算 (qmod x^n) 总时间为 (12E(n)=2M(n))。
- 计算 (g_0) ( o) (9E(n))
- 计算 (mathcal F_{2n}(g_0),q_0) ( o) (6E(n))
- 计算 (fq_0) ( o) (M(n)=6E(n))
- 计算 ((fq_0-h)g_0)，已知 (mathcal F_{2n}(g_0)) ( o) (4E(n))。
因此计算 (qmod x^{2n}) 总时间是 (25E(n))，计算 (qmod x^n) 总时间是 (12.5E(n))。
观察 (q=g_0h_0mod x^n,(fq_0-h)g_0)，可以发现符合上文描述的技巧的使用条件，其中 ((fq_0-h)g_0) 可以视为计算 (((fq_0-h)/x^n)g_0)，再考虑到相同的 DFT 只用计算一次，需要 (9E(n)) 时间计算，总时间就是 (12E(n)=2M(n))。
Code
poly poly::quo(const poly&rhs)const{
if(rhs.deg()==0)return 1ll*a[0]*nt.inv(rhs[0])%P;
poly g=nt.inv(rhs[0]);
int t=0,n;
for(n=1;(n<<1)<=rhs.deg();++t,n<<=1)nit.inv(rhs.slice((n<<1)-1),g,t);
poly nttg=g;
nttg.redeg((n<<1)-1);
ntt.dft(nttg.base(),t+1,1);
poly eps1=rhs.slice((n<<1)-1);
ntt.dft(eps1.base(),t+1,1);
rep(i,0,(n<<1)-1)eps1[i]=1ll*eps1[i]*nttg[i]%P;
ntt.dft(eps1.base(),t+1,-1);
memcpy(eps1.base(),eps1.base()+n,sizeof(int)<<t);
memset(eps1.base()+n,0,sizeof(int)<<t);
ntt.dft(eps1.base(),t+1,1);
poly h0=slice(n-1);
h0.redeg((n<<1)-1);
ntt.dft(h0.base(),t+1,1);
poly h0g0=zeroes((n<<1)-1);
rep(i,0,(n<<1)-1)h0g0[i]=1ll*h0[i]*nttg[i]%P;
ntt.dft(h0g0.base(),t+1,-1);
poly h0eps1=zeroes((n<<1)-1);
rep(i,0,(n<<1)-1)h0eps1[i]=1ll*h0[i]*eps1[i]%P;
ntt.dft(h0eps1.base(),t+1,-1);
rep(i,0,n-1)h0eps1[i]=reduce(operator[](i+n)-h0eps1[i]);
memset(h0eps1.base()+n,0,sizeof(int)<<t);
ntt.dft(h0eps1.base(),t+1,1);
rep(i,0,(n<<1)-1)h0eps1[i]=1ll*h0eps1[i]*nttg[i]%P;
ntt.dft(h0eps1.base(),t+1,-1);
memcpy(h0eps1.base()+n,h0eps1.base(),sizeof(int)<<t);
memset(h0eps1.base(),0,sizeof(int)<<t);
return (h0g0+h0eps1).slice(rhs.deg());
}
</details>