【学习笔记】之多项式使人头秃-白红宇

【学习笔记】之多项式使人头秃

阅读量：6069 次

发布时间：2019-06-20

本文共 6339 字，大约阅读时间需要 21 分钟。

真的自闭= =

多项式是什么鬼哦

首先介绍 FFT

~~我才不想写那么多柿子呢~~

大体说一下FFT干了啥

我们对两个多项式进行卷积（即多项式乘法）

$f=a*b$ 也就是 $f_i =\sum_{j=0}^i a_j * b_{i-j}$

暴力计算的话是n^2的

我们考虑把它变成点值[即(x,y)表示f(x)=y] 点值相乘就快了嘛但是变成点值了以后咋变回来呢

有个叫傅里叶的nb的人他发明了一个nb的东西叫傅里叶变换= =

也就是通过虚数中的单位根来计算就可以变回来了

单位根是什么东西呢就是在复平面上的一个单位圆将其弧等分成若干份第一个点位于(0,1)的n个点把这些数带进去就可以做啦

说起来很奇特对不对他其实就很奇特= =

具体详细的东西右转百度吧我实在是懒得写QAQ

实现上可以直接使用模板库里的complex（~~虽然我用起来非常不习惯~~）

扔个代码跑路。

#include
    
     #include
     
      #include
      
       #include
       
        #define maxn 3000010using namespace std;const double Pi = acos(-1.0);struct complex{    double x,y;    complex(double xx=0.0,double yy=0.0){x=xx,y=yy;}}A[maxn],B[maxn];int l,r[maxn],limit=1;complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}void FFT(complex *a,int type){    for(int i=0;i
        
         >1]>>1)|((i&1)<<(l-1));    FFT(A,1);FFT(B,1);    for(i=0;i<=limit;i++)        A[i]=A[i]*B[i];    FFT(A,-1);    for(i=0;i<=N+M;i++)        printf("%d ",(int)(A[i].x/limit+0.5));    return 0;}

然后我们就遇到了一个神奇的模数 998244353 ~~才不是1XXXXXX7~~

为什么是这个模数呢因为他是一个2^x* ...+1的一个素数具有一些优美的性质

我们就可以进行NTT[快速数论变换] /斜眼笑/

我们刚刚FFT中用的复平面中的单位根所以是有小数的

这个样子可不大好因为我们要取模所以我们有了一个很nb的东西叫做原根

原根有一些优美的性质就跟单位根一样 G ^ ((p-1)/i) 就是可以当成单位根使用的数啦 [i|(p-1)]这也就是p为什么要是 2^x *... +1的原因啦

小姿势：998244353的原根是3

其他详细的细节还是右转百度吧【大雾

扔个代码。

#include
    
     #include
     
      #include
      
       #define maxn 300005#define modn 998244353#define G 3#define ll long longusing namespace std;int q_pow(ll base,ll pow){	ll ans=1;	while(pow)	{		if(pow&1){ans*=base;ans%=modn;}		base*=base;base%=modn;pow>>=1;	}	return (int)ans;}int A[maxn],B[maxn],C[maxn];int l,r[maxn],limit,inv;void FFT(int *a,int type){	int i,j,k;	for(i=0;i
       
        i)			swap(a[r[i]],a[i]);	for(i=2;i<=limit;i<<=1)	{		int mid=i>>1;		int Wn=q_pow(G,(modn-1)/i);		if(type)	Wn=q_pow(Wn,(modn-2));		for(j=0;j
        
         =modn)	a[j+k]-=modn;				a[j+k+mid]=x-(ll)w*y%modn;				if(a[j+mid+k]<0)	a[j+mid+k]+=modn;				w=(ll)w*Wn%modn;			}		}	}	if(type)	{		for(i=0;i
         
          >1]>>1)|((i&1)<<(l-1)); inv=q_pow(limit,modn-2); FFT(A,0);FFT(B,0); for(i=0;i<=limit;i++) C[i]=(ll)A[i]*B[i]%modn; FFT(C,1); for(i=0;i<=(n+m);i++) printf("%d\n",C[i]); return 0;}

于是我的姿势还只停留在 FFT/NTT 只是能求个卷积【大雾

然后就有一些奇奇怪怪的东西了=.=+

【奇奇怪怪的东西一】多项式求逆

我们现在有一个多项式f 我们要求一个多项式g满足 $f * g \equiv 1\ (mod\ x^n)$

这玩意看上去是不是非常奇怪啊~~【明明就是非常奇怪！~~

假设我们现在已知一个多项式h 满足 $f * h \equiv 1\ (mod\ x^{\lceil \frac{n}{2}\rceil})$

我们可以得到 $g-h\equiv 0\ (mod\ x^{\lceil \frac{n}{2}\rceil})$

平方 $g^2 -2g*h +h^2 \equiv 0 \ (mod\ x^n)$

卷上f $g - 2h\ + fh^2 \equiv 0\ (mod\ x^{n})$

移项 $g = 2h -fh^2 \ (mod \ x ^n)$

递归求解就好啦

边界当然是n=1的时候 g直接取f的常数项的逆元啦qwq。

附代码。

#include
    
     #include
     
      #include
      
       #include
       
        #define inf 20021225#define ll long long#define mdn 998244353#define G 3#define mxn 300100using namespace std;int rev[mxn];int ksm(int bs,int mi){	int ans=1;	while(mi)	{		if(mi&1)	ans=(ll)ans*bs%mdn;		bs=(ll)bs*bs%mdn;mi>>=1;	}	return ans;}int inv;void NTT(int *a,int lim,int f){	for(int i=0;i
        
         i)	swap(a[i],a[rev[i]]);	for(int k=2;k<=lim;k<<=1)	{		int Wn=ksm(G,(mdn-1)/k),mid=k>>1;		if(f)	Wn=ksm(Wn,mdn-2);		for(int w=1,i=0;i
         
          >1); int lim=1,l=0; while(lim<(n<<1)) lim<<=1,l++; for(int i=0;i
          
           >1]>>1)|((i&1)<<(l-1)); inv=ksm(lim,mdn-2); for(int i=0;i

【奇奇怪怪的东西二】多项式对数函数

看上去是不是很高大上！~~【实则蠢得一批~~

对于多项式 f 求g=ln f

这之前我们科普一点求导和积分的小姿势

对于一个普通多项式

求导 $f'(x) = \sum_{i=0}^{n-1} (i+1)a_{i+1} x^i$

积分 $\int x^a dx = \frac{1}{a+1} x^{a+1}$

两个过程都很像哒就是反过来做而已233

ln x求导是 1/x

复合函数求导 $(g(f(x)))'= g'(f(x))f(x)$

然后ln f的求导

$(\textup{ln} f(x))' =\frac{1}{f(x)}*f'(x)$

直接多项式求逆然后求导再积分回去就好啦qwq

代码等我一哈【不咕不咕必定不可能咕

update：真的没有咕！

#include
    
     #include
     
      #include
      
       #include
       
        #define inf 20021225#define ll long long#define mdn 998244353#define G 3#define mxn 300100using namespace std;int ksm(int bs,int mi){	int ans=1;	while(mi)	{		if(mi&1)	ans=(ll)ans*bs%mdn;		bs=(ll)bs*bs%mdn;mi>>=1;	}	return ans;}int inv,rev[mxn];void NTT(int *a,int lim,int f){	for(int i=0;i
        
         i)	swap(a[i],a[rev[i]]);	for(int k=2;k<=lim;k<<=1)	{		int mid=k>>1,Wn=ksm(G,(mdn-1)/k);		if(f)	Wn=ksm(Wn,mdn-2);		for(int w=1,i=0;i
         
          >1]>>1)|((i&1)<
          
           >1);int lim=1,l=0; while(lim<(n<<1)) lim<<=1,l++; for(int i=0;i

【奇奇怪怪的东西三】多项式指数函数

对于给定 f(x) 求h(x)= e^(f(x)) mod(x^n)

有点鬼畜是不是= =||

还是先讲一些前缀姿势

1. 泰勒展开

对于一个函数 f(x) 我们可以使用高阶导数对其无限逼近

$f(x) = f(a) +f'(a)\frac{(x-a)}{1!}+f''(a)\frac{(x-a)^2}{2!}+...$

2.牛顿迭代

我们现在要求 $g(f(x)) \equiv 0\ (mod\ x^n)$ 其中g已知

假设我们知道原式 $mod\ x^{\lceil \frac{n}{2} \rceil}$ 的答案为 $f_0(x)$

根据泰勒展开

$0= g(f)=g(f_0)+g'(f_0)(f-f_0)+g''(f_0)\frac{(f-f_0)^2}{2!}+...$

显然 $(f-f_0)^2 \equiv 0 (mod\ x^n)$ 所以从第三项开始都是模x^n 为0的我们可以不用考虑

那么就是 $g(f_0)+g'(f_0)(f-f_0) = 0 (mod\ x^n)$

拆开移项得 $f = f_0- \frac{g(f_0)}{g'(f_0)}$

那么我们就可以递归求解啦= =+

诶你突然发现问题了对不对= =+

g我们好像还不知道是什么呢

g当然是我们自己构造的啦由于 $h(x) \equiv e^{f(x)} (mod \ x^n)$ 两边同时取ln

得到 $\textup{ln} h -f =0$

根据牛顿迭代的柿子 $g= \textup{ln} h - f$

所以把式子带进去 $h = h_0 - f =h_0 - \frac{\textup{ln} h_0 - f}{\textup{ln} h_0}= h_0*(1- \textup{ln}h_0 +f)$

这次就没问题啦= =+

~~【喝了口开水差点被烫死】~~

#include
    
     #include
     
      #include
      
       #include
       
        #define inf 20021225#define ll long long#define mxn 400100#define mdn 998244353#define G 3using namespace std;int ksm(int bs,int mi){	int ans=1;	while(mi)	{		if(mi&1)	ans=(ll)ans*bs%mdn;		bs=(ll)bs*bs%mdn;mi>>=1;	}	return ans;}int rev[mxn],inv;int init(int n){	int lim=1,l=0;	while(lim
        
         >1]>>1)|((i&1)<
         
          i) swap(a[rev[i]],a[i]); for(int k=2;k<=lim;k<<=1) { int mid=k>>1,Wn=ksm(G,(mdn-1)/k); if(f) Wn=ksm(Wn,mdn-2); for(int w=1,i=0;i
          
           >1); int lim=init(n<<1); for(int i=0;i
           
            >1;poly_exp(a,mid); for(int i=0;i<(n<<1);i++) e[i]=g[i]=0; poly_ln(s,n); for(int i=0;i

【奇奇怪怪的东西五】多项式除法

听起来是不是十分酷炫实际上用到一个神奇的转化就可以轻轻松松~~松松松~~（大雾）的通过啦

我们对于给定 n次多项式f(x) 和 m次多项式g(x) 求 $f(x)=g(x)*d(x)+r(x)$ 中的 d 和 r 没错就是平常见过的多项式除法其中（n>m）

这个玩意如何实现呢？直接做的话肯定是O(nm) 十分不优秀况且既然有了 FFT这种nb的东西怎么能让这么不优美的复杂度存在呢

下面来介绍这个神奇的操作

我们将一个多项式的系数翻转【没错就是前后倒过来那个翻转】

这个应该怎么在数学中实现呢就是这个样子

$f^R(x) = x^n f(\frac{1}{x})$

然后我们来干一些有趣的事情把前面的所有多项式中的x替换成 $\frac{1}{x}$ 然后等式两边同时乘 $x^n$ 于是就有了

$x^nf(\frac{1}{x})=x^mg(\frac{1}{x})x^{n-m}d(\frac{1}{x})+x^{n-m+1}x^{m-1}r(\frac{1}{x})$

然后我们发现这个玩意很优美可以化成

$f^R(x)=g^R(x)d^R(x)+x^{n-m+1}r^R(x)$

我们发现余数项的最低次都是n-m+1 所以我们可以让整个柿子对 $x^{n-m+1}$ 取模来消除r对答案的影响。

然后我们就可以鱼块的多项式求逆来求出 $d^R(x)$ 再翻转回来得到 $d(x)$ 带回原式把 $r(x)$ 求出来就做完啦

所以其实比前面的还要好写= =+

代码。

#include
    
     #include
     
      #include
      
       #include
       
        #define inf 20021225#define ll long long#define mdn 998244353#define mxn 400100#define G 3using namespace std;int ksm(int bs,int mi){    int ans=1;    while(mi)    {        if(mi&1)    ans=(ll)ans*bs%mdn;        bs=(ll)bs*bs%mdn;mi>>=1;    }    return ans;}int inv,rev[mxn];void NTT(int *a,int lim,int f){    for(int i=1;i
        
         i)    swap(a[rev[i]],a[i]);    for(int k=2;k<=lim;k<<=1)    {        int mid=k>>1,Wn=ksm(G,(mdn-1)/k);        if(f)   Wn=ksm(Wn,mdn-2);        for(int i=0,w=1;i
         
          >1]>>1)|((i&1)<
          
           >1;poly_inv(a,mid); int lim=init(n<<1); for(int i=0;i