多项式与快速傅里叶变换（FFT）

注意：以下内容将使用专用符号$i$表示$sqrt{-1}$

多项式

 一个以$x$为变量的多项式定义在一个代数域$F$上，将函数$A(x)$表示为形式和：

$$A(x)=sum_{j=0}^{n-1} a_{j} x^{j}$$

 我们称$a_{0}, a_{1}, cdots, a_{n-1}$为如上多项式的系数，所有系数都属于域$F$，典型的情形是复数集合$C$，如果一个多项式$A(x)$的最高次的非零系数是$a$，则称$A(x)$的次数是$k$，记$degree(A)=k$。任何严格大于一个多项式次数的整数都是该多项式的次数界，因此,对于次数界为$n$的多项式，其次数可以是$0$~$n-1$之间的任何整数,包括$0$和$n-1$。 对于多项式乘法,如果

$A(x)$和$B(x)$皆是次数界为$n$的多项式，则它们的乘积$C(x)$是一个次数界为$2n-1$的多项式，对所有属于定义域的$x$，都有$C(x)=A(x)B(x)$。有一种表示乘积$C(x)$的方法是：

$$C(x)=sum_{j=0}^{2 n-2} c_{j} x^{j}$$

其中，

$$c_{j}=sum_{k=0}^{j} a_{k} b_{j-k}$$

注意,，$degree(C)= degree(A)+ degree(B)$，意味着如果$A$是次数界为$n_a$的多项式，$B$是次数界为$n_b$的多项式，那么$C$是次数界为$n_a+n_b-1$的多项式。因为一个次数界为$k$的多项式也是次数界为$k+1$的多项式，所以通常称乘积多项式$C$是一个次数界为$n_a+n_b$的多项式。

多项式的不同表示方式

系数表示

对一个次数界为$n$的多项式$A(x)=sum_{j=0}^{n-1} a_{j} x^{j}$而言，其系数表达是一个由系数组成的向量$a=left(a_{0}, a_{1}, cdots, a_{n-1} ight)$。

某些运算在系数表示下非常方便，比如对多项式$A(x)$在给定点$x$的求值运算，使用霍纳法则，我们可以在$Theta(n)$时间复杂度内完成求值运算:

$$Aleft(x_{0} ight)=a_{0}+x_{0}left(a_{1}+x_{0}left(a_{2}+cdots+x_{0}left(a_{n-2}+x_{0}left(a_{n-1} ight) ight) cdots ight) ight)$$

类似地,对两个分别用系数向量$left(a_{0}, a_{1}, cdots, a_{n-1} ight)$和$b=left(b_{0}, b_{1}, cdots, b_{n-1} ight)$做乘法运算时，在$Thetaleft(n^{2} ight)$下得到向量$c=left(c_{0}, c_{1}, cdots, c_{n-1} ight)$，$c$称为$a$和$b$的卷积，记作$c=a otimes b$。

点值表示

一个次数界为$n$的多项式$A(x)$的点值表达就是一个由$n$个点值对所组成的集合

$$
left{left(x_{0}, y_{0} ight),left(x_{1}, y_{1} ight), cdots,left(x_{n-1}, y_{n-1} ight) ight}
$$

使得对$k=0,1, cdots, n-1$，所有$s_k$互不相同，

$$
y_{k}=Aleft(x_{k} ight)
$$

求值计算的逆(从一个多项式的点值表达确定其系数表达形式)称为插值。

 (插值多项式的唯一性)对于任意n个点值对组成的集合$left{left(x_{0}, y_{0} ight),left(x_{1}, y_{1} ight), cdots ight.left.left(x_{n-1}, y_{n-1} ight) ight}$其中所有的$x$都不同，那么存在唯一的次数界为$n$的多项式$A(x)$，满足$y_k=Aleft(x_{k} ight), k=0,1, cdots, n-1$。

证明主要是根据某个矩阵存在逆矩阵。原式等价于矩阵方程：

$$
left[egin{array}{ccccc}
{1} & {x_{0}} & {x_{0}^{2}} & {cdots} & {x_{0}^{m}} \
{1} & {x_{1}} & {x_{1}^{2}} & {cdots} & {x_{1}^{n-1}} \
{vdots} & {vdots} & {vdots} & {ddots} & {vdots} \
{1} & {x_{n-1}} & {x_{n-1}^{2}} & {cdots} & {x_{n-1}^{n-1}}
end{array} ight]left[egin{array}{c}
{a_{0}} \
{a_{1}} \
{vdots} \
{a_{n-1}}
end{array} ight]=left[egin{array}{c}
{y_{0}} \
{y_{1}} \
{vdots} \
{y_{n-1}}
end{array} ight]
$$

左边的矩阵表示为$Vleft(x_{0}, x_{1}, cdots, x_{n-1} ight)$，称为范德蒙德矩阵，该矩阵的行列式值为

$$
prod_{0 leq j<k leq n-1}left(x_{k}-x_{j} ight)
$$

因此，给定点值表达，我们能够唯一确定系数

$$
a=Vleft(x_{0}, x_{1}, cdots, x_{n-1} ight)^{-1} y
$$

一种更快的基于$n$个点的插值算法是基于拉格朗日公式

$$
A(x)=sum_{k=0}^{n-1} y_{k} frac{prod_{i eq k}left(x-x_{j} ight)}{prod_{j eq k}left(x_{k}-x_{j} ight)}
$$

给定$A$的扩展点值表达

$$
left{left(x_{0}, y_{0} ight),left(x_{1}, y_{1} ight), cdots,left(x_{2 n-1}, y_{2 n-1} ight) ight}
$$

和$B$的对应扩展点值表达

$$
left{left(x_{0}, y_{0}^{prime} ight),left(x_{1}, y_{1}^{prime} ight), cdots,left(x_{2 n-1}, y_{2 n-1}^{prime} ight) ight}
$$

则$C$的点值表达

$$
left{left(x_{0}, y_{0} y_{0}^{prime} ight),left(x_{1}, y_{1} y_{1}^{prime} ight), cdots,left(x_{2 n-1}, y_{2 n-1} y_{2 n-1}^{prime} ight) ight}
$$

DFT与FFT

如果使用单位复数根，可以在$Theta(n lg n)$完成求值与插值运算。

单位复数根

$n$次单位复数根是满足$omega ^{n}=1$的复数$omega$，$n$次单位复数根恰好有$n$个，对于$k=0,1, cdots ,n-1$，这些根是$e^{2 pi ik/n}$。为了解释这个表达式，我们使用复数的指数形式的定义：

$$e^{iu}= cos (u)+i sin (u)$$

$n$个单位复数根均匀地分布在以复平面的原点为圆心的单位半径的圆周上。值$omega _n=e^{2 pi i/n}$称为主$n$次单位根，所有其他$n$次单位负数根都是$omega_n$的幂次。

单位负数根在乘法意义下形成一个群，该群具有与加法群相同的结构，因为

$$omega _{n}^{j}omega _{n}^{k}=omega _{n}^{j+k}=omega _{n}^{left ( j+k ight ) mod n}$$

类似地，$omega _{n}^{-1}=omega _{n}^{n-1}$

下面将给出$n$次单位根的一些基本性质：

消去性质：$omega _{dn}^{dk}=omega _{n}^{k}$

对任意偶数$n>0$，有$omega _{n}^{n/2}=omega _{2}=-1$

折半性质：如果$n>0$为偶数，那么$n$个$n$次单位复数根的平方的集合就是$n/2$个$n/2次单位复数根的集合

求和性质：对任意整数$ngeq 1$和不能被$n$整除的非负整数$k$，有$sum _{j=0}^{n-1}left ( omega _{n}^{k} ight )^{j}=0$

DFT

我们希望计算次数界为$n$的多项式在$omega _{n}^{0},omega _{n}^{1},omega _{n}^{2},cdots ,omega _{n}^{n-1}$处的值

假设$A$以系数形式给出:$a=left ( a_{0},a_{1},cdots ,a_{n-1} ight )$，接下来对$k=0,1,cdots ,n-1$定义结果$y_k$:

$$y_{k}=Aleft ( omega _{n}^{k} ight )=sum _{j=0}^{n-1}a_{j}omega _{n}^{kj}$$

向量$y=left ( y_{0},y_{1},cdots ,y_{n-1} ight )$就是系数向量$a=left ( a_{0},a_{1},cdots ,a_{n-1} ight )$的离散傅里叶变换（DFT），记为$y=DFT_{n}left ( a ight )$

FFT

如果利用复数单位根的特殊性质，通过快速傅里叶变换的方法，就可以在$Theta left ( n lg n ight )$时间内计算出$DFT_{n}left ( a ight )$。

（以下的内容中假设$n$是$2$的整数幂）

$FFT$采用分治的策略，利用$A(x)$中偶数下标的系数与奇数下标的系数，分别定义两个新的次数界为$frac{n}{2}$的多项式$A^{left [ 0 ight ]}left ( x ight )$和$A^{left [ 1 ight ]}left ( x ight )$

$$A^{left [ 0 ight ]}left ( x ight )=a_{0}+a_{2}x+a_{4}x^{2}+cdots a_{n-2}x^{n/2-1}$$

$$A^{left [ 1 ight ]}left ( x ight )=a_{1}+a_{3}x+a_{5}x^{2}+cdots a_{n-1}x^{n/2-1}$$

注意到，$A^{left [ 0 ight ]}left ( x ight )$包含$A$中所有偶数下标的系数，以及$A^{left [ 1 ight ]}left ( x ight )$中包含$A$中所有奇数下标的系数。于是有，

$$Aleft ( x ight )=A^{left [ 0 ight ]}left ( x^{2} ight )+xA^{left [ 1 ight ]}left ( x^{2} ight )$$

所以，求$Aleft ( x ight )$在$omega _{n}^{0},omega _{n}^{1},omega _{n}^{2},cdots ,omega _{n}^{n-1}$处的值的问题转换为：

求次数界为$frac{n}{2}$的多项式$A^{left [ 0 ight ]}left ( x ight )$和$A^{left [ 1 ight ]}left ( x ight )$在点

$$left (omega _{n}^{0} ight )^{2},left (omega _{n}^{1} ight )^{2},cdots ,left (omega _{n}^{n-1} ight )^{2}$$

的取值。

于是原来的问题转化为了两个规模为原来一半的子问题。

在单位复数根处插值

 把$DFT$写成一个矩阵方程，然后再观察其逆矩阵的形式。

$$egin{bmatrix}y_{0}\ y_{1}\ y_{2}\ y_{3}\ vdots \ y_{n-1}end{bmatrix} = egin{bmatrix}1 & 1 & 1 & 1 & cdots & 1\ 1 & omega _{n} & omega _{n}^{2} & omega _{n}^{3} & cdots & omega _{n}^{n-1}\ 1 & omega _{n}^{2} & omega _{n}^{4} & omega _{n}^{6} & cdots & omega _{n}^{2left ( n-1 ight )}\ 1 & omega _{n}^{3} & omega _{n}^{6} & omega _{n}^{9} & cdots & omega _{n}^{3left ( n-1 ight )}\ vdots & vdots & vdots & vdots & ddots & vdots \ 1 & omega _{n}^{n-1} & omega _{n}^{2left ( n-1 ight )} & omega _{n}^{3left ( n-1 ight )} & cdots & omega _{n}^{left ( n-1 ight )left ( n-1 ight )}end{bmatrix} egin{bmatrix}a_{0}\ a_{1}\ a_{2}\ a_{3}\ vdots \ a_{n-1}end{bmatrix}$$

我们可以把$DFT$写成矩阵乘积$y=V_{n}a$，其中$V_n$是一个由$omega_n$适当幂次填充成的范德蒙德矩阵，对$j,k=0,1,cdots ,n-1,V_{n}$的$(k,j)$处元素为$omega _{n}^{kj}$。对于逆运算$a=DFT_{n}^{-1}left ( y ight )$，我们把$y$乘以$V_n$的逆矩阵$V_{n}^{-1}$来进行处理。

对$j,k=0,1,cdots ,n-1,V_{n}$，$V_n$的$(k,j)$处元素为$omega _{n}^{-kj}/n$，显然成立

这样快速傅里叶的内容就结束了：

void fft(complex<double>*a,int n,int inv)//快速傅里叶变换 
{
    for(int i=0;i<n;i++)//只有i<rev_i时才交换，防止交换两次 
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i*=2)//i表示要合并的序列长度的1/2 
    {
        complex<double>wn=exp(complex<double>(0,inv*pie/i));//求单位根w_n^1 
        for(int j=0;j<n;j+=i*2)//j表示合并到第几个序列 
        {
            complex<double> w(1,0);//求单位根 
            for(int k=j;k<j+i;k++)//k表示序列左半部分的指针，计算左半部分的同时计算右半部分 
            {
                complex<double>x=a[k];//蝴蝶变换 
                complex<double>y=w*a[k+i];
                a[k]=x+y;
                a[k+i]=x-y;
                w*=wn;//求w_n^k 
            }
        }
    }
    if(inv==-1)//若inv为-1表示要进行快速傅里叶逆变换 
        for(int i=0;i<n;i++)
            a[i]/=n;
}
void init()//读入并预处理 
{
    int lena=strlen(sa),lenb=strlen(sb);
    len=max(lena,lenb);
    for(int i=0;i<lena;i++)//逆序存储 
        an[i]=(double)(sa[lena-i-1]-'0');
    for(int i=0;i<lenb;i++)
        bn[i]=(double)(sb[lenb-i-1]-'0');
    while((1<<tot)<(len*2))//将多项式的次数界上调至2的n次幂 
    {
        tot++;
        s<<=1;
    }
    for(int i=0;i<s;i++)//预处理多项式中每一项的系数变换之后的位置：位逆序置换 
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(tot-1));
    fft(an,s,1);//对两个多项式求值 
    fft(bn,s,1);
    for(int i=0;i<s;i++)//点值乘法 
        an[i]*=bn[i];
    fft(an,s,-1);//对乘积插值 
    for(int i=0;i<s;i++)
    {
        ans[i]+=(int)(an[i].real()+0.5);//防止double类型精度损失 
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    while((!ans[s])&&s)
        s--;
    s++;
    while(s--)
        printf("%d",ans[s]);
}