多项式小结

Preface

• 蒟蒻的我最近才会多项式exp……
• 在这里总结一下多项式的各种芝士。顺便挂个模板，以备未来之需。

万恶之源——FFT/NTT

• 我一个图像分析算法，怎么跑到OI来了
• 一次dft，要求解\(A(x)\)在\(\omega_n^i(i\in[0,n))\)的点值；idft要根据这些点值插出一个\(n-1\)次多项式。
• FFT核心思想：记\(A(x)=\sum_i a_ix^i,A^0(x)=\sum_i a_{2i}x^i,A^1(x)=\sum_i a_{2i+1}x^i\)。则有：

\[A(\omega_n^i)=A^0(\omega_n^{2i})+A^1(\omega_n^{2i})*\omega_n^i \]

\[A(\omega_n^{i+\frac n2})=A^0(\omega_n^{2i})+A^1(\omega_n^{2i})*\omega_n^{i+\frac n2}=A^0(\omega_n^{2i})-A^1(\omega_n^{2i})*\omega_n^i \]

• 普通FFT的时候，\(\omega_n^j=e^{2\pi ij/n}=\cos\frac{2\pi j}n+i\sin\frac{2\pi j}n\)；NTT的时候，\(\omega_n^i=g^{i(mo-1)/n}\)。

求逆

• 求\(B(x)\equiv A^{-1}(x)(mod\ x^n)\)。
• 已知\(B_0(x)\equiv A^{-1}(x)(mod\ x^{n/2})\)，则\(B(x)-B_0(x)\equiv 0(mod\ x^{n/2})\)。
• 两边平方、移项，再乘上个\(A(x)\)，得到\(B(x)\equiv 2B_0(x)-A(x)B_0^2(x)(mod\ x^n)\)。

除法、取模

• 已知\(n\)次\(A(x)\)，\(m(m≤n)\)次\(B(x)\)，求\(D(x),R(x)\)，使\(A(x)=D(x)B(x)+R(x)\)。
• \(x^nA(\frac 1x)=x^{n-m}D(\frac 1x)x^mB(\frac 1x)+x^{n-m+1}x^{m-1}R(\frac 1x)\)。
• 记\(A^R(x)=x^nA(\frac 1x)\)，可以发现这个操作就是将其系数reverse。把\(R(x)\)看成\(m-1\)次（不足补0），则\(A^R(x)=D^R(x)B^R(x)+x^{n-m+1}R^R(x)\)。
• \(A^R(x)\equiv D^R(x)B^R(x)(mod\ x^{n-m+1})\)。

开根

• 已知\(B_0^2\equiv A(mod\ x^{n/2})\)，求\(B^2\equiv A(mod\ x^n)\)。
• 移项、平方，得到\((B_0-B)^2\equiv 0(mod\ x^n)\)。
• 拆开来、移项，得到\(B_0^2+B^2\equiv 2B_0B(mod\ x^n)\)。
• \(B\equiv B_0/2+A/(2B_0)(mod\ x^n)\)。
• 注意如果常数项不为1，可能要打个二次剩余。（坑点：\(x^2\equiv n(mod\ 998244353)\)，同时\((998244353-x)^2\equiv n(mod\ 998244353)\))

ln

• \(G(x)=\ln F(x)=\int\frac{F'(x)}{F(x)}dx\)。
• 注意这个当\(A(x)\)的常数项不为1时是求不了的。如果想求多项式幂函数，要先提出一个因式，使得\(A(x)\)的常数项为1。

exp

• 求\(B(x)=e^{A(x)}(mod\ x^n)\)。
• 设\(F(B)=\ln B-A\)，我们要求它\(mod\ x^n\)意义下的零点。（注意，这里是以\(B\)为自变量，\(A\)对于\(F(B)\)来说是个常数，因此\(F'(B)=\frac 1B\)）
• 已知\(F(B_0)\equiv 0(mod\ x^{n/2})\)。泰勒展开：\(F(B)=F(B_0)+\sum_i \frac{F^{(i)}(B_0)(B-B_0)^i}{i!}\)。
• \(B-B_0\equiv 0(mod\ x^{n/2})\)，二次方后\(mod\ x^n=0\)，因此可以不管\(i>1\)的项。于是得到\(F(B)\equiv F'(B_0)(B-B_0)(mod\ x^n)\)。
• 因此\(B\equiv B_0-\frac{F(B_0)}{F'(B_0)}\equiv B_0(1-\ln B_0+A)(mod\ x^n)\)。
• 注意这个当\(A(x)\)的常数项不为0时，应该是求不了的。

常系数齐次线性递推

• 求一个满足\(k\)阶齐次线性递推数列\(a_i\)的第\(n\)项，即：\(a_n=\sum_{i=1}^k f_i\times a_{n-i}\)。
• 看这篇题解秒懂：https://www.luogu.com.cn/blog/BJpers2/solution-p4723
• 不需要任何高级线代知识，立即可知我们只需求\(x^n\)对多项式\(f(x)=x^k-\sum_{i=0}^{k-1}f_{k-i}x^i\)的结果。
• 数学追求简洁性，“譬如从北京城里走到颐和园那样，可有许多条路，要选择一条最准确无错误，又最短最好的道路”；我觉得信息学里亦然。不需要矩阵特征值，行列式，C-H定理这些东西，就不必把它们搬出来。

三角函数

• 一个无聊的东西。
• 根据欧拉公式有：\(e^{ix}=\cos x+i\sin x\)，解方程可得\(\sin x=\frac{e^{ix}-e^{-ix}}{2i},\cos x=\frac{e^{ix}+e^{-ix}}2\)。
• 把多项式代入\(x\)即可。在模\(998244353\)意义下，\(i=\sqrt{998244352}=±86583718\)。

---

• 注意，以下皆为口胡……

---

反三角函数

• 一个更无聊的东西。
• 对于\(y=\arcsin x\)，有\(y'=\frac 1{\sqrt{1-x^2}}\)。那么对于\(G(x)=\arcsin F(x)\)，有\(G=\int\frac{F'}{\sqrt{1-F^2}}dx\)。这个exp都不用了，直接求逆+开根。
• 对于\(y=\arctan x\)，有\(y'=\frac 1{1+x^2}\)。那么对于\(G(x)=\arctan F(x)\)，有\(G=\int\frac {F'}{1+F^2}dx\)。这个开根都不用了。

多点求值

• 已知多项式\(A(x)\)和\(n\)个点\(x_1,x_2,...,x_n\)，求\(A(x_1),A(x_2),...,A(x_n)\)。
• 考虑分治，将\(n\)个点分为前后两半。记\(B_0(x)=\prod_{i=1}^{n/2}(x-x_i)\)，则有\(\forall i\in[1,n/2]:B_0(x_i)=0\)。
• 若我们将\(A(x)\)表示成\(D(x)B_0(x)+R(x)\)的形式，可得\(\forall i\in[1,n/2]:A(x)=R(x)\)。
• 那么就可以用分治NTT+多项式取模解决。另一半同理。

快速插值

• 普通的朗格朗日插值长这样：\(F(x)=\sum_{i=1}^n y_i\prod_{j≠i}\frac{x-x_j}{x_i-x_j}\)。

• 记\(M(x)=\prod_{i=1}^n(x-x_i)\)，洛一洛可得：\(\prod_{j\neq i} (x_i - x_j) = \lim_{x\rightarrow x_i} \frac{\prod_{j=1}^n (x - x_j)}{x - x_i} = M'(x_i)\)。

• 因此\(F(x)=\sum_{i = 1}^n \frac{y_i}{M'(x_i)}\prod_{j \neq i}(x - x_j)\)。

• 再考虑分治。令\(F_0(x)=\sum_{i=1}^{n/2}\frac{y_i}{M'(x_i)}\prod_{j \neq i\wedge j≤n/2}(x - x_j),M_0=\prod_{i=1}^{n/2}(x-x_i)\)（另一半同理），则可得\(F(x)=F_0(x)M_1(x)+F_1(x)M_0(x)\)。

Code

• 注意：主程序是常系数齐次线性递推。

#include<ctime>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define P(x,y) x=(x+y)%mo
#define T(x,y) x=x*(y)%mo
#define clr(a,n) memset(a,0,8*(n))
#define go(i,a,b) for(int i=a;i< b;i++)
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define CON(a,b,n) DFT(a,n,0);go(i,0,n)T(a[i],b[i]);DFT(a,n,1);
using namespace std;
typedef long long ll;
typedef unsigned long long ul;
const int N=1<<18,L=1<<17;
const ll mo=998244353,_1=86583718,_2=mo+1>>1;
int n,k,n1,n2;
ll a[N],b[N],p[N],t[N],s,ans;
ll ksm(ll x,ll y) {ll a=1; for(;y;y/=2,T(x,x))if(y&1)T(a,x); return a;}
ll iv(ll x) {return ksm(x,mo-2);}
ll R() {return (rand()*19260817ll+rand())%mo;}
namespace NTT
{
	int tf[N]; ll w[N];
	void build(int n)
	{
		for(int i=1;i<n;i*=2)
		{
			w[i]=1; ll v=ksm(114514,(mo-1)/2/i);
			go(j,1,i) w[i+j]=w[i+j-1]*v%mo;
		}
	}
	void DFT(ll a[N],int n,bool f)
	{
		go(i,1,n) if(i<(tf[i]=tf[i/2]/2+(i&1)*n/2)) swap(a[i],a[tf[i]]);
		ll v;
		for(int i=1;i<n;i*=2) for(int j=0;j<n;j+=2*i) go(k,0,i)
			v=a[i+j+k]*w[i+k]%mo, a[i+j+k]=(a[j+k]-v+mo)%mo, P(a[j+k],v);
		if(f)
		{
			reverse(a+1,a+n);
			v=mo-(mo-1)/n;
			go(i,0,n) T(a[i],v);
		}
	}
	void mtp(ll a[N],ll b[N],int n) {DFT(b,n,0); CON(a,b,n)}
	void inv(ll a[N],ll b[N],int m)
	{
		static ll c[N]; int n1;
		for(n1=1;n1<m;n1*=2);
		clr(b,n1); b[0]=iv(a[0]);
		for(int n=2;n<=n1;n*=2)
		{
			DFT(b,n*2,0);
			clr(c,n*2); memcpy(c,a,8*n); DFT(c,n*2,0);
			go(i,0,n*2) T(b[i],2-c[i]*b[i]%mo+mo);
			DFT(b,n*2,1); clr(b+n,n);
		}
		clr(b+m,n1-m);
	}
	void mul(ll a[N],ll b[N])
	{
		static ll c[N];
		memcpy(c,b,8*n1); DFT(c,n1,0);
		CON(a,c,n1)
		go(i,0,n1) c[i]=i>=n1-k?0:a[n1-i-1];
		CON(c,t,n2)
		reverse(c,c+n1-k); clr(c+n1-k,n1+k);
		CON(c,p,n1)
		go(i,0,n1) P(a[i],mo-c[i]);
	}
	ll _;
	struct pl
	{
		ll x,y;
		pl operator*(const pl&a)const{return {(x*a.x+y*a.y%mo*_)%mo,(x*a.y+y*a.x)%mo};}
	};
	bool ck(ll x) {return ksm(x,mo-1>>1)==1;}
	pl ksm(pl x,ll y) {pl a={1,0}; for(;y;y/=2,x=x*x)if(y&1)a=a*x; return a;}
	ll Sqrt(ll n)
	{
		if(!n) return 0;
		ll a;
		do a=R();
		while(ck(_=(a*a%mo-n+mo)%mo));
		return ksm({a,1},mo+1>>1).x;
	}
	void Sqrt(ll a[N],ll b[N],int n1)
	{
		static ll c[N],d[N];
		clr(b,n1); b[0]=Sqrt(a[0]);
		if(b[0]>mo-b[0]) b[0]=mo-b[0];
		for(int m=2;m<=n1;m*=2)
		{
			clr(c,m*2);
			go(i,0,m) c[i]=b[i]*2%mo;
			inv(c,d,m);
			clr(c,m*2); memcpy(c,a,8*m);
			mtp(c,d,m*2);
			go(i,m/2,m) b[i]=(b[i]*_2+c[i])%mo;
		}
	}
	void ln(ll a[N],ll b[N],int n)
	{
		static ll c[N];
		clr(b,n*2); inv(a,b,n);
		clr(c,n*2);
		go(i,1,n) c[i-1]=a[i]*i%mo;
		mtp(b,c,n*2);
		fd(i,n-1,1) b[i]=b[i-1]*iv(i)%mo;
		b[0]=0;
	}
	void exp(ll a[N],ll b[N],int n1)
	{
		static ll c[N];
		clr(b,n1); b[0]=1;
		for(int n=2;n<=n1;n*=2)
		{
			ln(b,c,n);
			go(i,0,n) c[i]=(a[i]-c[i]+mo)%mo;
			++c[0]%=mo;
			mtp(b,c,n*2); clr(b+n,n);
		}
	}
	void sin_cos(ll a[N],ll b[N],int n,bool f)
	{
		static ll c[N],b1[N];
		go(i,0,n) c[i]=a[i]*_1%mo;
		exp(c,b,n);
		go(i,0,n) T(c[i],mo-1);
		exp(c,b1,n);
		go(i,0,n1) b[i]=(f?b[i]+b1[i]:(b[i]-b1[i]+mo)*(mo-_1)%mo)*_2%mo;
	}
}
using namespace NTT;
int main()
{
	scanf("%d%d",&n,&k); p[0]=1;
	fo(i,1,k) scanf("%lld",p+i), p[i]=(mo-p[i])%mo;
	for(n1=1;n1<2*k;n1*=2); n2=n1*2;
	build(L);
	inv(p,t,n1-k);
	DFT(t,n2,0);
	reverse(p,p+k+1); DFT(p,n1,0);
	a[0]=b[1]=1;
	for(;n;n/=2,mul(b,b)) if(n&1)mul(a,b);
	go(i,0,k) scanf("%lld",&s), P(ans,a[i]*(mo+s));
	printf("%lld",ans);
}