实现卷积的算法

First.FFT

两个多项式做乘法，暴力的想法是系数依次考虑，复杂度O(n^2)

然而如果我们观察出来了一个性质，不妨设我们现在想求解：

F(x)*G(x)

那么显然有乘出来的多项式P(x)中第k项的系数c_k满足：

c_k=\sum_{i+j=k}a_i*b_j

这是一个加法卷积。

所以FFT说他用来做加法卷积也可以，说他用来搞多项式乘法也可以

接下来是FFT的核心思想：

将两个多项式转成点值表达。

将点值两两乘起来。

将点值换成系数表达。

前面那一部分我们叫做DFT，后面那一部分叫做IDFT（逆）

前方高能：

当我们在做求点值表达的时候，会发现复杂度很高，原因是因为我们没有代入合适的值。

所以我们考虑代入合适的值进去。

单位根\omega_{n}^k表示单位根

他是定义在复平面的单位圆上的。

n$次单位根会将此圆$n$等分，然后取k块作为$\omega_{n}^k

显然\omega_n^{0}=1

以及n次单位根的n次方=1

复数乘法有个好玩的，就是它的乘法遵循模长相乘幅角相加。

因为单位根定义在单位圆上所以模长是1，而他们的幅角就加起来好了。

我们尝试给一个n-1次多项式代入n个单位根进去（以下n均视为2的整数次幂）

会由于下列两条性质而发生翻天覆地的变化：

第一： $$\omega_{n}^{2k} = \omega_{n/2}^{k}$$ 这个还好吧，显然有$\omega_{n}^k=(\omega_n^1)^k

且有：\omega_{n/2}^1=\omega_{n}^2

第二：

\omega_{n}^{k+n/2}=-\omega_{n}^k

原因很简单

这里我们定义的负数的关于原点对称了的。

且因为：\omega_{n}^{n/2}=-1

所以\omega_{n}^{n/2}*\omega_{n}^k=-\omega_{n}^k

$$\sum_{j=0}^{n-1}(\omega_{n}^k)^j=0/n$$ 简单证明一下： 如果$k=0$，那么显然答案为$n

若k\ne0，那么我们可以记s_n=\sum_{j=0}^{n-1}(\omega_{n}^k)^j

于是有：

s_n*\omega_{n}^k=s_n - 1+(\omega_{n}^k)^{n-1+1}

稍微化简一下得：s_n=\dfrac{(\omega_{n}^k)^n-1}{\omega_{n}^k-1}

显然有分母不为0，而分子=0

于是我们得到了求和性质。

前方真高能！

我们有了折半性质后怎么做DFT?

对于一个多项式F(x)

不妨设它为F(x)=a_0+a_1x+a_2x^2+...a_{n-1}x^{n-1}

不妨设：

Fl(x)=a_0+a_2x+a_4x^2+...a_{n-2}x^{\frac{n}{2}-1} Fr(x)=a_1+a_3x+a_5x^2+...a_{n-1}x^{\frac{n}{2}-1}

于是有：

F(x)=Fl(x^2)+xFr(x^2)

接下来我们代入\omega_n^k

得到：

F(\omega_n^k)=Fl((\omega_n^k)^2)+\omega_n^k*Fr((\omega_n^k)^2)

于是：

F(\omega_n^k)=Fl(\omega_{n/2}^k)+\omega_n^k*Fr(\omega_{n/2}^k)

这里假设k< n/2

所以我们想要求出k取[0,n/2)的所有值的单位根的点值都只要递归处理Fl,Fr即可。

接下来我们考虑k\ge n/2

不妨设它为k+n/2(k<n/2)

于是有：

F(\omega_n^{k+n/2})=Fl((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}*Fr((\omega_n^{k+n/2})^2)

得到：

F(\omega_n^{k+n/2})=Fl((-\omega_n^{k})^2)-\omega_n^{k}*Fr((-\omega_n^{k})^2)

即：

F(\omega_n^k)=Fl(\omega_{n/2}^k)-\omega_n^k*Fr(\omega_{n/2}^k)

你会惊人的发现这个东东和前面的一样只是一个符号不同而已

这同时表明我们只要知道Fl,Fr的n/2个点值就可以求出F的n个点值。

于是复杂度就是T(x)=x+2*T(\dfrac{x}{2})=x\log x

即n\log n

接下来我们考虑插值，即知道n个点值会逆推得到系数。

这里设我们得到的点值为A

我们发现实际上我们要做的就是求解这样的一组方程组：

(\omega_{n}^0)^0*a_0+(\omega_{n}^0)^1*a_1+...+(\omega_{n}^0)^{n-1}*a_{n-1}=A(\omega_{n}^0) (\omega_{n}^1)^0*a_0+(\omega_{n}^1)^1*a_1+...+(\omega_{n}^1)^{n-1}*a_{n-1}=A(\omega_{n}^1) ...... (\omega_{n}^{n-1})^0*a_0+(\omega_{n}^{n-1})^1*a_1+...+(\omega_{n}^{n-1})^{n-1}*a_{n-1}=A(\omega_{n}^{n-1})

我们可以将之写作矩阵的形式。

于是我们就只要对左边的矩阵求逆了。

假设现在的式子是：V*X=F

那么我们需要求出V'使得V'*V*X=X=F*V'

这个矩阵可以构造出来，我们利用一下求和性质。

于是新矩阵的第i,j项为：

c_{i,j}=\sum_{k=0}^{n-1}a_{i,k}*b_{k,j}

如果我们令V'中的每一项都恰好为V的相反数（雾）

即(i,j)为(\omega_{n}^{-i})^j

于是可以得到：

c_{i,j}=\sum_{k=0}^{n-1}(\omega_{n}^{-i})^k*(\omega_{n}^k)^j

也就是：

c_{i,j}=\sum_{k=0}^{n-1}(\omega_{n}^{-i})^k*(\omega_{n}^j)^k c_{i,j}=\sum_{k=0}^{n-1}(\omega_{n}^{j-i})^k

由求和性质

当i=j时，其为n

否则其为0

所以我们构造出来的V'还要除一个n才是单位矩阵

也就有：V'*F/n=X

我们考虑求解V'*F

因为F固定，所以我们尝试将之看作一个新的系数

首先V'*F肯定是一个大小为1*n的矩阵

所以有：V'*F后得到的第i项其实就是：

A_0*(\omega_{n}^{-i})^0+A_1*(\omega_{n}^{-i})^1+...+A_{n-1}*(\omega_{n}^{-i})^{n-1}

你惊人的发现它变成了我们想要求出当x取\omega_{n}^{-0},\omega_{n}^{-1}...\omega_{n}^{-(n-1)}时的点值了

于是就再做一遍FFT，不过是每次取\omega_{n}^{-1}为基本增量。

复杂度O(n\log n)

Second.NTT

我们尝试利用整数int类型来代替浮点数进行计算。

首先是阶的定义：（下面的等于(=)是同余）

如果有a^n=1(mod~p)

则称满足此条件最小的n为a模p的阶

若一个数g的阶为\phi(p)，则称g为p的原根

我们尝试使用原根代替单位根进行计算。

首先要有单位根的如下几条性质：

其次是折半性质$(DFT)$和求和性质$(IDFT)

我们尝试证明原根亦满足如下三条性质。

1.$若$g$为$p$的原根，则$g^i\%p$互不相等$(0\le i <\phi(p))

证明：

若存在两者相等，不妨设k>j且g^k=g^j(mod~p)

于是有：g^{k-j}=1(mod~p)

又k<\phi(p),j<\phi(p)

故k-j<\phi(p)

与定义矛盾，得证。

2.$我们定义原根$g_{n}^1$为$g^{\frac{p-1}{n}}

那么g_{n}^k=(g_n^1)^k，由性质1，可以得到原根互不相等。

不难得到：g_{n}^{2k}=g_{n/2}^k(mod~p)

把两边拆开，代入g_{n}^1即可。

其次：

我们接下来考虑将之指数看作一个角度为p-1的大圆，那么里面的n就是将之n等分，其余与单位根相似。

那么其乘法就是幅角不断累加。

于是可以得到：-g_n^k=p-g_n^k=g_n^{k+n/2}

这里就感性理解吧，我实在是不会证了，按照单位根的方法理解。

性质3:求和性质

求和性质利用了等比数列求和的定理，我们套用即可，因为是在\%p意义下，乘法操作仍然满足。

故g满足求和性质。

有了这三个性质，我们就能愉快的把原根代入单位根运算了。

然后因为FFT中的DFT和IDFT的n均为2^n

所以我们取的p要p-1整除足够大的n

比较常用的p是998244353，它的原根为3

且998244353=119*8388608(2^{23})+1

Third.非递归版FFT和NTT

我们尝试优化它，从递归变成循环 仔细划分$FFT$会发现其呈现出树的结构 我们假设我们通过某种手段快速的得到了其最后的序列（叶子节点） 那么就可以两两往上合并来求答案 比如原序列： $0~ 1~ 2~ 3~ 4~ 5~ 6~ 7

最后得到的序列应该是:

0~4~2~6~2~5~3~7

我们打个表，发现：

000,001,010,011,100,101,110,111

变成了：

000,100,010,110,001,101,011,111

其实就是二进制反转了

考虑求出这个东东

貌似有一个递推式，不妨记R[i]为i反转后的值

于是有：R[i]=(R[i>>1]>>1)|((i\&1)<<(L-1))

简单讲讲 我们假设$[0,x)$的反转值都已经求出，现在想求出$x$的反转值 如果$x$为奇数，那么$x$可以写作$2*(x>>1)|1

那么有我们将x>>1反转后，可以将之右移一位，然后x还要在其最高位补上一个1

偶数类似，不用补1而已

这样就O(n)

然后因为这个关系双向，所以我们做swap即可。

#include<bits/stdc++.h>
using namespace std;
#define rep( i, s, t ) for( register int i = s; i <= t; ++ i )
#define re register
#define LL long long
int read() {
    char cc = getchar(); int cn = 0, flus = 1;
    while(cc < '0' || cc > '9') {  if( cc == '-' ) flus = -flus;  cc = getchar();  }
    while(cc >= '0' && cc <= '9')  cn = cn * 10 + cc - '0', cc = getchar();
    return cn * flus;
}
const int P = 998244353 ;
const int G = 3 ; 
const int Gi = 332748118 ;
const int N = 4e6 + 5 ; 
int n, m, R[N], L, limit ;
LL A[N], B[N] ;
LL fpow( LL x, int k ) {
    LL ans = 1, base = x ;
    while( k ) {
        if( k & 1 ) ans = ans * base % P ;
        k >>= 1, base *= base, base %= P ; 
    }
    return ans ; 
}
void NTT( LL *a, int type ) {
    for( int i = 0; i < limit; ++ i )
    if( i < R[i] ) swap( a[i], a[R[i]] ) ;
    for( int k = 1; k < limit; k <<= 1 ) {
        LL dg = fpow( ( type == 1 ) ? G : Gi, ( P - 1 ) / ( k << 1 ) ); //这里要注意
        for( int i = 0; i < limit; i += ( k << 1 ) ) //注意*2
        for( LL j = i, g = 1; j < i + k; ++ j, g = ( g * dg ) % P ) {
            LL Nx = a[j], Ny = ( a[j + k] * g ) % P;
            a[j] = ( Nx + Ny ) % P, a[j + k] = ( Nx - Ny + P ) % P ; 
        }
    }
}
void init() {
    while( limit <= n + m ) limit <<= 1, ++ L ; 
    for( int i = 0; i < limit; ++ i ) R[i] = ( R[i >> 1] >> 1 ) | ( ( i & 1 ) << ( L - 1 ) ) ;
}
signed main()
{
    n = read(), m = read(), limit = 1 ; 
    rep( i, 0, n ) A[i] = read() ; 
    rep( i, 0, m ) B[i] = read() ; 
    init(), NTT( A, 1 ), NTT( B, 1 ) ;
    rep( i, 0, limit ) A[i] = ( A[i] * B[i] ) % P ; 
    NTT( A, -1 ) ; int inv = fpow( limit, P - 2 ) ;
    rep( i, 0, n + m ) printf("%lld ", ( A[i] * inv ) % P ) ; 
    return 0;
}

Fourth.MTT/任意模数FFT

貌似有很多做法，先记着一个我会的，最慢的，7FFT

将每一项的系数拆成k1*m+b1=a,k2*m+b2=b

于是显然有：m*K1(x)+B1(x)=A(x)，另一个同理。

我们就分别对K1,B1,K2,B2做一遍DFT

然后考虑合并

其在x处的取值，我们做乘法本来是A(x)*B(x)

变成(m*K1+B1)*(m*K2+B2)

然后暴力展开，发现还要对K1*K2,K1*B2+K2*B1,B1*B2做一遍IDFT

最后合并系数即可，7FFT

以后学了别的会慢慢补充.......

MTT

考虑对于两个多项式A,B

构造P=A+iB,Q=A-iB

显然有，记P_k表示P在\omega_{n}^k处的取值

记录X=\dfrac{2\pi jk}{n}

P_k=A(\omega_{n}^k)+i*B(\omega_n^k) P_k=\sum_{j=0}^kA_j(\omega_n^{jk})+i*B_j(\omega_n^{jk}) P_k=\sum_{j=0}^k(A_j+i*B_j)*\omega_{n}^{jk} P_k=\sum_{j=0}^k(A_j+i*B_j)*(cosX+isinX) Q_k=\sum_{j=0}^kA_j(\omega_n^{jk})-i*B_j(\omega_n^{jk}) Q_k=\sum_{j=0}^k(A_j-i*B_j)*\omega_{n}^{jk} Q_k=\sum_{j=0}^k(A_j-i*B_j)(cosX+isinX) Q_k=\sum_{j=0}^k(A_jcosX+B_jsinX)+i*(A_jsinX-B_jcosX) Q_k=conj(\sum_{j=0}^k(A_jcosX+B_jsinX)-i*(A_jsinX-B_jcosX)) Q_k=conj(\sum_{j=0}^kA_j(cosX-i*sinX)+i*B_j(cosX-i*sinX)) Q_k=conj(\sum_{j=0}^k(A_j+i*B_j)(cosX-i*sinX)) Q_k=conj(\sum_{j=0}^k(A_j+i*B_j)\omega_{n}^{-jk}) Q_k=conj(\sum_{j=0}^k(A_j+i*B_j)(\omega_{n}^{n-k})^j) Q_k=conj(P_{n-k})

所以我们只要DFT出P，就能还原Q了

然后求出了P和Q后

就有：

A=\dfrac{P+Q}{2},B=\dfrac{P-Q}{2}*(-i)

然后对A*B做一遍IDFT即可

3$次$FFT->2FFT

上述优化过程后得到的\rm FFT叫做\rm MTT

接下来考虑优化\rm IDFT

假设我们有A的点值和B的点值

根据求和性质，我们只需要求出其在\omega_{n}^{-k}的点值就可以get到系数

不妨构造P=\dfrac{A+iB}{2},Q=\dfrac{A-iB}{2}

P_k=\sum_{j=0}^k(A_j+i*B_j)*(\omega_{n}^{n-k})^j Q_k=\sum_{j=0}^kA_j(\omega_n^{n-k})^j-i*B_j(\omega_n^{n-k})^j

令X=\dfrac{2\pi*(-kj)}{n}

同上可得：

Q_k=conj(P_{n-k})

貌似还有MTT2，可以做到1.5FFT搞出结果，不会告辞

这样大概可以将两次算点值做到一次同时get

然后就可以将7FFT这种丢人的东西变成4FFT...

我寻思大家为什么都用鱼写的那个奇怪的东西叫做三次变两次啊...真正的三次变两次不应该是这玩意儿吗...这是真正意义上同时DFT出两个点值啊...

然后就要开始上全家桶了qwq

华丽的分割线

多项式全家桶

First.多项式求逆

给定n-1次多项式F(x)，求解G(x)使得：

G(x)*F(x)=1(mod~x^n)

首先显然那些什么n次幂及以上的部分对答案没有影响

于是我们考虑递推

不妨设已知：H(x)*F(x)=1(mod~x^{n/2})

则显然：G(x)*F(x)=1(mod~x^{n/2})

故做差，得到：

G(x)-H(x)=0(mod~x^{n/2})

平方得到：

(G(x)-H(x))^2=0(mod~x^n)

于是有：

G^2+H^2-2G*H=0(mod~x^n)

两边同时乘F得到：

G+F*H^2-2H=0(mod ~x^n)

于是有：

G=2H-F*H^2

如此递推即可。

从1递推到2^k>n比较合适

复杂度O(n\log n)

关于证明：

T(x)=x\log x+T(x/2)

于是T'(x)=\log x*(O'(x))>T(x)

O'(x)=x+O(x/2)=O(n)

所以T'(x)=\log x *O'(x)=n\log n

Second. 【模板】分治 FFT

考虑构造一个函数F(x)=\sum_{i=0}^{\infty}f_i*x^{i}

当然这样的f我们无限的生成下去

类似的构造一个G(x)，没输入的地方都取0

于是显然有：F(x)*G(x)+f_0=F(x)

然后就有：F(x)=f_0*(1-G(x))^{-1}

因为我们只要求出n-1项，所以有：

F(x)=f_0*(1-G(x))^{-1}(mod~x^n)

因为f_0=1，所以我们直接跑多项式求逆就好了...

Third.【模板】多项式除法

给定F,G

求出Q,R使得：

F=G*Q+R

注意到如下性质：

f(\dfrac{1}{x})*x^n=f_r(x) 于是有： $$F(x)=G*Q(x)+R(x)$$ $$F(\dfrac{1}{x})*x^n=G*Q(\dfrac{1}{x})*x^n+R(\dfrac{1}{x})*x^n$$ $$F(\dfrac{1}{x})*x^n=G(\dfrac{1}{x})*x^m*Q(\dfrac{1}{x})*x^{n-m}+R(\frac{1}{x})*x^{m-1}*x^{n-m+1}$$ $$F_r(x)=G_r(x)*Q_r(x)+R_r(x)*x^{n-m+1}$$ 同时取模得到： $$F_r(x)=G_r(x)*Q_r(x)(~mod~ x^{n-m+1}~)$$ $$F_r(x)*G_r^{-1}(x)=Q_r(x)(~mod~ x^{n-m+1}~)$$ 于是就求个逆，然后做乘法$......

Sixth.多项式\ln/多项式\exp(4\&5为泰勒展开+牛顿迭代)

这个东西要前置芝士...

1.求导

定义式：

f'(x)=\lim_{dx\to 0}\dfrac{f(x+dx)-f(x)}{dx}

关于求导运算的一些结论：

1.(f(x)*g(x))'=f'(x)*g(x)+f(x)*g'(x)

证明：

(f*g)'(x)=\dfrac{(f*g)(x+dx)-(f*g)(x)}{dx} \dfrac{(f*g)(x+dx)-f(x+dx)g(x)+f(x+dx)g(x)-(f*g)(x)}{dx} \dfrac{f(x+dx)(g(x+dx)-g(x))+g(x)*(f(x+dx)-f(x))}{dx}

于是就有：

因为\lim_{dx\to 0}f(x+dx)=f(x)

所以就有：

原式=

f(x)*g'(x)+g(x)*f'(x)

于是：

(f*g)(x)=f(x)*g'(x)+g(x)*f'(x) 2.f(g(x))'=f'(g(x))*g'(x)

证明：

f(g(x))'=\dfrac{f(g(x+dx))-f(g(x))}{dx}

可以构造dg=\dfrac{}{dx}

咕咕咕...

(\dfrac{f(x)}{g(x)})'=\dfrac{f'(x)g(x)-g'(x)f(x)}{g(x)^2} (f/g)'=f\times (g^{-1})'=f\times (g^{-1})'+f'/g

注意到 f^k 是一个复合导，所以 (f^k)'=f'\times k\times f^{k-1}

所以：

(f/g)'=(-1)fg'/g^2+f'/g=(f'g-fg')/g^2

一些比较常用的结论：

1.\ln'(x)=\dfrac{1}{x} 2.e'^x=e^x 3.f(x)=x^k,f'(x)=k*x^{k-1}

3.式比较有用

2.积分

积分和求导为逆运算，定义式：

\int dx*f(x)

然后积分可加...(面积可加)

于是同理求导可加

然后积分可以直接乘以倍数/系数，反过来求导也可以qwq

Fourth.泰勒展开

首先是一个函数的泰勒展开：

F(x)=f(x_0)+\sum_{i=1}^{\infty}\dfrac{f^{i}(x_0)(x-x_0)^i}{i!}

其中f^i表示i阶导数

然而在精度允许的条件下由于阶乘非常的大所以只需要把n开到12,13左右即可

然后由于我们可以随便选择一个x_0进行泰勒展开，所以一般选择0

于是我们就可以将原式转化为

F(x)=f(0)+\sum_{i=1}^{\infty}\dfrac{f^i(0)*x^i}{i!}

Firth.牛顿迭代

用于求解0点

与二分不同的是牛顿迭代的求法是对一个函数求导，降次，求解降次后的0点以逼近原函数的0点

多项式牛顿迭代即求解满足G(F(x))=0(mod~~x^n)的解F(x)

n=1$的时候满足条件的多项式可以直接求出来，这个时候只有常数项，令常数项$=0

假设现在已经求出了：

G(F_0(x))=0(mod ~~x^{n/2})

考虑拓展到x^n次意义下，可以在G(F_0(x))处将G(F(x))泰勒展开变成：

G(F(x))=G(F_0(x))+\sum_{i=1}^n\dfrac{G^{i}(F_0(x))*(F(x)-F_0(x))^i}{i!}

因为我们有：G(F(x))=0(mod~~x^n)

所以就有：

G(F(x))=G(F_0(x))+\sum_{i=1}^n\dfrac{G^{i}(F_0(x))*(F(x)-F_0(x))^i}{i!}=0(mod~~x^n)

又G(F_0(x))=0(mod~~x^{n/2})

然后又由于理论上应该有：F(x)和F_0(x)的n/2项应该相等

所以那个奇怪的减法会导致这个函数没有最后n/2项

于是就发生了神奇的结论！

F(x)-F_0(x)

没有最后n/2项

所以(F(x)-F_0(x))^2没有最后n-1项

于是令人惊喜的是我们泰勒展开的第二项就是(F(x)-F_0(x))^2

所以我们居然只需要展开两项即可！！！woc

所以就得到了神奇的结论：

G(F(x))=G(F_0(x))+G'(F_0(x))*(F(x)-F_0(x))=0(mod~~x^n)

于是就有：

F(x)=F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}(mod~~x^n)

于是就可以倍增处理了...

问题的关键主要是对G'(F_0(x))求导才对吧...

下面给出两个对于G'(F_0(x))求导的例子：

Sixth.多项式\ln/多项式\exp

1.多项式\ln

然后接下来推一下多项式\rm ln

B(x)=\ln(A(x))

两边同时求导得：

令u=A(x)

B'(x)=\ln'(u)*A'(x)

又\ln'(u)=\dfrac{1}{u}

所以得：

B'(x)=\dfrac{A'(x)}{A(x)}

所以就只需要计算出A'(x)即可

因为求导与积分互逆，所以求导可加=积分可加，个人感觉

首先有对一个数求导是0

所以我们知道

A'(x)=\sum_{i=1}^na_i*ix^{i-1}

然后在求一下A(x)的逆A^{-1}(x)

然后将A'(x)和A^{-1}(x)乘起来就好了qwq

然后还原回去就套用上面求导那里的公式

2.多项式\exp

现在想要求解：

B(x)=e^{A(x)}

于是两边同时取\ln得：

\ln B(x)-A(x)=0

可以设G(B(x))=\ln B(x)-A(x)于是问题变成求解它的0点

这个多项式里面A(x)为定值，而B(x)为变量

我们假设有：

G(H(x))=0(mod~~x^{n/2})

问题是求解B(x)

于是代入牛顿迭代的公式得到：

B(x)=H(x)-\dfrac{G(H(x))}{G'(H(x))}

因为G'(H(x))=\ln H(x)-A(x)

而其中A(x)为常量，由于求导的可加性，A'(x)=0，简单求和可以得到：

G'(H(x))=\dfrac{1}{H(x)}

所以我们可以得到：

B(x)=H(x)-(\ln H(x)-A(x))*H(x) B(x)=H(x)(1-\ln H(x)+A(x))

然而因为求解\ln H(x)的复杂度是O(n\log n)的

于是求解\exp的复杂度也是T(x)=x\log x+T(x/2)

即O(n\log n)

不过常数比较大就是了......

Seventh.多项式开根

与上面类似

假设有

B(x)=\sqrt {A(x)}\quad (mod~~x^n)

两边同时平方得到：

B^2(x)-A(x)=0\quad (mod~~x^n)

可以设一个函数G(B(x))=B^2(x)-A(x)

代入牛顿迭代的式子就是：

B(x)=H(x)-\dfrac{G(H(x))}{G'(H(x))}

然后类似的，将A(x)看作常量，那么可以得到：G'(B(x))=2B(x)

所以得到：

B(x)=H(x)-\dfrac{H(x)^2-A(x)}{2H(x)} B(x)=\dfrac{1}{2}H(x)+\dfrac{A(x)}{2H(x)}

于是要求逆，复杂度O(n\log n)

#include<bits/stdc++.h>
using namespace std ;
#define rep( i, s, t ) for( register int i = s; i <= t; ++ i )
#define re register
#define LL long long 
int gi() {
    char cc = getchar() ; int cn = 0, flus = 1 ;
    while( cc < '0' || cc > '9' ) {  if( cc == '-' ) flus = - flus ; cc = getchar() ; }
    while( cc >= '0' && cc <= '9' )  cn = cn * 10 + cc - '0', cc = getchar() ;
    return cn * flus ;
}
const int P = 998244353 ;
const int Gi = 332748118 ;
const int G = 3 ; 
const int N = 5e5 + 5 ; 
int n, R[N], L, limit ;
LL A[N] ; 
LL fpow( LL x, LL k ) {
    LL ans = 1, base = x ;
    while( k ) {
        if( k & 1 ) ans = ( ans * base ) % P ;
        base = ( 1ll * base * base ) % P, k >>= 1 ; 
    }
    return ans % P ; 
} 
LL Gx[N], H[N], Inv ;
void Init( int x ) {
    limit = 1, L = 0 ;
    while( limit < x ) limit <<= 1, ++ L ; //因为都是2的整数次幂，实际上并不需要 <=
    //这样会造成8倍长度的数组，最后TLE 
    rep( i, 0, limit - 1 ) R[i] = ( R[i >> 1] >> 1 ) | ( ( i & 1 ) << ( L - 1 ) ) ;
    Inv = fpow( limit, P - 2 ) ;
}
void NTT( LL *a, int type ) {
    rep( i, 0, limit - 1 ) if( R[i] > i ) swap( a[i], a[R[i]] ) ;
    for( re int k = 1; k < limit; k <<= 1 ) {
        int d = fpow( ( type == 1 ) ? G : Gi, ( P - 1 ) / ( k << 1 ) ) ;
        for( re int i = 0; i < limit; i += ( k << 1 ) )
            for( re LL j = i, g = 1; j < i + k; ++ j, g = ( d * g ) % P ) {
                LL Nx = a[j], Ny = ( a[j + k] * g ) % P ;
                a[j] = ( Nx + Ny ) % P, a[j + k] = ( Nx - Ny + P ) % P ; 
            }
    }
    if( type != 1 ) rep( i, 0, limit - 1 ) a[i] = ( a[i] * Inv ) % P ;
}
void GetInv( LL *a, int x ) {
    int k = 1 ; memset( H, 0, sizeof(H) ), memset( Gx, 0, sizeof(Gx) ), H[0] = fpow( A[0], P - 2 ) ;
    while( k < x ) { //这里同理 
        k <<= 1, Init( k * 2 ) ;
        rep( i, 0, k - 1 ) Gx[i] = a[i] ;
        NTT( Gx, 1 ), NTT( H, 1 ) ;
        rep( i, 0, limit ) H[i] = ( 2ll - Gx[i] * H[i] % P + P ) * H[i] % P ; 
        NTT( H, -1 ) ; rep( i, k, limit ) H[i] = 0 ;
    } 
}
LL Gt[N], B[N] ; 
void Sqrt( LL *a, int x ) {
    int k = 1 ; LL Iv = 499122177 ; B[0] = 1 ; 
    while( k <= x ) {
        k <<= 1 ; rep( i, 0, k - 1 ) Gt[i] = a[i] ; 
        GetInv( B, k ), Init( k * 2 ) ;
        NTT( Gt, 1 ), NTT( H, 1 ), NTT( B, 1 ) ;
        rep( i, 0, limit ) B[i] = 1ll * Iv * ( B[i] % P + Gt[i] * H[i] % P ) % P ;
        NTT( B, -1 ) ; rep( i, k, limit ) B[i] = 0 ; 
    }
}
signed main()
{
    n = gi() ;
    rep( i, 0, n - 1 ) A[i] = gi() ; 
    Sqrt( A, n ) ;
    rep( i, 0, n - 1 ) printf("%d ", B[i] ) ;
    return 0 ;
}