离散傅里叶变换

定义:

是将 $\omega_n^1,\omega_n^2,…,\omega_n^n$ 带入多项式所得到的一种特殊的多项式点值表示

性质:

考虑将 $A(x)$ 的每一项按下标的奇偶分成两个多项式之和:

$A(x) = (a_0 + a_2x^2 + … + a_{n – 2}x^{n – 2}) + (a_1x + a_3x^3 + … + a_{n-1}x^{n-1})$

如果设:

$A_1(x) = a_0 + a_2x + … + a_{n – 2}x^{\frac{n}{2} – 1}$

$A_2(x) = a_1 + a_3x + … + a_{n – 1}x^{\frac{n}{2} – 1}$

那么:

$A(x) = A_1(x^2) + xA_2(x^2)$

于是我们把 $x = \omega_n^k$ 带入$A(x)$:

当 $k\subseteq \lbrace 0,1,2,… ,\frac{n}{2}-1\rbrace $ 时:

$A(\omega_{n}^{k})=A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$

$A(\omega_{n}^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^{k}) – \omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$

这样如果我们知道了 $A_1(x),A_2(x)$ 的离散傅里叶变换表示就可以 $O(n)$ 地求出 $A(x)$ 的离散傅里叶变换。

快速傅里叶变换

根据上面所提的离散傅里叶变换的性质,我们就可以递归地在 $O(nlogn)$的时间内将多项式的系数表示转换为点值表示。下面就介绍一下优化和实现上的细节:

蝴蝶操作:

有一个很神奇的性质:将一个多项式按下标的奇偶分为两部分,并将奇数部分放在偶数部分之后,然后不断递归下去。那么多项式下标为 $i$ 的数在这个操作之后的下标就变成了「将 $i$ 二进制翻转所得的数」。

for(rint i=0; i<lim; ++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(l-1));

根据这个性质,我们可以先预处理出多项式中每个数操作之后的位置,并将它放到这个位置。然后类似于递归时回溯的操作,不断向上还原求出点值表示。

for(rint i=0; i<lim; ++i)
    if(i<pos[i]) swap(a[i], a[pos[i]]);
for(rint i=1; i<lim; i<<=1)
{
    Complex del(cos(PI/i), sin(PI/i));
    for(rint j=0; j<lim; j+=(i<<1))
    {
        Complex t(1, 0);
        for(rint k=0; k<i; ++k, t=t*del)
        {
            Complex x=a[j+k], y=t*a[j+k+i];
            a[j+k]=x+y; a[j+k+i]=x-y;
        }
    }
}

傅里叶逆变换

按照上述方法我们在 $O(nlogn)$ 内将多项式的系数表示转换成了点值表示,然后我们在 $O(n)$ 的时间内将两个多项式的点值表示相乘,得到了乘积的点值表示。利用傅里叶逆变换就可以同样在 $O(nlogn)$ 内将点值表示转换成系数表示。

离散傅里叶变换的另一个性质:把多项式 $A(x)$ 的离散傅里叶变换作为另一个多项式 $B(x)$ 的系数,并将 $\omega_{n}^{0}, \omega_{n}^{-1}, \omega_{n}^{-2}, …, \omega_{n}^{-(n – 1)}$ 代入多项式 $B(x)$ ,得到的每个数除 $n$,这样的 $n$ 个数就恰好是多项式 $A(x)$ 的系数。

根据上面的性质,我们在求离散傅里叶变换的代码上稍作修改,就能在 $O(nlogn)$ 内将点值表示转换成系数表示。

最终代码

下面是用 FFT 优化的高精度乘法代码(BZOJ2179

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#include <vector>
#include <map>
#include <set>
#define MAXN 200005
#define INF 0x3f3f3f3f
#define rint register int
#define LL long long
#define LD long double
using namespace std;

const double PI=acos(-1);

int n, m, len, pos[MAXN], ans[MAXN];
char sa[MAXN], sb[MAXN];

struct CP
{
    double x, y;
    CP(double a=0, double b=0)
    {
        x=a; y=b;
    }
    friend CP operator * (CP a, CP b)
    {
        return CP(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
    }
    friend CP operator + (CP a, CP b)
    {
        return CP(a.x+b.x, a.y+b.y);
    }
    friend CP operator - (CP a, CP b)
    {
        return CP(a.x-b.x, a.y-b.y);
    }
} a[MAXN], b[MAXN];

void fft(CP a[], int op)
{
    for(rint i=0; i<n; ++i)
        if(i<pos[i]) swap(a[i], a[pos[i]]);
    for(rint i=1; i<n; i<<=1)
    {
        CP omg(cos(PI/i), op*sin(PI/i));
        for(rint j=0; j<n; j+=(i<<1))
        {
            CP t(1, 0);
            for(rint k=0; k<i; ++k, t=t*omg)
            {
                CP x=a[j+k], y=t*a[j+k+i];
                a[j+k]=x+y; a[j+k+i]=x-y;
            }
        }
    }
}

int main()
{
    scanf("%d", &len);
    scanf("%s%s", sa, sb);
    for(rint i=0; i<len; ++i)
    {
        a[i].x=sa[len-1-i]-'0';
        b[i].x=sb[len-1-i]-'0';
    }
    for(n=1; n<(len*2); ) n<<=1, m++;
    for(rint i=0; i<=n; ++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(m-1));
    fft(a, 1);
    fft(b, 1);
    for(rint i=0; i<n; ++i) a[i]=a[i]*b[i];
    fft(a, -1);
    for(rint i=0; i<n; ++i)
    {
        ans[i]+=(int)(a[i].x/n+0.5);
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    len=len*2-1;
    while(!ans[len] && len>=1) len--;
    for(rint i=len; i>=0; i--) printf("%d", ans[i]);
    return 0;
}

例题

1.BZOJ3771 Triple
2.BZOJ4259 残缺的字符串
3.AGC005F Many Easy Problems

1 对 “快速傅里叶变换小结”的想法;

发表评论

邮箱地址不会被公开。 必填项已用*标注