HDU 1402 快速傅里叶变换FFT

最新推荐文章于 2020-12-25 22:30:26 发布

韩小侠

最新推荐文章于 2020-12-25 22:30:26 发布

阅读量1.1k

点赞数

分类专栏：数论

本文链接：https://blog.csdn.net/prime7/article/details/9470957

版权

数论专栏收录该内容

69 篇文章 1 订阅

订阅专栏

题意：大数乘法。

看了算法导论（第2版）第30章的《多项式与快速傅里叶变换》

多项式有系数表示法和点值表示法。

两个次数界为n的多项式A（x）和B（x）相乘，输入输出均采用系数表示法。（假定n为2的幂）

1）使次数界增加一倍：A（x）和B（x）扩充为次数界为2n的多项式，并构造起系数表示

2）求值：两次应用2n阶FFT，计算出A（x）和B（x）的长度为2n的点值表示

3）点乘：计算多项式C(x)=A（x）*B（x）的点值表示

4）插值：对2n个点值对应用一次FFT计算出其逆DFT，就可以构造出多项式C（x）的系数表示

第1步和第3步需要时间为O(n)，第2步和第4步运用FFT需要时间为O(nlgn)。

第2步求取n个点的值所需时间为O(n^2)，如果我们巧妙地选取x（k）则其运行时间变为O(nlgn)，FFT主要利用了单位复根（w^n=1的w就是n次单位复根，他们是e^(2PI*k/n),k=0,1,……n-1）的特殊性质。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define N 200005
#define pi acos(-1.0) // PI值
using namespace std;
struct complex
{
    double r,i;
    complex(double real=0.0,double image=0.0)
    {
        r=real;
        i=image;
    }
    // 以下为三种虚数运算的定义
    complex operator + (const complex o)
    {
        return complex(r+o.r,i+o.i);
    }
    complex operator - (const complex o)
    {
        return complex(r-o.r,i-o.i);
    }
    complex operator * (const complex o)
    {
        return complex(r*o.r-i*o.i,r*o.i+i*o.r);
    }
} x1[N],x2[N];
int sum[N]; // 结果存在sum里
void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
{
    register int i,j,k;
    for(i=1,j=l/2; i<l-1; i++)
    {
        if(i<j)    swap(y[i],y[j]); // 交换互为下标反转的元素
        // i<j保证只交换一次
        k=l/2;
        while(j>=k) // 由最高位检索，遇1变0，遇0变1，跳出
        {
            j-=k;
            k/=2;
        }
        if(j<k)    j+=k;
    }
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
// 其中on==1时为DFT，on==-1为IDFT
{
    register int h,i,j,k;
    complex u,t;
    brc(y,l); // 调用反转置换
    for(h=2; h<=l; h<<=1) // 控制层数
    {
        // 初始化单位复根
        complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
        for(j=0; j<l; j+=h) // 控制起始下标
        {
            complex w(1,0); // 初始化螺旋因子
            for(k=j; k<j+h/2; k++) // 配对
            {
                u=y[k];
                t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn; // 更新螺旋因子
            } // 据说上面的操作叫蝴蝶操作…
        }
    }
    if(on==-1)    for(i=0; i<l; i++)    y[i].r/=l; // IDFT
}
char a[N/2],b[N/2];
int main()
{
    int i;
    while(~scanf("%s%s",a,b))
    {
        memset(sum,0,sizeof(sum));
        int l1=strlen(a),l2=strlen(b),l=1;
        while(l<l1*2||l<l2*2) l<<=1;
        for(i=0; i<l1; i++)
            x1[i].r=a[l1-i-1]-'0',x1[i].i=0;
        for(; i<l; i++)
            x1[i].r=x1[i].i=0;
        for(i=0; i<l2; i++)
            x2[i].r=b[l2-i-1]-'0',x2[i].i=0;
        for(; i<l; i++)
            x2[i].i=x2[i].r=0;
        fft(x1,l,1);
        fft(x2,l,1);
        for(i=0; i<l; i++) x1[i]=x1[i]*x2[i];
        fft(x1,l,-1);
        for(i=0; i<l; i++) sum[i]=x1[i].r+0.5;
        for(i=0; i<l; i++) sum[i+1]+=sum[i]/10,sum[i]%=10;
        i=l1+l2-1;
        while(sum[i]<=0&&i>0)
            i--;
        for(; i>=0; i--)
            printf("%d",sum[i]);
        puts("");
    }
    return 0;
}