多项式小结

Preface

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

万恶之源——FFT/NTT

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

[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^{2pi ij/n}=cosfrac{2pi j}n+isinfrac{2pi 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^2equiv A(mod x^{n/2}))，求(B^2equiv A(mod x^n))。
• 移项、平方，得到((B_0-B)^2equiv 0(mod x^n))。
• 拆开来、移项，得到(B_0^2+B^2equiv 2B_0B(mod x^n))。
• (Bequiv B_0/2+A/(2B_0)(mod x^n))。
• 注意如果常数项不为1，可能要打个二次剩余。（坑点：(x^2equiv n(mod 998244353))，同时((998244353-x)^2equiv n(mod 998244353)))

ln

• (G(x)=ln F(x)=intfrac{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_0equiv 0(mod x^{n/2}))，二次方后(mod x^n=0)，因此可以不管(i>1)的项。于是得到(F(B)equiv F'(B_0)(B-B_0)(mod x^n))。
• 因此(Bequiv 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 imes 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+isin 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=intfrac{F'}{sqrt{1-F^2}}dx)。这个exp都不用了，直接求逆+开根。
• 对于(y=arctan x)，有(y'=frac 1{1+x^2})。那么对于(G(x)=arctan F(x))，有(G=intfrac {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 iin[1,n/2]:B_0(x_i)=0)。
• 若我们将(A(x))表示成(D(x)B_0(x)+R(x))的形式，可得(forall iin[1,n/2]:A(x)=R(x))。
• 那么就可以用分治NTT+多项式取模解决。另一半同理。

快速插值

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

• 记(M(x)=prod_{i=1}^n(x-x_i))，洛一洛可得：(prod_{j eq i} (x_i - x_j) = lim_{x ightarrow 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 eq i}(x - x_j))。

• 再考虑分治。令(F_0(x)=sum_{i=1}^{n/2}frac{y_i}{M'(x_i)}prod_{j eq iwedge 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);
}