多项式

多项式

前置知识

多项式的度：

多项式$f(x)$的最高次项的次数为该多项式的度，记为$degf$。

多项式的逆元：

对于多项式$f(x)$，若存在$g(x)$满足：
$f(x)g(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
$degg\le degf$
则称$g(x)$为$f(x)$在模$x^n$意义下的逆元，记作$f^{-1}(x)$

多项式的余数和商：

对于多项式$f(x),g(x)$，存在唯一的$Q(x),R(x)$满足：
$f(x)=Q(x)g(x)+R(x)$
$degQ=degf-degg$
$degR<degg$
$f(x)\equiv R(x)(\mathrm{m}\mathrm{o}\mathrm{d}g(x))$
$Q(x)$为商，$R(x)$为余数

多项式的多点求值和插值

多点求值：给出一个多项式$f(x)$和$n$个点$x_1,x_2,\dots ,x_n$,求$$f(x_1),f(x_2),\dots ,f(x_n)$$。
插值：给出$n+1$个点$$(x_0,y_0),(x_1,y_1),\dots ,(x_n,y_n)$$
求一个$n$次多项式$f(x)$使得这$n+1$个点都在$f(x)$上。

拉格朗日插值

给出$n$个点$P_i(x_i,y_i)$，将过这$n$个点的最多$n-1$次多项式记为$f(x)$,求$f(k)$的值。
由这$n$个点可以构建$n$个简单多项式$g_i(x)$,使得$g_i(x_i)=y_i\&\&g_i(x_j)=0,j\ne i$
那么$f(x)={\sum }_{i=1}^n{g}_i(x)$
很明显，$g_i(x)=y_i\times (\prod _{j\ne i}\frac{x-x_j}{x_i-x_j})$
所以$$f(x)=\sum _{i=1}^n{y}_i\times (\prod _{j\ne i}\frac{x-x_j}{x_i-x_j})$$
代入$k$,$$f(k)=\sum _{i=1}^n{y}_i\times (\prod _{j\ne i}\frac{k-x_j}{x_i-x_j})$$

快速傅里叶变换（FFT）

系数表示法

f ( x ) = a 2 x 2 + a 1 x + a 0 ⇔ f ( x ) = { a 2 , a 1 , a 0 } f(x)=a_2x^2+a_1x+a_0\Leftrightarrow f(x)=\{a_2,a_1,a_0\} f(x)=a2​x2+a1​x+a0​⇔f(x)={a2​,a1​,a0​}

点值表示法

从多项式$f(x)$上取$n+1$个有效点，来唯一表示这个函数。
如下：
$$f_1(x)=y_1=a_0+a_1{x}_1+a_2{x}_1^2+\dots +a_n{x}_1^n$$
$$f_2(x)=y_2=a_0+a_1{x}_2+a_2{x}_2^2+\dots +a_n{x}_2^n$$
$$f_3(x)=y_3=a_0+a_1{x}_3+a_2{x}_3^2+\dots +a_n{x}_3^n$$
$$⋮$$
$$f_{n+1}(x)=y_1=a_0+a_1{x}_{n+1}+a_2{x}_{n+1}^2+\dots +a_n{x}_{n+1}^n$$

DFT

多项式由系数表示法转为点值法的过程

IDFT

多项式的点值表示法转化为系数表示法的过程

FFT

通过取某些特殊的$x$的点值来加速DFT和IDFT的过程

单位复根

两个点值表示的多项式乘积复杂度为$O(n)$，如下：
假设我们DFT过程对于两个多项式选取的$x$序列相同，那么可以得到：
$$f(x)=(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),\dots ,(x_n,f(x_n))$$
$$g(x)=(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),\dots ,(x_n,g(x_n))$$
如果设$F(x)=f(x)\times g(x)$:
则$F(x)$的点值表达式为：$$F(x)=(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),\dots ,(x_n,f(x_n)g(x_n))$$
如果能快速得到$F(x)$的系数表达式，则这个问题就完美解决了。
我们可以观察到，为了计算$F(x)$，就需要计算很多$x^i$，为了能快速得到$x^i$，引入${\omega }_n^k,0\le k<n$,表示复数域上单位圆被$n$等分后的第$k$个，从0开始。由复数的性质可知$({\omega }_n^k{)}^{n\ast i}=1$,${\omega }_n^k=({\omega }_n^1{)}^k$.

FFT的流程

DFT

例子：
$f(x)=a_0+a_{1}x+a_2{x}^2+\dots +a_7{x}^7$
$f(x)=(a_0+a_2{x}^2+a_4{x}^4+a_6{x}^6)+x(a_1+a_3{x}^2+a_5{x}^4+a_7{x}^6)$
令：$G(x)=a_0+a_{2}x+a_4{x}^2+a_6{x}^3$
$H(x)=a_1+a_{3}x+a_5{x}^2+a_7{x}^3$
则$F(x)=G(x^2)+xH(x^2)$
即有：$DFT(f({\omega }_n^k))=DFT(G(({\omega }_n^k{)}^2))+{\omega }_n^{k}DFT(H(({\omega }_n^k{)}^2))$
这个函数能处理的多项式长度只能是$2^m$,高次系数补0即可。这个代码还需要继续优化。
设初始序列为 { x 0 , x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 } \{x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7\} {x0​,x1​,x2​,x3​,x4​,x5​,x6​,x7​}
一次二分后为 { x 0 , x 2 , x 4 , x 6 } , { x 1 , x 3 , x 5 , x 7 } \{x_0,x_2,x_4,x_6\},\{x_1,x_3,x_5,x_7\} {x0​,x2​,x4​,x6​},{x1​,x3​,x5​,x7​}
二次二分后为 { x 0 , x 4 } , { x 2 , x 6 } , { x 1 , x 5 } , { x 3 , x 7 } \{x_0,x_4\},\{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\} {x0​,x4​},{x2​,x6​},{x1​,x5​},{x3​,x7​}
三次二分后为 { x 0 } , { x 4 } , { x 2 } , { x 6 } , { x 1 } , { x 5 } , { x 3 } , { x 7 } \{x_0\},\{x_4\},\{x_2\},\{x_6\},\{x_1\},\{x_5\},\{x_3\},\{x_7\} {x0​},{x4​},{x2​},{x6​},{x1​},{x5​},{x3​},{x7​}
其实就是原来的那个序列，每个数用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标。

IDFT

由于写成矩阵形式后，
$\left[\begin{array}{c}y[0]\\ y[1]\\ y[2]\\ y[3]\\ \dots \\ y[n-1]\end{array}\right]\begin{array}{c}=\\ =\\ =\\ =\\ \\ =\end{array}\left[\begin{array}{cccccc}1& 1& 1& 1& \dots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \dots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \dots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \dots & {\omega }_n^{3(n-1)}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \dots & {\omega }_n^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}a[0]\\ a[1]\\ a[2]\\ a[3]\\ \dots \\ a[n-1]\end{array}\right]$
这个矩阵是范德蒙德矩阵，其逆矩阵是每一项求倒数后除以$n$，然后发现${\omega }_n^i$倒数是就是${\omega }_n^{-i}$。
所以把fft多加一个参数，$on==1$时是FFT，$on==-1$时是IFFT，只要构造单位复根时取${\omega }_n=\cos (on\ast 2\pi /n)+i\sin (on\ast 2\pi /n)$,就完成了FFT和IFFT的构造。

///FFT
///高精度乘法
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;

const double PI = acos(-1.0);
//复数
struct Complex{
    double x,y;
    Complex(double _x = 0.0,double _y = 0.0){
        x=_x;y=_y;
    }
    Complex operator-(const Complex &b) const {
        return Complex(x-b.x,y-b.y);
    }
    Complex operator+(const Complex &b) const {
        return Complex(x+b.x,y+b.y);
    }
    Complex operator*(const Complex &b) const {
        return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
    }
}; 
/*
 *进行FFT和IFFT前的反置变换
 *位置i和i的二进制翻转后的位置互换
 *len必须为2的幂次
 */
void change(Complex y[],int len){
    int i,j,k;
    for(int i=1,j=len/2;i<len-1;i++){
        if(i<j)swap(y[i],y[j]);
        k=len/2;
        while(j>=k){
            j=j-k;
            k=k/2;
        }
        if(j<k) j+=k;
    }
}
/*
 *FFT
 *len必须是2的幂次
 *on==1时是DFT，on==-1时是IDFT
 */
void fft(Complex y[],int len ,int on){
    change(y,len);//翻转系数
    //倍增法
    for(int h=2;h<=len;h<<=1){
        //取单位复根
        Complex wn(cos(2*PI/h),sin(on*2*PI/h));
        for(int j=0;j<len;j+=h){
            Complex w(1,0);
            for(int k=j;k<j+h/2;k++){
                Complex u=y[k];
                Complex t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    //IFFT时，需要将结果/len
    if(on==-1){
        for(int i=0;i<len;i++){
            y[i].x/=len;
        }
    }
}

const int maxn = 200020;
Complex x1[maxn],x2[maxn];
char str1[maxn/2],str2[maxn/2];
int sum[maxn];

int main(){
    while(scanf("%s%s",str1,str2)==2){
        int len1=strlen(str1);
        int len2=strlen(str2);
        int len=1;
        while(len<len1*2||len<len2*2)len<<=1;
        for(int i=0;i<len1;i++)
            x1[i]=Complex(str1[len1-1-i]-'0',0);
        for(int i=len1;i<len;i++)
            x1[i]=Complex(0,0);
        for(int i=0;i<len2;i++)
            x2[i]=Complex(str2[len2-1-i]-'0',0);
        for(int i=len2;i<len;i++)
            x2[i]=Complex(0,0);
        fft(x1,len,1);
        fft(x2,len,1);
        //转化为点值表示法后计算乘积
        for(int i=0;i<len;i++)
            x1[i]=x1[i]*x2[i];
        //IDFT
        fft(x1,len,-1);
        for(int i=0;i<len;i++)
            sum[i]=int(x1[i].x + 0.5);
        for(int i=0;i<len;i++){
            sum[i+1]+=sum[i]/10;
            sum[i]%=10;
        }
        len=len1+len2-1;
        while(sum[len]==0&&len>0)len--;
        for(int i=len;i>=0;i--)
            printf("%c",sum[i]+'0');
        printf("\n");
    }
    return 0;
}

快速数论变换（NTT）

生成子群

群：一个满足封闭性、满足结合律、有单位元、有逆元的二元运算的代数结构。
子群：群$(S,\oplus ),S(S^{\prime },\oplus )$，满足$S^{\prime }\subset S$,则$(S^{\prime },\oplus )$是$(S,\oplus )$的子群
拉格朗日定理：$|S^{\prime }|$是$|S|$的因子
生成子群：$a\in S$的生成子群 < a > = { a ( k ) , k ≥ 1 } \left<a\right>=\{a^{(k)},k\ge 1\} ⟨a⟩={a(k),k≥1}，$a$是$⟨a⟩$的生成元
阶：群$S$中$a$的阶是满足$a^r=e$,($e$是单位元)，的最小的$r$,记为$ord(a)=|⟨a⟩|$
考虑群 Z n × = { [ a ] , n ∈ Z n : g c d ( a , n ) = 1 } , ∣ Z n × ∣ = φ ( n ) Z_n^{\times}=\{[a],n\in Z_n:gcd(a,n)=1\},|Z_n^{\times}|=\varphi(n) Zn×​={[a],n∈Zn​:gcd(a,n)=1},∣Zn×​∣=φ(n)

原根

$g$满足$or{d}_n(g)=|Z_n^{\times }|=\phi (n)$，对于质数$p$，也就是说$g^i\mathrm{m}\mathrm{o}\mathrm{d}p,0\le i<p$时的结果互不相同。
模$n$有原根的充要条件：$n=2,4,p^a,2\ast p^a$
离散对数：$g^t\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}n),in{d}_{n,g}(a)=t$
因为$g$是原根，所以$gt$每$\phi (n)$一个周期，可以取到$|Z\times n|$的所有元素对于$n$是质数时，就是得到$[1,n-1]$的所有数，就是$[0,n-2]$到$[1,n-1]$的映射，离散对数满足对数的相关性质。
求原根可以证明满足$g^r\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$的最小的$r$一定是$p-1$的约数，对于质数$p$，质因子分解$p-1$，若$g^{(p-1)/pi}\ne 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$恒成立，$g$为$p$的原根

NTT

对于质数$p=qn+1,(n=2^m)$，原根$g$满足$g^{qn}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$看作${\omega }_n$的等价，则其满足相似的性质，比如$g_n^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p),g_n^{n/2}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$
可以将$g^{qn}$等价为$e^{2\pi n}$，迭代到长度$l$时，${\omega }_n=g_l=g_{p-1}^{\frac{p-1}{l}}$。

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#include<iomanip>
#include<vector>
#include<string>
#include<bitset>
#include<queue>
#include<map>
#include<set>
using namespace std;
 
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=10*x+ch-'0';ch=getchar();}
	return x*f;
}
void print(int x)
{if(x<0)putchar('-'),x=-x;if(x>=10)print(x/10);putchar(x%10+'0');}
 
const int N=300100,P=998244353;
 
inline int qpow(int x,int y)
{
	int res(1);
	while(y)
	{
		if(y&1) res=1ll*res*x%P;
		x=1ll*x*x%P;
		y>>=1;
	}
	return res;
}
 
int r[N];
 
void ntt(int *x,int lim,int opt)
{
	register int i,j,k,m,gn,g,tmp;
	for(i=0;i<lim;++i)
		if(r[i]<i)
			swap(x[i],x[r[i]]);
	for(m=2;m<=lim;m<<=1)
	{
		k=m>>1;
		gn=qpow(3,(P-1)/m);
		for(i=0;i<lim;i+=m)
		{
			g=1;
			for(j=0;j<k;++j,g=1ll*g*gn%P)
			{
				tmp=1ll*x[i+j+k]*g%P;
				x[i+j+k]=(x[i+j]-tmp+P)%P;
				x[i+j]=(x[i+j]+tmp)%P;
			}
		}
	}
	if(opt==-1)
	{
		reverse(x+1,x+lim);
		register int inv=qpow(lim,P-2);
		for(i=0;i<lim;++i)
			x[i]=1ll*x[i]*inv%P;
	}
}
 
int A[N],B[N],C[N];
 
char a[N],b[N];
 
int main()
{
	register int i,lim(1),n;
	scanf("%s",&a);
	n=strlen(a);
	for(i=0;i<n;++i) A[i]=a[n-i-1]-'0';
	while(lim<(n<<1)) lim<<=1;
	scanf("%s",&b);
	n=strlen(b);
	for(i=0;i<n;++i) B[i]=b[n-i-1]-'0';
	while(lim<(n<<1)) lim<<=1;
	for(i=0;i<lim;++i)
		r[i]=(i&1)*(lim>>1)+(r[i>>1]>>1);
	ntt(A,lim,1);ntt(B,lim,1);
	for(i=0;i<lim;++i)
		C[i]=1ll*A[i]*B[i]%P;
	ntt(C,lim,-1);
	int len(0);
	for(i=0;i<lim;++i)
	{
		if(C[i]>=10)
			len=i+1,
			C[i+1]+=C[i]/10,C[i]%=10;
		if(C[i]) len=max(len,i);
	}
	while(C[len]>=10)
		C[len+1]+=C[len]/10,C[len]%=10,len++;
	for(i=len;~i;--i)
		putchar(C[i]+'0');
	puts("");
	return 0;
}

快速沃尔什变换（FWT）