卷积

$$ f[k]=\sum_{i\oplus j=k} a[i]\times b[j] $$

系数表达式 $\to$ 频域

点值表达式 $\to$ 时域

时域卷积，频域乘积，频域卷积，时域乘积。

$O(Nlog N)$ 计算卷积。

快速傅里叶变换

用主 $n$ 次复数单位根 $\omega_n^k$ 插值。

$$ \omega_n^k=(e^{\frac{2\pi i}{n}})^k=\cos k\frac{2\pi}{n}+i \sin k\frac{2\pi}{n} $$ 消去引理：

$$ \omega_{dn}^{dk}=(e^{\frac{2\pi i}{dn}})^{dk}=(e^{\frac{2\pi i}{n}})^k=\omega_n^k $$ 折半引理：

$$ \omega_n^{k+\frac{n}{2}}=(\cos(\frac{n}{2}\times \frac{2\pi}{n})+i\sin(\frac n 2\times \frac{2\pi} 2))\times \omega_n^k=(\cos\pi+i\sin\text{ }\pi)\times \omega_n^k=-\omega_n^k $$ 求和引理：

$$ \sum_{i=0}^{n-1}(\omega_n^k)^i=0(k\ne0) $$ 对系数奇偶分类

$$ \begin{align} &A(x)=(a_0+a_2x^2+a_4x^4+⋯+a_{n−2}x^{n−2})+(a_1x+a_3x^3+a_5x^5+⋯+a_{n−1}x^{n−1})\\ &A_1(x)=a_0+a_2x+a_4x^2+⋯+a_{n−2}x^{\frac n 2−1}\\ &A_2(x)=a_1+a_3x+a_5x^2+⋯+a_{n−1}x^{\frac n 2−1} \end{align} $$ $k<\frac n 2$

$$ A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2)=A_1(\omega_\frac n 2 ^k)+\omega_n^kA_2(\omega_\frac n 2 ^k) $$ $k>\frac n 2$

$$ \begin{align} &A(\omega_n^{k+\frac n 2})\\ =&A_1((\omega_n^{k+\frac n 2})^2)+\omega_n^{k+\frac n 2}A_2((\omega_n^{k+\frac n 2})^2)\\ =&A_1(\omega_n^{2k}\times\omega_n^k)-\omega_n^kA_2(\omega_n^{2k}\times\omega_n^n)\\ =&A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\\ =&A_1(\omega_\frac n 2^k)-\omega_n^kA_2(\omega_\frac n 2^k) \end{align} $$ 在外层循环 $i$ ， $F[p]$ 表示分成的段中 $\omega_i^k$ 在当前这个点插值的值，刚开始是多项式系数，即为 $\omega_1^1$ 在各个点插值的值，最后就是点值表达式，即为 $\omega_n^p$ 在各个点插值的值。

逆 $FFT$ ：

已知 $\omega_n^k$ 插值得到的点值表达式 $B$ ，求系数表达式 $A$ 使得 $A(\omega_n^k)=B_k$ 。

用 $-\omega_n^k$ 插值：

$$ C_k=\sum_{i=0}^{n-1}B_k(\omega_n^{-k})^i=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}A_j(\omega_n^i)^j)(\omega_n^{-k})^i=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}A_j(\omega_n^{j-k})^i=\sum_{j=0}^{n-1}A_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i $$ 由求和引理得：当 $k\ne 0$ 时： $$ \sum_{i=0}^{n-1}(\omega_n^k)^i=0 $$ 否则： $$ \begin{align}\sum_{i=0}^{n-1}(\omega_n^k)^i=n\end{align} $$ 所以 $C_k=n\times A_k$ ，即 $A_k=C_k/n$ 。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

const double PI = acos(-1);

int n,m;

#define MAXN 3100013

int rev[MAXN];

struct complex

{

	double r,i;

	complex(double r_ = 0.0,double i_ = 0.0){r = r_;i = i_;}

}a[MAXN],b[MAXN];

complex operator + (complex a,complex b){return complex(a.r + b.r,a.i + b.i);}

complex operator - (complex a,complex b){return complex(a.r - b.r,a.i - b.i);}

complex operator * (complex a,complex b){return complex(a.r * b.r - a.i * b.i,a.r * b.i + a.i * b.r);}

void FFT(complex *f,int l,int type)

{

	for(int i = 0;i < l;++i)

    {

        if(i < rev[i])swap(f[i],f[rev[i]]);

    }

	for(int i = 2;i <= l;i = i << 1)

	{

		complex wn(cos(-type * 2 * PI / i),sin(-type * 2 * PI / i));

		for(int j = 0;j < l;j += i)

		{

			complex w(1,0);

			for(int k = j;k < j + i / 2;++k)

			{

				complex u = f[k],t = w * f[k + i / 2];

				f[k] = u + t;

				f[k + i / 2] = u - t;

				w = w * wn;

			}

		}

	}

	if(type == -1)

	{

		for(int i = 0;i < l;++i)

		{

			f[i].r /= l;

		}

	}

	return;

}

int main()

{

	scanf("%d%d",&n,&m);

	for(int i = 0;i <= n;++i)scanf("%lf",&a[i].r);

	for(int i = 0;i <= m;++i)scanf("%lf",&b[i].r);

	int len,l;

	for(len = 1,l = 0;len <= n + m;len = len << 1,++l);

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	FFT(a,len,1);

	FFT(b,len,1);

	for(int i = 0;i < len;++i)

	{

		a[i] = a[i] * b[i];

	}

	FFT(a,len,-1);

	for(int i = 0;i <= n + m;++i)

	{

		printf("%d ",(int)(a[i].r + 0.5));

	}

	return 0;

}

快速数论变换

计算卷积对 $P=r\times2^k+1$ 取模的值。

原根的幂在 $mod\text{ }p$ 意义下有循环性质。

将 $g^{\frac{P-1}{n}}$ 看作 $e^{-\frac{2\pi i}{n}}$ 的等价

设 $\omega_n^k=(g^{\frac{P-1}{n}})^k$

消去引理：

$\omega_{dn}^{dk}=(g^{\frac{P-1}{dn}})^{dk}=(g^{\frac{P-1}{n}})^k=\omega_n^k$

折半引理：

$\omega_n^{k+\frac{n}{2}}=(g^{\frac{P-1}{n}})^{k+\frac n 2}=(g^{\frac{P-1}{n}})^k\times g^{\frac{P-1}{2}}$

由于 $g^{\frac{P-1} 2}\equiv-1(mod\text{ }P)$

$\omega_n^{k+\frac n 2}=-g^{\frac{P-1} 2}=-\omega_n^k$

求和引理：

$\begin{align}\sum_{i=0}^{n-1}(\omega_n^k)^i=\frac{(\omega_n^k)^n-1}{\omega_n^k-1} =\frac{((g^\frac {P-1}{n})^k)^n-1}{(g^{\frac{P-1}{n}})^k-1}=\frac{(g^{P-1})^k-1}{\omega_n^k-1}\end{align}$

由费马小定理得： $g^{P-1}\equiv1(mod\text{ }P)$

$\begin{align}= \frac{1-1}{\omega_n^k-1}=0\end{align}$

即 $\begin{align}\sum_{i=0}^{n-1}(\omega_n^k)^i=0\end{align}$

于是就和 $FFT$ 一样了。

逆 $NTT$ 时也应除掉 $n$ ，计算一下 $n$ 对 $MOD$ 的逆元，即 $n^{MOD-2}$ 。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

typedef long long ll;

const ll P = 998244353;

int n,m;

#define MAXN 4000010

ll a[MAXN],b[MAXN];

int rev[MAXN];

ll ww[MAXN << 1],*g = ww + MAXN;

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = (res * a) % P;

		a = (a * a) % P;

		b = b >> 1;

	}

	return res;

}

void NTT(ll *f,int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		ll wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			ll w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				ll u = f[k],t = (w * f[k + i / 2]) % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = (w * wn) % P;

			}

		}

	}

	if(type == -1)

	{

		ll ni = power(l,P - 2);

		for(int i = 0;i < l;++i)

		{

			f[i] = (f[i] * ni) % P;

		}

	}

	return;

}

int main()

{

	scanf("%d%d",&n,&m);

	for(int i = 0;i <= n;++i)scanf("%lld",&a[i]);

	for(int i = 0;i <= m;++i)scanf("%lld",&b[i]);

	int len,l;

	for(len = 1,l = 0;len <= n + m;len = len << 1,++l);

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	NTT(a,len,1);

	NTT(b,len,1);

	for(int i = 0;i < len;++i)

	{

		a[i] = a[i] * b[i] % P;

	}

	NTT(a,len,-1);

	for(int i = 0;i <= n + m;++i)printf("%lld ",a[i]);

	return 0;

}

多项式求逆

给出 $n-1$ 次多项式 $F(x)$ ，求一个 $n-1$ 次多项式 $G(x)$ 使得 $F(x)\times G(x)\equiv 1(mod\ x^n)$

若 $n=1$ ，则 $\%x=1$ ，那么要求的就是 $F[0]$ 的逆元，即 $F[0]^{P-2}$

若 $n>1$ ，使用倍增方法。

假设已经求出了 $F(x)G'(x)\equiv1\ (mod\ x^{\lceil \frac n 2\rceil})$

由于 $F(x)G(x)\equiv1(mod\ x^n)$ ，则 $F(x)G(x)\equiv 1\ (mod\ x^{\lceil \frac n 2 \rceil})$

两式相减得： $F(x)(G'(x)-G(x))\equiv 0\ (mod\ x^{\lceil \frac n 2 \rceil})$

由于 $F(x)$ 不为 $0$ ，则 $G'(x)-G(x)\equiv 0\ (mod\ x^{\lceil\frac n 2\rceil})$

平方，得： $G'(x)^2+G(x)^2-2G'(x)G(x)\equiv0\ (mod\ x^n)$

为什么模 $x^n$ 呢？因为根据卷积的定义，有

$\begin{align}(G'(x)-G(x))^2[i]=\sum_{j=0}^i((G'(x)-G(x))[i]\times(G'(x)-G(x))[j-i])\end{align}$

由于 $i$ 和 $j-i$ 至少有一项小于 $\lceil\frac n 2\rceil$ ，而 $(x^{\lceil\frac n 2\rceil})^2 \ge x^n$ ，所以乘完之后为零。

由于 $F(x)G(x)\equiv1\ (mod\ x^n)$ ，所以乘上一个 $F(x)$ 得： $F(x)G'(x)^2+G(x)-2G'(x)\equiv0\ (mod\ x^n)$

于是 $G(x)\equiv2G'(x)-F(x)G'(x)^2\ (mod\ x^n)$

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cstring>

using namespace std;

int n;

#define MAXN 200010

typedef long long ll;

const ll P = 998244353;

ll w[MAXN << 2],*g = w + (MAXN << 1);

ll A[MAXN << 1],B[MAXN << 1],c[MAXN << 1];

int rev[MAXN << 1];

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % P;

		a = a * a % P;

		b = b >> 1;

	}

	return res;

}

void NTT(ll f[],int l,int type)

{

	for(int i = 0;i < l;++i)

		if(i < rev[i])

			swap(f[i],f[rev[i]]);

	for(int i = 2;i <= l;i = i << 1)

	{

		ll wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			ll w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				ll u = f[k],t = w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		ll ni = power(l,P - 2);

		for(int i = 0;i < l;++i)

			f[i] = f[i] * ni % P;

	}

	return;

}

void work(int deg,ll *a,ll *b)

{

	if(deg == 1)

	{

		b[0] = power(a[0],P - 2);

		return;

	}

	work(((deg + 1) >> 1),a,b);

	int len,l;

	for(len = 1,l = 0;len < (deg << 1);len = len << 1,++l);

	for(int i = 1;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	for(int i = 0;i < deg;++i)c[i] = a[i];

	for(int i = deg;i < len;++i)c[i] = 0;

	NTT(b,len,1);NTT(c,len,1);

	for(int i = 0;i < len;++i)

	{

		b[i] = ((2 * b[i] % P) - (c[i] * b[i] % P * b[i] % P) + P) % P;

	}

	NTT(b,len,-1);

	for(int i = deg;i < len;++i)b[i] = 0;

	return;

}

int main()

{

	scanf("%d",&n);

	for(int i = 0;i < n;++i)scanf("%lld",&A[i]);

	work(n,A,B);

	for(int i = 0;i < n;++i)printf("%lld ",B[i]);

	return 0;

}

分治FFT

大概要求的就是给定长度为 $n-1$ 的数组 $g[1],g[2],\dots,g[n-1]$ ，求 $f[0],f[1],\dots,f[n-1]$ ，其中 $$ f[i]=\sum_{j=1}^if[i-j]g[j] $$ 边界为 $f[0]=1$ 。答案模 $998244353$ 。（洛谷分治FFT模板）

首先可以类似 $cdq$ 分治那样的做法，左边的位置对右边的贡献是独立的，假如求出了 $[l,mid]$ 的答案，对右边 $x$ 位置的贡献为： $$ w_x=\sum_{i=l}^{mid}f[i]\times g[x-i] $$ 这本质上是个卷积，可以用 $FFT$ 加速，时间复杂度 $O(n\log^2n )$ 。

但是还有更优秀的做法。

考虑多项式 $F$ ，把它写成形式幂级数或者说是生成函数的形式，即： $$ F(x)=\sum_{i=0}^\infty f[i]x^i $$ 由于 $g_0=0$ ，所以有 $F(x)G(x)=F(x)-f_0$ ，这个也很好理解，因为卷积的过程就是上面那个式子求的过程，所以自然 $F$ 除了最后一位别的都是不变的。

那么就有 $F(x)\equiv(1-G(x))^{-1}(\mod x^n)$ ，于是用多项式求逆解决就好了。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n;

#define MAXN 400010

#define P 998244353

typedef long long ll;

ll F[MAXN],G[MAXN];

ll ww[MAXN << 1],*g = ww + MAXN;

int rev[MAXN];

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % P;

		a = a * a % P;

		b = b >> 1;

	}

	return res;

}

void NTT(ll f[],int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		ll wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			ll w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				ll u = f[k],t = w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		for(int i = 0;i < l;++i)f[i] = f[i] * power(l,P - 2) % P;

	}

	return;

}

ll c[MAXN];

void calc_inv(int deg,ll a[],ll b[])

{

	if(deg == 1)

	{

		b[0] = power(a[0],P - 2) % P;

		return;

	}

	calc_inv((deg + 1) >> 1,a,b);

	int l = 0,len = 1;

	for(;len < (deg << 1);len = len << 1,++l);

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	for(int i = 0;i < deg;++i)c[i] = a[i];

	for(int i = deg;i < len;++i)c[i] = 0;

	NTT(c,len,1);NTT(b,len,1);

	for(int i = 0;i < len;++i)

	{

		b[i] = (2 * b[i] % P - c[i] % P * b[i] % P * b[i] % P + P) % P;

	}

	NTT(b,len,-1);

	for(int i = deg;i < len;++i)b[i] = 0;

	return;

}

int main()

{

	scanf("%d",&n);

	for(int i = 1;i < n;++i)scanf("%lld",&G[i]);

	for(int i = 1;i < n;++i)G[i] = P - G[i];

	G[0] = (G[0] + 1) % P;

	int len;

	for(len = 1;len < n;len = len << 1);

	calc_inv(len,G,F);

	for(int i = 0;i < n;++i)printf("%lld ",F[i]);

	cout << endl;

	return 0;

}

另一种做法是 $CDQ$ 分治，也就是考虑左边对右边的影响，用 $FFT$ 加速，注意下标的变化。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n;

#define MAXN 400010

typedef long long ll;

ll g[MAXN];

ll f[MAXN];

#define P 998244353

ll power(ll a,ll b)

{

    ll res = 1;

    while(b > 0)

    {

        if(b & 1)res = res * a % P;

        a = a * a % P;

        b = b >> 1;

    }

    return res;

}

namespace poly

{

    int rev[MAXN];

    ll ww[MAXN << 1],*g = ww + MAXN;

    void NTT(ll f[],int l,int type)

    {

        for(int i = 0;i < l;++i)

        {

            if(i < rev[i])swap(f[i],f[rev[i]]);

        }

        for(int i = 2;i <= l;i = i << 1)

        {

            ll wn = g[type * l / i];

            for(int j = 0;j < l;j += i)

            {

                ll w = 1;

                for(int k = j;k < j + i / 2;++k)

                {

                    ll u = f[k],t = w * f[k + i / 2] % P;

                    f[k] = (u + t) % P;

                    f[k + i / 2] = (u - t + P) % P;

                    w = w * wn % P;

                }

            }

        }

        if(type == -1)

        {

            ll ni = power(l,P - 2);

            for(int i = 0;i < l;++i)f[i] = f[i] * ni % P;

        }

        return;

    }

    void mul(ll a[],ll b[],ll res[],int deg)

    {

        int l = 0,len = 1;

        for(;len <= (deg << 1);len = len << 1)++l;

        g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

        for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

        for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

        NTT(a,len,1);NTT(b,len,1);

        for(int i = 0;i < len;++i)res[i] = a[i] * b[i] % P;

        NTT(res,len,-1);

        for(int i = 0;i < len;++i)a[i] = b[i] = 0;

        return;

    }

}

ll tmp1[MAXN],tmp2[MAXN],tmp3[MAXN];

void cdq(int l,int r)

{

    if(l == r)return;

    int mid = ((l + r) >> 1);

    cdq(l,mid);

    for(int i = l;i <= mid;++i)tmp1[i - l] = f[i];

    for(int i = 0;i <= r - l + 1;++i)tmp2[i] = g[i];

    int d = r - l + 1 + 1;

    poly::mul(tmp1,tmp2,tmp3,d);

    for(int i = mid + 1;i <= r;++i)f[i] = (f[i] + tmp3[i - l]) % P;

    for(int i = 0;i < d;++i)tmp1[i] = tmp2[i] = tmp3[i] = 0;

    cdq(mid + 1,r);

    return;

}

int main()

{

    scanf("%d",&n);

    for(int i = 1;i < n;++i)scanf("%lld",&g[i]);

    f[0] = 1;

    cdq(0,n);

    for(int i = 0;i < n;++i)printf("%lld ",f[i]);

    cout << endl;

    return 0;

}

多项式除法

给定一个 $n$ 次多项式 $F(x)$ 和一个 $m$ 次多项式 $G(x)$ ，请求出多项式 $Q(x)$ ， $R(x)$ ，满足以下条件：

$Q(x)$ 次数为 $n-m$ ， $R(x)$ 次数小于 $m$ ， $F(x) = Q(x) \times G(x) + R(x)$

所有的运算在模 $998244353$ 意义下进行。

首先定义多项式的反转，设一个次数界为 $n$ 的多项式 $F(x)$ ， $F^R(x)=x^nF(\frac 1 x)$ ，因为观察一下就会发现这个实际上是把原多项式高低次数的位置依次交换。

那么我们可以进行这样一个操作，把上面那个式子两边同成 $x^n$ ，并把 $\frac 1 x$ 当作 $x$ 带入，即： $$ \begin{align} &x^nF(\frac 1 x)=x^{n-m}Q(\frac 1 x)\times x^mG(\frac 1 x)+x^{n-m+1}x^{m-1}R(\frac 1 x)\\ \Longrightarrow&F^R(x)=Q^R(x)\times G^R(x)+x^{n-m+1}R^R(x) \end{align} $$ 由于我们最终要求的 $Q(x)$ 是一个 $n-m$ 次多项式，所以有： $$ F^R(x)=Q^R(x)\times G^R(x)(\mod x^{n-m+1}) $$ 所以 $Q^R=F^R(x)\times (G^R)^{-1}(\mod x^{n-m+1})$ ，多项式求逆就可以了。

多项式取模

多项式除法后把 $Q$ 带回去解出 $R$ 即可。

除了优化递推好像没啥用。

多项式除法和取模的代码：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,m;

#define MAXN 300010

typedef long long ll;

const ll P = 998244353;

ll ww[MAXN << 1],*g = ww + MAXN;

int rev[MAXN];

ll F[MAXN],G[MAXN],Q[MAXN];

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % P;

		a = a * a % P;

		b = b >> 1;

	}

	return res;

}

void NTT(ll f[],int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		ll wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			ll w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				ll u = f[k],t = w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		ll inv = power(l,P - 2);

		for(int i = 0;i < l;++i)f[i] = f[i] * inv % P;

	}

	return;

}

ll c[MAXN];

void calc_inv(int deg,ll a[],ll b[])

{

	if(deg == 1)

	{

		b[0] = power(a[0],P - 2);

		return;

	}

	calc_inv((deg + 1) >> 1,a,b);

	int l = 0,len = 1;

	for(;len < (deg << 1);len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	for(int i = 0;i < deg;++i)c[i] = a[i];

	for(int i = deg;i < len;++i)c[i] = 0;

	NTT(c,len,1);NTT(b,len,1);

	for(int i = 0;i < len;++i)

	{

		b[i] = (2 * b[i] % P - c[i] * b[i] % P * b[i] % P + P) % P;

	}

	NTT(b,len,-1);

	for(int i = deg;i < len;++i)b[i] = 0;

	return;

}

int main()

{

	scanf("%d%d",&n,&m);

	for(int i = 0;i <= n;++i)scanf("%lld",&F[n - i]);

	for(int i = 0;i <= m;++i)scanf("%lld",&G[m - i]);

	calc_inv(n - m + 1,G,Q);

	++n;++m;

	int l = 0,len = 1;

	for(;len < (n << 1);len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	NTT(F,len,1);NTT(Q,len,1);

	for(int i = 0;i < len;++i)Q[i] = Q[i] * F[i] % P;

	NTT(Q,len,-1);NTT(F,len,-1);

	for(int i = n - m + 1;i < len;++i)Q[i] = 0;

	for(int i = 0,j = n - m;i < j;++i,--j)swap(Q[i],Q[j]);

	for(int i = 0;i < n - m + 1;++i)printf("%lld ",Q[i]);puts("");

	for(int i = 0,j = n - 1;i < j;++i,--j)swap(F[i],F[j]);

	for(int i = 0,j = m - 1;i < j;++i,--j)swap(G[i],G[j]);

	l = 0,len = 1;

	for(;len < n + 1;len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = g[i - 1] * g[1] % P;

	NTT(G,len,1);NTT(Q,len,1);

	for(int i = 0;i < len;++i)Q[i] = Q[i] * G[i] % P;

	NTT(Q,len,-1);

	for(int i = 0;i < len;++i)F[i] = (F[i] - Q[i] + P) % P;

	for(int i = 0;i < m - 1;++i)printf("%d ",F[i]);puts("");

	return 0;

}

##多项式牛顿迭代法

假设我们有一个多项式的函数 $G(x)$ ，我们想求它的一个零点，换一种方法来说就是求解： $$ G(F(x))=0(\mod x^n) $$ 首先当 $n=1$ 时， $G(F(x))=0(\mod x)$ ，这个需要单独求出。

当 $n\ne 1$ 时，假设已经求出了： $$ G(F_0(x))=0(\mod x^{\lfloor\frac n 2\rfloor}) $$ 我们可以把 $G(x)$ 在 $F_0(x)$ 这里泰勒展开，即： $$ G(F(x))=G(F_0(x))+\frac{G'(F_0(x))}{1!}(F(x)-F_0(x))+\frac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+\cdots $$ 由于 $F(x)$ 和 $F_0(x)$ 最后 $\lfloor\frac n 2\rfloor$ 项一定相同，因为他们实际上只是模的东西不一样，本质是一个东西，所以 $(F(x)-F_0(x))^2$ 后 $n$ 项一定都是 $0$ ，于是就可以模掉了，于是有： $$ 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)-\frac{G(F_0(x))}{G'(F_0(x))} $$ 所以如果支持把 $F_0(x)$ 快速代入 $G()$ ，那么我们就有了一个复杂度为 $T(n)=T(\frac n 2)+O(代入G+n\log n)$ 的求解多项式函数 $G(x)$ 的零点的做法。

多项式开方

给出多项式 $H(x)$ ，要求一个多项式 $F(x)$ 满足： $$ F^2(x)=H(x)(\mod x^n) $$ 做法还是利用牛顿迭代，设函数 $G(x)$ 为： $$ G(F(x))=F^2(x)-H(x) $$ 带入牛顿迭代公式为： $$ F=F_0-\frac{G(F_0)}{G'(F_0)}=F_0-\frac{F_0^2-H}{2F_0} $$ 时间复杂度 $O(n\log n)$ 。

注意当 $deg=0$ 的时候要求模意义下的二次剩余，二次剩余有一个叫做 $Cipolla$ 的算法，但是如果 $0$ 次项的值已知那么就可以在下面把二次剩余枚举得出。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,m;

#define MAXN 400010

int ww[MAXN << 1],*g = ww + MAXN;

int rev[MAXN];

#define P 998244353

int power(int a,int b)

{

	int res = 1;

	while(b > 0)

	{

		if(b & 1)res = 1ll * res * a % P;

		a = 1ll * a * a % P;

		b = b >> 1;

	}

	return res;

}

void NTT(int f[],int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		int wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			int w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				int u = f[k],t = 1ll * w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = 1ll * w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		int ni = power(l,P - 2);

		for(int i = 0;i < l;++i)f[i] = 1ll * f[i] * ni % P;

	}

	return;

}

int invtmp[MAXN];

void calc_inv(int deg,int a[],int b[])

{

	if(deg == 1)

	{

		b[0] = power(a[0],P - 2);

		return;

	}

	calc_inv((deg + 1) >> 1,a,b);

	int l = 0,len = 1;

	for(;len <= (deg << 1);len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = 1ll * g[i - 1] * g[1] % P;

	for(int i = 0;i < deg;++i)invtmp[i] = a[i];

	for(int i = deg;i < len;++i)invtmp[i] = 0;

	NTT(invtmp,len,1);NTT(b,len,1);

	for(int i = 0;i < len;++i)

	{

		b[i] = (2 * b[i] % P - 1ll * invtmp[i] * b[i] % P * b[i] % P + P) % P;

	}

	NTT(b,len,-1);

	for(int i = deg;i < len;++i)b[i] = 0;

	//for(int i = 0;i < deg;++i)cout << b[i] * 2 % P << " ";cout << endl;

	return;

}

int multmp[MAXN];

void mul(int deg,int a[],int b[])

{

	int l = 0,len = 1;

	for(;len <= (deg << 1);len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = 1ll * g[i - 1] * g[1] % P;

	for(int i = 0;i < len;++i)multmp[i] = b[i];

	NTT(a,len,1);NTT(multmp,len,1);

	for(int i = 0;i < len;++i)

	{

		a[i] = 1ll * a[i] * multmp[i] % P;

	}

	NTT(a,len,-1);

	for(int i = deg;i < len;++i)a[i] = 0;

	return;

}

#define REM 1

int tmp1[MAXN],tmp2[MAXN];

void calc_sqr(int deg,int a[],int b[])

{

	if(deg == 1)

	{

		b[0] = REM;

		return;

	}

	calc_sqr((deg + 1) >> 1,a,b);

	for(int i = 0;i < deg;++i)tmp1[i] = b[i];

	mul(deg,tmp1,b);

	for(int i = 0;i < deg;++i)tmp1[i] = (tmp1[i] - a[i] + P) % P;

	memset(tmp2,0,sizeof(tmp2));

	calc_inv(deg,b,tmp2);

	for(int i = 0;i < deg;++i)tmp2[i] = 1ll * tmp2[i] * power(2,P - 2) % P;

	mul(deg,tmp1,tmp2);

	for(int i = 0;i < deg;++i)b[i] = (b[i] - tmp1[i] + P) % P;

	for(int i = 0;i < deg;++i)cout << b[i] << " ";cout << endl;

	return;

}

int H[MAXN],F[MAXN];

int main()

{

	scanf("%d",&n);

	for(int i = 0;i < n;++i)scanf("%d",&H[i]);

	calc_sqr(n,H,F);

	for(int i = 0;i < n;++i)printf("%d ",F[i]);

	cout << endl;

	return 0;

}

多项式取模优化常系数齐次线性递推

常系数齐次线性递推，就是给你一个线性递推关系： $$ f_n=\sum_{i=1}^ka_if_{n-i} $$ 让你求第 $n$ 项的值。

显然可以矩阵乘法，但是矩阵乘法复杂度是 $O(k^3\log n)$ ，不难发现复杂度的瓶颈主要在 $k^3$ 这部分，因为对于绝大多数问题， $\log n$ 都比较小。

特征多项式：

特征值，特征向量：

如果有一个常数 $\lambda$ 和一个向量 $\overrightarrow v$ ，满足 $\lambda\overrightarrow v=A\overrightarrow v$ ，则称 $\overrightarrow v$ 为矩阵 $A$ 的一组特征向量， $\lambda$ 为矩阵 $A$ 的一组特征值。

秩为 $k$ 的矩阵有 $k$ 组线性不相关的特征向量。

特征多项式，特征方程：

对关系式进行变换： $(\lambda I-A)\overrightarrow v=0$ ，这个等式有解的充要条件是 $\det(\lambda I-A)=0$ ，那么我们可以把 $\det(\lambda I-A)$ 看作一个 $k$ 次多项式，记为： $f(\lambda)=\det(\lambda I-A)$ ，这个称为矩阵的特征多项式， $f(\lambda)$ 称为矩阵的特征方程，这个方程有 $k$ 个解可以写作： $$ f(\lambda)=\prod_{i=1}^k(\lambda_i-\lambda) $$ 每个解是一个特征值。

Cayley-Hamilton定理：

根据：Cayley-Hamilton定理，有 $f(A)=0$ ，证明的话： $$ \begin{align} &f(A)=\prod_{i=1}^k(\lambda_iI-A)\\ &首先有：(\lambda_iI-A)\times(\lambda_j I-A)=\lambda_i\lambda_jI^2-(\lambda_i+\lambda_j)IA+A^2=(\lambda_jI-A)\times(\lambda_i I-A)\\ &我们只要证明f(A)乘任何特征向量都会得到0矩阵，就可以说明矩阵A=0\\ &f(A)\times \overrightarrow {v_i}=\prod_{j\ne i}(\lambda_jI-A)(\lambda_i I-A)\overrightarrow {v_i}=\prod_{j\ne i}(\lambda_jI-A)\times 0=0\\ &\therefore f(A)=0 \end{align} $$

优化递推：

设转移矩阵为 $M$ ，根据矩阵快速幂，我们只要能求出来 $M^n$ 就行了。

$x^n$ 可以看作一个 $n$ 次多项式，那么我们用 $f(x)$ 来表示就是 $x^n=f(x)\times g(x)+h(x)$ 。

即 $M^n=f(M)\times g(M)+h(M)$ ，但是由于 $f(M)=0$ ，因此： $M^n=h(M)$ ，于是我们可以倍增：也就是每次 $H(n)=h(n)\times h(n)$ ，但是 $H(n)$ 是 $2k-2$ 次的，因此我们可以手动多项式取模，复杂度 $O(k^2\log n)$ ，如果用 $FFT$ 卷积 $+$ 多项式除法可以做到 $O(k\log k\log n)$ 。

还有一个问题，就是 $f(x)$ 怎么求，手动模拟一下求行列式的过程会发现 $[x^i]f(x)=-a[k-i],[x^k]f(x)=1$ 。

总之当个板子记住就行了。

Code：

暴力取模：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,k;

#define MAXK 4010

#define MOD 1000000007

inline int rd()

{

	register int res = 0,f = 1;register char c = getchar();

	while(!isdigit(c)){if(c == '-')f = -1;c = getchar();}

	while(isdigit(c))res = (res << 1) + (res << 3) + c - '0',c = getchar();

	return res * f;

}

int f[MAXK],h[MAXK];

struct poly

{

	int a[MAXK];

	poly(){memset(a,0,sizeof(a));}

	int& operator [](const int &i){return a[i];}

	int operator [](const int &i)const{return a[i];}

	inline poly operator * (const poly &b)

	{

		poly res;

		for(register int i = 0;i < k;++i)

			for(register int j = 0;j < k;++j)

				(res[i + j] += 1ll * a[i] * b[j] % MOD) %= MOD;

		for(register int i = 2 * k - 2;i >= k;res[i--] = 0)

			for(register int j = 1;j <= k;++j)

				(res[i - j] += 1ll * res[i] * f[j] % MOD) %= MOD;

		return res;

	}

}m;

poly power(poly a,int b)

{

	register poly res;

	res[0] = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a;

		a = a * a;

		b = b >> 1;

	}

	return res;

}

int main()

{

	scanf("%d%d",&n,&k);

	for(register int i = 1;i <= k;++i)f[i] = rd(),f[i] = (f[i] > 0 ? f[i] : f[i] + MOD);

	for(register int i = 0;i < k;++i)h[i] = rd(),h[i] = (h[i] > 0 ? h[i] : h[i] + MOD);

	if(n < k){cout << h[n] << endl;return 0;}

	m[1] = 1;int ans = 0;

	m = power(m,n);

	for(register int i = 0;i < k;++i)(ans += 1ll * m[i] * h[i] % MOD) %= MOD;

	cout << ans << endl;

	return 0;

}

$FFT$ 卷积 $+$ 多项式取模：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,k;

#define MAXK 248000

#define MOD 998244353

inline int rd()

{

	register int res = 0,f = 1;register char c = getchar();

	while(!isdigit(c)){if(c == '-')f = -1;c = getchar();}

	while(isdigit(c))res = (res << 1) + (res << 3) + c - '0',c = getchar();

	return res * f;

}

int f_[MAXK],h[MAXK];

namespace polynomial

{ 

	int rev[MAXK];

	int ww[MAXK << 1],*g = ww + MAXK;

	int power(int a,int b)

	{

		int res = 1;

		while(b > 0)

		{

			if(b & 1)res = 1ll * res * a % MOD;

			a = 1ll * a * a % MOD;

			b = b >> 1;

		}

		return res;

	}

	int init(int n)

	{

		int l = 0,len = 1;

		for(;len <= n;len = len << 1)++l;

		for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

		g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(MOD - 1) / len);

		for(int i = 2;i < len;++i)g[i] = g[i - len] = 1ll * g[i - 1] * g[1] % MOD;

		return len;

	}

	struct poly

	{

		int a[MAXK];

		int d;

		poly(){memset(a,0,sizeof(a));d = 0;}

		int& operator [](int x){return a[x];}

		int operator [](int x)const{return a[x];}

	};

	void NTT(poly &f,int l,int type)

	{

		for(int i = 0;i < l;++i)

		{

			if(i < rev[i])swap(f[i],f[rev[i]]);

		}

		for(int i = 2;i <= l;i = i << 1)

		{

			int wn = g[type * l / i];

			for(int j = 0;j < l;j += i)

			{

				int w = 1;

				for(int k = j;k < j + i / 2;++k)

				{

					int u = f[k],t = 1ll * w * f[k + i / 2] % MOD;

					f[k] = (u + t) % MOD;

					f[k + i / 2] = (u - t + MOD) % MOD;

					w = 1ll * w * wn % MOD;

				}

			}

		}

		if(type == -1)

		{

			int ni = power(l,MOD - 2);

			for(int i = 0;i < l;++i)f[i] = 1ll * f[i] * ni % MOD;

		}

		return;

	}

	poly operator * (poly a,poly b)

	{

		poly res;

		int len = init(a.d + b.d);res.d = a.d + b.d;

		NTT(a,len,1);NTT(b,len,1);

		for(int i = 0;i < len;++i)res[i] = 1ll * a[i] * b[i] % MOD;

		NTT(res,len,-1);

		return res;

	}

	poly invtmp;

	void calc_inv(poly &a,poly &b,int deg)

	{

		if(deg == 1){b[0] = power(a[0],MOD - 2);return;}

		calc_inv(a,b,((deg + 1) >> 1));

		int len = init(deg << 1);

		for(int i = 0;i < deg;++i)invtmp[i] = a[i];

		for(int i = deg;i < len;++i)invtmp[i] = 0;

		NTT(invtmp,len,1);NTT(b,len,1);

		for(int i = 0;i < len;++i)

		{

			b[i] = (2ll * b[i] % MOD - 1ll * invtmp[i] * b[i] % MOD * b[i] % MOD + MOD) % MOD;

		}

		NTT(b,len,-1);

		for(int i = deg;i < len;++i)b[i] = 0;

		b.d = deg - 1;

		return;

	}

	poly operator + (poly a,poly b)

	{

		poly res;

		res.d = max(a.d,b.d);

		for(int i = 0;i < res.d;++i)res[i] = (a[i] + b[i]) % MOD;

		return res;

	}

	poly operator - (poly a,poly b)

	{

		poly res;

		res.d = max(a.d,b.d);

		for(int i = 0;i < res.d;++i)res[i] = (a[i] - b[i] + MOD) % MOD;

		return res;

	}

	poly divtmp;

	poly operator / (poly a,poly b)

	{

		for(int i = 0;i <= divtmp.d;++i)divtmp[i] = 0;

		for(int i = 0,j = a.d;i < j;++i,--j)swap(a[i],a[j]);

		for(int i = 0,j = b.d;i < j;++i,--j)swap(b[i],b[j]);

		calc_inv(b,divtmp,b.d + 1);

		divtmp = divtmp * a;

		for(int i = a.d - b.d + 1;i <= divtmp.d;++i)divtmp[i] = 0;

		divtmp.d = a.d - b.d;

		for(int i = 0,j = divtmp.d;i < j;++i,--j)swap(divtmp[i],divtmp[j]);

		return divtmp;

	}

	poly operator % (poly a,poly b)

	{

		poly res = a - a / b * b;

		res.d = b.d - 1;

		return res;

	}

}

using namespace polynomial;

poly f;

struct matrix

{

	poly a;

	int& operator [](int x){return a[x];}

	int operator [](int x)const{return a[x];}

	friend inline matrix operator * (matrix a,matrix b)

	{

		matrix res;

		res.a = a.a * b.a;

		res.a = res.a % f;

		return res;

	}

}m;

matrix power(matrix a,int b)

{

	matrix res;

	res[0] = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a;

		a = a * a;

		b = b >> 1;

	}

	return res;

}

int main()

{

	scanf("%d%d",&n,&k);

	for(register int i = 1;i <= k;++i)f_[i] = rd(),f_[i] = (f_[i] > 0 ? f_[i] : f_[i] + MOD);

	for(register int i = 0;i < k;++i)h[i] = rd(),h[i] = (h[i] > 0 ? h[i] : h[i] + MOD);

	if(n < k){cout << h[n] << endl;return 0;}

	f[k] = 1;

	for(int i = 0;i < k;++i)f[i] = MOD - f_[k - i];

	f.d = k;

	m[1] = 1;m.a.d = k - 1;int ans = 0;

	m = power(m,n);

	ans = 0;

	for(int i = 0;i < k;++i)(ans += 1ll * m[i] * h[i] % MOD) %= MOD;

	cout << ans << endl;

	return 0;

}

多项式多点求值

多项式快速插值

多项式辗转相除法

多项式求导和积分

求导： $$ (x^a)'=ax^{a-1} $$

for(int i = 1;i < n;++i)f[i - 1] = 1ll * i * f[i] % P;

f[n - 1] = 0;

不定积分： $$ \int^x_0 t^adt=\frac1{a+1}x^{a+1} $$

for(int i = n - 1;i >= 1;--i)f[i] = 1ll * f[i - 1] * power(i,P - 2) % P;

多项式对数函数

求 $g(x)=\ln f(x)$ 。

对两边同时求导，得： $$ g'(x)=(\ln f(x))'=\ln'f(x)f'(x)=\frac{f'(x)}{f(x)} $$ 所以就先对 $f$ 求导除以 $f$ 的逆然后再不定积分就可以了。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n;

#define MAXN 400010

int f[MAXN],f_[MAXN];

#define P 998244353

int ww[MAXN << 1],*g = ww + MAXN;

int rev[MAXN];

int a[MAXN],b[MAXN],c[MAXN];

int power(int a,int b)

{

	int res = 1;

	while(b > 0)

	{

		if(b & 1)res = 1ll * res * a % P;

		a = 1ll * a * a % P;

		b = b >> 1;

	}

	return res;

}

void NTT(int f[],int l,int type)

{

	for(int i = 0;i < l;++i)

	{

		if(i < rev[i])swap(f[i],f[rev[i]]);

	}

	for(int i = 2;i <= l;i = i << 1)

	{

		int wn = g[type * l / i];

		for(int j = 0;j < l;j += i)

		{

			int w = 1;

			for(int k = j;k < j + i / 2;++k)

			{

				int u = f[k],t = 1ll * w * f[k + i / 2] % P;

				f[k] = (u + t) % P;

				f[k + i / 2] = (u - t + P) % P;

				w = 1ll * w * wn % P;

			}

		}

	}

	if(type == -1)

	{

		int ni = power(l,P - 2);

		for(int i = 0;i < l;++i)f[i] = 1ll * f[i] * ni % P;

	}

	return;

}

void calc_inv(int deg,int a[],int b[])

{

	if(deg == 1)

	{

		b[0] = power(a[0],P - 2);

		return;

	}

	calc_inv((deg + 1) >> 1,a,b);

	int l = 0,len = 1;

	for(;len < (deg << 1);len = len << 1,++l);

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i <= len;++i)g[i] = g[i - len] = 1ll * g[i - 1] * g[1] % P;

	for(int i = 0;i < deg;++i)c[i] = a[i];

	for(int i = deg;i < len;++i)c[i] = 0;

	NTT(c,len,1);NTT(b,len,1);

	for(int i = 0;i < len;++i)

	{

		b[i] = (2 * b[i] % P - 1ll * c[i] * b[i] % P * b[i] % P + P) % P;

	}

	NTT(b,len,-1);

	for(int i = deg;i < len;++i)b[i] = 0;

	return;

}

int main()

{

	scanf("%d",&n);

	for(int i = 0;i < n;++i)scanf("%d",&f[i]);

	calc_inv(n,f,f_);

	//for(int i = 0;i < n;++i)cout << f_[i] << " ";cout << endl;

	for(int i = 1;i < n;++i)f[i - 1] = 1ll * i * f[i] % P;

	//for(int i = 0;i < n;++i)cout << f[i] << " ";cout << endl;

	f[n - 1] = 0;

	int l = 0,len = 1;

	for(;len <= (n << 1);len = len << 1)++l;

	for(int i = 0;i < len;++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

	g[0] = g[-len] = 1;g[1] = g[1 - len] = power(3,(P - 1) / len);

	for(int i = 2;i < len;++i)g[i] = g[i - len] = 1ll * g[i - 1] * g[1] % P;

	NTT(f,len,1);NTT(f_,len,1);

	for(int i = 0;i < len;++i)f[i] = 1ll * f[i] * f_[i] % P;

	NTT(f,len,-1);

	for(int i = n - 1;i >= 1;--i)f[i] = 1ll * f[i - 1] * power(i,P - 2) % P;

	f[0] = 0;

	for(int i = 0;i < n;++i)printf("%d ",f[i]);cout << endl;

	return 0;

}