多项式入门

在数学中，由若干个单项式相加组成的代数式叫做多项式（若有减法：减一个数等于加上它的相反数）。多项式中的每个单项式叫做多项式的项，这些单项式中的最高项次数，就是这个多项式的次数。其中多项式中不含字母的项叫做常数项。
在数学中，多项式（polynomial）是指由变量、系数以及它们之间的加、减、乘、幂运算（非负整数次方）得到的表达式。对于比较广义的定义，一个或零个单项式的和也算多项式。按这个定义，多项式就是整式。实际上，还没有一个只对狭义多项式起作用，对单项式不起作用的定理。$0$ 作为多项式时，次数定义为负无穷大（或 $0$）。单项式和多项式统称为整式。

多项式的定义初中就有，没必要多说吧。本文将会积累所有学了的有关多项式的知识，不定期更新内容。

多项式求积

对于一个 $n$ 项和一个 $m$ 项的多项式 $A,B$，如何求它们的积 $C=A\ast B$？考虑如果暴力求解，也就是卷积：

$C = \sum_{i=0}^{n+m-2} \sum_{j=0}^i a_jb_{i-j}x^i$

复杂度是 $n^2$ 的。不能接受。直接抛弃系数表示法，尝试用点值表示法，显然如果变成了点值表示法，卷积就是把每个点值乘起来。复杂度变成 $\mathcal O(n)$ 了。

随意取点把多项式的系数表达式转化成点值表达式是 $n^2$ 复杂度的。这就和最开始没有区别了，因此我们需要一种方式加速这样的转换。

FFT

快速傅里叶变换，是一种把系数表达式转化为点值表达式的方式。首先，为了方便，把 $n$ 转化为 $2$ 的幂，如果 $n$ 不是二的幂，则补充 $0$ 作为高次项系数补成 $2$ 的幂。

考虑选出一些 $x$ 求出多项式的点值的时候，如果能满足 $\exists p,x^p = 1$，那么乘方运算的常数将会有下降。而且这样的数还有一些特殊性质。如果放眼纯虚数或者实数，只能找到四个这样的数 $\pm1,\pm i$。然而如果放在复平面上，半径为 $1$ 圆心为原点的圆上所有点都满足。

于是把这个圆 $n$ 等分掉：

从 $X_0 = (1,0)$ 开始逆时针标号，记 $\omega_n^k$ 表示 $X_k$ 对应的点值。不难发现 $(\omega_n^1)^k = \omega_n^k$。定义 $w_n^1$ 为单位根。

考虑 $\omega$ 的性质：

• $(\omega_n^p)^k = \omega_n^{pk}$。
• $\omega_n^n = \omega_n^0 = 1$。
• $\omega_k^n = \omega_{2n}^{2k}$。
• $\omega_n^{k+\frac{n}{2}} = -\omega_n^k$。
• $\omega_n^k\omega_n^p = \omega_n^{k+p}$。

前两个显然，第三个相当于把 $n$ 等分圆变成 $2n$ 等分，第四个相当于旋转 $180^\circ$，第五个相当于辐角相加。

把 $n$ 个 $\omega$ 当作 $x$ 代进去。用分治的思想，按照指数奇偶性分成两个部分：

$(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2}) + x (a_1+a_3x^2+\cdots + a_{n-1}x^{n-2})$

这两个部分长得差不多，拆成两个多项式 $F_1,F_2$：

$F_1(x) = a_0 + a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}$

$F_2(x) = a_1 + a_3x + a_5x^2 + \cdots + a_{n-1}x^{\frac{n}{2}-1}$

于是有 $F=F_1(x^2)+xF_2(x^2)$。对于 $k<\frac{n}{2}$，把 $x=\omega_n^k$ 代入：

$\begin{aligned} A_1(x^2) + xA_2(x^2) &= A_1(\omega_n^{2k}) + \omega_n^k A_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^k) + \omega_n^k A_2(\omega_{\frac{n}{2}}^k) \end{aligned}$

另一半是 $x=\omega_n^{k+\frac{n}{2}}$：

$\begin{aligned} A_1(x^2) + xA_2(x^2) &= A_1((-\omega_n^k)^2) - \omega_n^k A_2(-(\omega_n^k)^2) \\ &= A_1(\omega_n^{2k}) - \omega_n^k A_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^k) - \omega_n^k A_2(\omega_{\frac{n}{2}}^k) \end{aligned}$

只有第二项的符号不一样。因此如果知道 $A_1(\omega_{\frac{n}{2}}^k),A_2(\omega_{\frac{n}{2}}^k),\omega_n^k$，就能知道 $A(\omega_n^k),A(\omega_n^{k+\frac{n}{2}})$。考虑递归，每次算出 $k\in[1,\frac{n}{2}]$ 的结果就能推出 $k\in(\frac{n}{2},n]$。边界条件 $k=1$，也就是单位根。递归时问题序列每次减半，时间复杂度是 $n\log n$。

只能转化成点值没有用啊，还要考虑怎么从点值转回来。如果只能拉插那不还是 $n^2$ 吗，废了废了。

下面不妨考虑 $\omega_n^k$ 中 $k<0$ 的情况。可以认为是顺时针编号，有 $\omega_n^{-k} = \omega_n^{n-k}$。

考虑怎么把一个点值表达式转化为系数表达式。假设我们点值表达式得到的每个点的点值提出来，组成 $g_i$，系数表达式的每一个系数也提出来，组成 $a_i$。

结论

$na_k = \sum_{p=0}^{n-1} g_p(\omega_n^{-k})^p$

证明

根据定义 $g_i = \sum_{k=0}^{n-1} a_k(\omega_n^i)^k$，代入：

$\begin{aligned}\sum_{p=0}^{n-1} g_p(\omega_n^{-k})^p &= \sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q(\omega_n^p)^q(\omega_n^{-k})^p \\&= \sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q\omega_n^{p(q-k)}\end{aligned}$

分类讨论 $q,k$ 的关系。$q=k$ 时结果为 $\sum_{p=0}^{n-1} a_k$，显然为 $na_k$。否则设 $r=q-k$：

$\begin{aligned}\sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q\omega_n^{pr} &= a_{r+k} \sum_{p=0}^{n-1} \omega_n^p\omega_n^r \\&= a_{r+k} \omega_n^r\sum_{p=0}^{n-1} \omega_n^p \\&= a_{r+k}\omega_n^r \frac{\omega_n^n-1}{\omega_n^1-1} \\&= 0\end{aligned}$

两部分加起来得证。

考虑这个性质，就是原来系数转点值时把 $\omega_n^i$ 换成 $\omega_n^{-i}$，最后再除上 $n$。

上面是 FFT 的理论，但是实践的时候，会发现递归实现常数很大。于是我们需要把递归转成循环。

我们考虑分治的时候怎么分治的：我们把下标（从 $0$ 开始）是偶数的提出来放在前面，下标是奇数的提出来放在后面，然后往下递归。考虑原数组下标（下标组一）和分治到最终位置后原数组的下标跑到的位置（下标组二），提出二进制：

下标组一	0	1	2	3	4	5	6	7
下标组二	0	4	2	6	1	5	3	7
下标组一二进制	000	001	010	011	100	101	110	111
下标组二二进制	000	100	010	110	001	101	011	111

每个位置分治后的最终位置为它的二进制倒过来得到的数字对应的位置。

预处理每个下标的最终位置 $rev$ 即可。一种简单的求法为：

rev[i]=(rev[i>>1]>>1) | ((i&1)<<(lim-1))
// lim 表示最高位

这样就不需要递归实现 FFT 了，这样的优化叫做蝴蝶变换优化。

继续卡常数。如果要求 $A\ast B = C$ 的话，需要对 $A$ FFT，$B$ FFT，求出 $C$ 之后再 FFT（转回系数表示）。需要三次 FFT。然而设复数 $z=a+bi$，有 $z^2 = a^2-b^2+2abi$。这后面有一个 $2ab$，因此把 B 的系数当作 A 的系数的虚部，得到多项式 $A'$，就只需要对 $A'$ FFT，然后让其自乘，取出系数的虚部除以 $2$ 就是 $C$ 的系数了。

这样就从三次 FFT 变成了两次，这个优化称为三次变两次优化。

NTT

FFT 非常强大，但是有一个显著的缺点：需要实现复数计算。然而如果需要结果对一个质数取模，FFT 就无法简单进行了。其次，FFT 的计算基于 double 类型，是拥有精度误差的，double 也比 int 之类的整数类型常数更大。

基于上面考虑，就有了专门在模意义下操作的 FFT：快速数论变换 NTT。这样的题目一般会给我们很好的模数：$998244353$。为什么选择它？

考虑 NTT 是如何避免复数计算的——它使用了模数的原根代替单位根。然而如果需要用原根代替单位根，则需要原根也满足单位根的那些在 FFT 里利用过的性质：

• $(\omega_n^1)^k = \omega_n^k$。
• $\omega_n^k = \omega_{2n}^{2k}$。
• $\omega_n^{k+n} = \omega_n^k$。
• $\omega_n^{0\sim(n-1)}$ 互不相同。

其他 FFT 里写的性质都可以用这四条推出。由第一点、第三点、第四点，综合起来就等价于：$\omega_n^1$ 在模意义下阶恰好为 $n$。设题目给出模数 $p$，$g$ 是其一个原根。显然 $g$ 的阶也就是 $p-1$ 不一定为 $n$，因此 $g$ 不能代替单位根。然而阶有这样的性质：$g^k$ 的阶为 $\frac{p-1}{\gcd(k,p-1)}$。考虑控制 $k$ 的值使得阶数取到 $n$，不难计算出 $k=\frac{p-1}{n}$。

于是 $g^k$ 就满足了性质的第一、三、四点，考虑证明第二点满足。注意到 $\omega_n^k = \omega_{2n}^{2k}$ 等价于求证 $(\omega_{2n}^1)^2 = \omega_n^1$，直接把 $g^k$ 代入就发现成立了。

因此一个题能 NTT，需要 $n\mid (p-1)$，此时可以用 $g^{\frac{p-1}{n}}$ 取代单位根。回到最开始的问题：$998244353$ 为什么适合 NTT，是因为 $998244352=2^{23}\times 7\times 17$，而我们最开始会把 $n$ 扩充成 $2$ 的幂。因此在正常数据范围里它就是 $n$ 的倍数。

注意因为没有了实数和虚数部分，所以没有三次变两次优化，只有蝴蝶变换优化。NTT 常数很小，有图为证：

MTT

如果模数并不是友好的 $998244353$ 或其他 $p-1$ 因数有 $2$ 的大整数幂的时候，就无法用 NTT 了。而且因为要求取模，FFT 也做不了。

这时可以选择三个友好的模数，跑九次 NTT，最后的答案用中国剩余定理合并。然而一旦没选好就要写高精，而且常数太大，无法接受。

考虑如何用 FFT 的方法去做，常见的方法是拆系数。根据值域范围拆。举个例子，如果值域 $10^9$，就把系数 $a$ 拆成 $a=2^{15}p+q$。可以发现 $p,q\leq 2^{15}$，计算中就不会爆范围了。然后考虑拆系数之后：

$a_1a_2 = (2^{15}p_1+q_1)(2^{15}p_2+q_2) = 2^{30}p_1p_2 + 2^{15}p_1q_2 + 2^{15} p_2q_1 + q_1q_2$

原本一次 FFT 要拆成四次 FFT，而且因为虚部被占用了不能三次变两次，总共要 $12$ 次 FFT，还是不可接受。

考虑现在有四个多项式 $A_1,A_2,B_1,B_2$，多项式系数都是实数。令多项式 $P,Q$，满足 $P$ 的实部是 $A_1$，虚部是 $A_2$， $Q$ 同理。

$C_1 = P\ast Q = A_1B_1 - A_2B_2 + A_1B_2i + A_2B_1 i$

再开一个 $R$，实部是 $A_1$，虚部是 $-A_2$。

$C_2 = R\ast Q = A_1B_1 + A_2B_2 + A_1B_2 i - A_2B_1 i$

然后有：

$C_1+C_2 = 2(A_1B_1 + A_1B_2 i)$

$C_1-C_2 = 2(A_2B_1 i - A_2B_2)$

分别取出这两个结果的实部和虚部就可以得到四个多项式两两的乘积分别是多少。代价是需要求出 $P,Q,R$ 三个多项式的 $P\ast Q,R\ast Q$，后面还需要还原 $C_1,C_2$。八次 FFT。

进一步，把 $A_1,B_1,A_2,B_2$ 当成拆系数后得到的四个多项式，然后利用虚实部转化为 $P,Q,R$，求出 $C_1,C_2$。然后提取 $a_1b_1,a_1b_2,a_2b_1,a_2b_2$。代回最前面的式子就求出答案了。

需要进行三次 FFT，两次 IFFT，常数不算那么难看了。

多项式求逆

设 $A$ 模 $x^n$ 意义下的逆元为 $B_1$，模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 意义下的逆元为 $B_2$，有：

$A\ast B_1 \equiv 1 \pmod {x^n}$

$A\ast B_2 \equiv \pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$

上面两个式子是这样的形式：$ab=xp^2+1,ac=yp+1$，于是有 $ab-ac=p(xp-y)$，即 $a(b-c)\equiv 0\pmod p$。于是 $b-c\equiv 0 \pmod p$：

$B_1-B_2\equiv 0 \pmod {x^{\left\lceil\frac{n}{2}\right\rceil}}$

如果 $p$ 能被 $a$ 整除，则 $p^2$ 能被 $a^2$ 整除。因此：

$A\ast (B_1^2+B_2^2-2B_1B_2) \equiv 0 \pmod{x^n}$

注意到 $A\ast B_1 \equiv 0\pmod{x^n}$：

$B_1\equiv 2B_2 - AB_2^2 \equiv 0 \pmod{x^n}$

也就是求 $A$ 的乘法逆 $B_1$，可以通过求出让 $A\ast B_2\equiv 1 \pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$ 的 $B_2$，然后多项式乘法求。

主体是 NTT 和倍增。复杂度是 $n\log n$ 的。

多项式求导积分

$x^{a'} = ax^{a-1}$

$\int x^a\, \mathrm{dx} = \frac{1}{a+1} x^{a+1}$

多项式 ln,exp

先考虑 $\ln$，设：

$G(x) = \ln F(x)$

两边求导：

$G'(x) = \ln'(F(x)) F'(x)$

$\ln'(x) = \frac{1}{x}$：

$G'(x) = \frac{F'(x)}{F(x)}$

然后就可以用求导积分和多项式逆元做了。下面考虑 exp，这个比较复杂：

$G(x) = e^{F(x)}$

则：

$\ln G(x) - F(x) = 0$

令 $A(G(x)) = \ln G(x) - F(x) = 0$，还是导：

$(A(G_0(x)))' = \frac{G_0'(x)}{G_0(x)}$

代入回去：

$G(x) = G_0(x) - \frac{\ln G_0(x) - F(x)}{\frac{G_0'(x)}{G_0(x)}}$

化简：

$G(x) = \frac{G_0(x)(1-\ln G_0(x) + F(x)}{G_0'(x)}$

然后每一次递归求一遍多项式求逆，做一次多项式 $\ln$，跑一遍多项式乘法。复杂度 $n\log n$，常数可想而知得大。

多项式快速幂

先讨论一个特殊情况，$F(0)=1$：

$G(x) = F^k(x)$

两边对数：

$\ln G(x) = k\ln F(x)$

然后 $\exp$ 回去：

$G(x) = e^{k\ln F(x)}$

对于 $F(0)\neq 1$ 的情况，先让 $F$ 降次再升次回来即可。复杂度 $n\log n$。

多项式开根

$G^2(x) - F(x) = 0$

设 $A(G(x)) = G^2(x) - F(x)$，则：

$A'(G(x)) = 2G(x)$

牛顿迭代：

$G(x) = G_0(x) - \frac{A(G_0(x))}{A'(G_0(x))}$

于是有：

$G(x) = \frac{F(x) + G_0^2(x)}{2G_0(x)}$

多项式除法

这里的多项式除法，不是乘逆元，而是求出商多项式和余数多项式，例如 $F = QG+R$，$F,G$ 为给定的 $n,m$ 次多项式。设存在一种操作 $r$ 满足：

$A^r(x) = x^n A(\frac{1}{x})$

不难有 $A^r_x = A^r_{n-x}$。因此这样的操作可以 $\mathcal O(n)$ 完成。因此满足原条件，可以视为满足：

$F(\frac{1}{x}) = Q(\frac{1}{x})\ast G(\frac{1}{x}) + R(\frac{1}{x})$

因此：

$x^nF(\frac{1}{x}) = x^{n-m} Q(\frac{1}{x}) \ast x^mG(\frac{1}{x}) + x^{n-m+1}\ast x^{m-1} R(\frac{1}{x})$

用 $r$ 来表示可以得到：

$F^r(x) = Q^r(x)\ast G^r(x) + x^{n-m+1}\ast R^r(x)$

于是有：

$F^r(x) = Q^r(x) \ast G^r(x)$

求出 $G^r$ 的逆，乘上 $F^r$ 得到 $Q^r$，就算出来 $Q$ 了，代回去可以算出 $R$。复杂度 $n\log n$。

代码

NTT 全家桶

包含 NTT，求逆，快速幂，开根，$\ln$，$\exp$，求导，积分：

namespace DXS{
    int rev[max_n*8],iv[max_n*8];
    int F[max_n*8],G[max_n*8],F1[max_n*8],G1[max_n*8],F2[max_n*8];
    const int _G=3,IG=332748118;

    inline void init(int n){
        for(register int i=0;i<=n*2;++i)
            iv[i]=mi(i);
    }
    inline void NTT(int *g,int n,bool op){
        for(register int i=0;i<n;++i){
            rev[i]=(rev[i>>1]>>1)|((i&1)?n>>1:0);
            if(i<rev[i])
                g[i]^=g[rev[i]]^=g[i]^=g[rev[i]];
        }
        for(register int mid=1;mid<n;mid<<=1){
            int omega=mi(op?_G:IG,(mod-1)/(mid+mid));
            for(register int i=0;i<n;i+=mid+mid){
                int w=1;
                for(register int j=0;j<mid;++j,w=(w*omega)%mod){
                    int x=g[i+j],y=w*g[i+j+mid]%mod;
                    g[i+j]=(x+y)%mod,
                    g[i+j+mid]=(x-y+mod)%mod;
                }
            }
        }
        if(!op) for(register int i=0;i<n;++i)
            g[i]=g[i]*iv[n]%mod;
    }
    inline void Mul(int *f,int *g,int n,int m){
        int tmp=1;
        while(tmp<=n+m) tmp<<=1;
        init(tmp),NTT(f,tmp,1),NTT(g,tmp,1);
        for(register int i=0;i<tmp;++i)
            f[i]=f[i]*g[i]%mod;
        NTT(f,tmp,0);
    }
    inline void Inv(int *f,int *g,int n){
        int len=1;g[0]=mi(f[0]);
        for(;len<(n<<1);len<<=1){
            int lim=len<<1;
            for(register int i=0;i<len;++i)
                F1[i]=f[i],G1[i]=g[i];
            NTT(F1,lim,1),NTT(G1,lim,1);
            for(register int i=0;i<lim;++i)
                g[i]=((2-F1[i]*G1[i])%mod*G1[i]%mod+mod)%mod;
            NTT(g,lim,0);
            for(register int i=len;i<lim;++i)
                g[i]=0;
        }
        for(register int i=0;i<len;++i)
            F1[i]=G1[i]=0;
        for(register int i=n;i<len;++i)
            g[i]=0;
    }
    inline void Dir(int *f,int *g,int n){
        for(register int i=1;i<n;++i)
            g[i-1]=f[i]*i%mod;
        g[n-1]=0;
    }
    inline void Int(int *f,int *g,int n){
        for(register int i=1;i<n;++i)
            g[i]=f[i-1]*iv[i]%mod;
        g[0]=0;
    }
    inline void Ln(int *f,int *g,int n){
        Dir(f,F,n),Inv(f,G,n);
        int len=n<<1;
        NTT(F,len,1),NTT(G,len,1);
        for(register int i=0;i<len;++i)
            F[i]*=G[i],F[i]%=mod;
        NTT(F,len,0),Int(F,g,len);
        for(register int i=0;i<len;++i)
            F[i]=G[i]=0;
    }
    inline void Exp(int *f,int *g,int n){
        if(n==1){
            g[0]=1;
            return;
        }
        int len=n<<1;
        Exp(f,g,n>>1),Ln(g,F2,n);
        F2[0]=(f[0]+1-F2[0]+mod)%mod;
        for(register int i=1;i<n;++i)
            F2[i]=(f[i]-F2[i]+mod)%mod;
        NTT(F2,len,1),NTT(g,len,1);
        for(register int i=0;i<len;++i)
            g[i]*=F2[i],g[i]%=mod;
        NTT(g,len,0);
        for(register int i=n;i<len;++i)
            g[i]=F2[i]=0;
    }
    inline void Pow(int *f,int *g,int n,int k){
        int tmp=1;
        while(tmp<=n) tmp<<=1;
        init(tmp),Ln(f,g,tmp);
        for(register int i=0;i<n;++i)
            f[i]=g[i]*k%mod;
        for(register int i=0;i<=tmp*2;++i)
            g[i]=0;
        Exp(f,g,tmp);
    }
    inline void Sqrt(int *f,int *g,int n){
        int len=1;g[0]=1;
        for(;len<(n<<1);len<<=1){
            int lim=len<<1;
            for(register int i=0;i<len;++i)
                F[i]=f[i];
            Inv(g,G,len);
            NTT(F,lim,1),NTT(G,lim,1);
            for(register int i=0;i<lim;++i)
                F[i]=F[i]*G[i]%mod;
            NTT(F,lim,0);
            for(register int i=0;i<len;++i)
                g[i]=(g[i]+F[i])*iv[2]%mod;
            for(register int i=len;i<lim;++i)
                g[i]=0;
        }
        for(register int i=0;i<len;++i)
            F[i]=G[i]=0;
        for(register int i=n;i<len;++i)
            g[i]=0;
    }
    inline void Div(int *f,int *g,int n,int m){
        int tmp=1,rmp=1;
        while(tmp<=n+m) tmp<<=1;
        while(rmp<=n-m) rmp<<=1;
        init(tmp);
        for(register int i=0;i<=n;++i)
            F[i]=f[n-i];
        for(register int i=0;i<=m;++i)
            F2[i]=g[m-i];
        for(register int i=n-m+1;i<=m;++i)
            F2[i]=0;
        Inv(F2,G,n-m+1),Mul(F,G,n,n-m);
        for(register int i=0;i<tmp;++i)
            G1[i]=f[i],F2[i]=G[i]=f[i]=0;

        for(register int i=0;i<=n-m;++i)
            F2[i]=f[i]=F[n-m-i];
        for(register int i=0;i<=m;++i)
            G[i]=g[i],g[i]=0;
        Mul(G,F2,m,n-m);
        for(register int i=0;i<m;++i)
            g[i]=(G1[i]-G[i]+mod)%mod;
        for(register int i=0;i<tmp;++i)
            G[i]=F1[i]=F[i]=F2[i]=G1[i]=0;
    }
}

乘法套餐

包含 FFT，NTT，MTT：

namespace DXS{
	inline int mi(int x,int p){
		int res=1;
		while(p){
			if(p&1) res*=x,res%=mod;
			x*=x,x%=mod,p>>=1;
		}
		return res;
	}
	struct com{
		double a,b; // a+bi
		com(double x=0,double y=0){a=x,b=y;}
		com operator + (const com b) const{
			com a=*this,res;
			res.a=a.a+b.a,res.b=a.b+b.b;
			return res;
		}
		com operator - (const com b) const{
			com a=*this,res;
			res.a=a.a-b.a,res.b=a.b-b.b;
			return res;
		}
		com operator * (const com b) const{
			com a=*this,res;
			res.a=a.a*b.a-a.b*b.b,
			res.b=a.a*b.b+a.b*b.a;
			return res;
		}
	};

	int rev[max_n*2];
	const int mod=998244353,G=3,IG=332748118;
	const double Pi=acos(-1.0);

	inline void NTT(int *g,int n,bool op){
		for(register int i=0;i<n;++i) if(i<rev[i])
			swap(g[i],g[rev[i]]);
		for(register int mid=1;mid<n;mid<<=1LL){
			int omega=mi(op?G:IG,(mod-1)/(mid+mid));
			for(register int i=0;i<n;i+=mid+mid){
				int w=1;
				for(register int j=0;j<mid;++j,w=(w*omega)%mod){
					int x=g[i+j],y=w*g[i+j+mid]%mod;
					g[i+j]=(x+y)%mod,
					g[i+j+mid]=(x-y+mod)%mod;
				}
			}
		}
		if(!op){
			int iv=mi(n,mod-2);
			for(register int i=0;i<n;++i)
				g[i]=g[i]*iv%mod;
		}
	}
	inline void FFT(com *g,int n,bool op){
		for(register int i=0;i<n;++i) if(i<rev[i])
			swap(g[i],g[rev[i]]);
		for(register int mid=1;mid<n;mid<<=1LL){
			com omega(cos(Pi/mid),sin(Pi/mid)*(op?1:-1));
			for(register int i=0;i<n;i+=mid+mid){
				com w(1,0);
				for(register int j=0;j<mid;++j,w=w*omega){
					com x=g[i+j],y=w*g[i+j+mid];
					g[i+j]=x+y,
					g[i+j+mid]=x-y;
				}
			}
		}
	}
	inline void mulNTT(int *f,int *g,int n,int m){
		int len=1;
		while(len<n+m) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
		NTT(f,len,1),NTT(g,len,1);
		for(register int i=0;i<len;++i)
			f[i]*=g[i],f[i]%=mod;
		NTT(f,len,0);
	}
	inline void mulFFT(int *f,int *g,int n,int m){
		static com F[max_n],G[max_n];
		for(register int i=0;i<n;++i)
			F[i].a=f[i],F[i].b=0;
		for(register int i=0;i<m;++i)
			G[i].a=g[i],G[i].b=0;
		int len=1;
		while(len<n+m) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
		FFT(F,len,1),FFT(G,len,1);
		for(register int i=0;i<len;++i)
			F[i]=F[i]*G[i];
		FFT(F,len,0);
		for(register int i=0;i<len;++i)
			f[i]=(int)(F[i].a/len+0.5);
	}
	const int K=32768,S=K*K;
	inline void mulMTT(int *f,int *g,int n,int m,int mod){
		static com P[max_n],Q[max_n],R[max_n];
		for(register int i=0;i<n;++i){
			P[i].a=f[i]/K,
			P[i].b=f[i]%K;
			R[i].a=f[i]/K,
			R[i].b=-f[i]%K;
		}
		for(register int i=0;i<m;++i){
			Q[i].a=g[i]/K,
			Q[i].b=g[i]%K;
		}
		int len=1;
		while(len<n+m) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);

		FFT(P,len,1),FFT(Q,len,1),FFT(R,len,1);
		for(register int i=0;i<len;++i){
			Q[i].a/=len,
			Q[i].b/=len;
		}
		for(register int i=0;i<len;++i){
			P[i]=P[i]*Q[i],
			R[i]=R[i]*Q[i];
		}
		FFT(P,len,0),FFT(R,len,0);
		for(register int i=0,a1b1=0,a1b2=0,a2b1=0,a2b2=0;i<len;++i){
			a1b1=(int)floor((P[i].a+R[i].a)/2+0.5)%mod,
			a1b2=(int)floor((P[i].b+R[i].b)/2+0.5)%mod,
			a2b1=((int)floor(P[i].b+0.5)-a1b2)%mod,
			a2b2=((int)floor(R[i].a+0.5)-a1b1)%mod;
			f[i]=(a1b1*S%mod+(a1b2+a2b1)%mod*K%mod+a2b2)%mod;
			while(f[i]<0) f[i]+=mod;
			f[i]%=mod;
		}
	}
}