2018-01-20 HDU 1402 Ver.A FFT 快速傅里叶变换

题目描述:

Calculate A * B.

输入:

Each line will contain two integers A and B. Process to end of file.

Note: the length of each integer will not exceed 50000.

输出:

For each case, output A * B in one line.

样例输入:

1

2

1000

2

样例输出:

2

2000

题目描述很简单的, 就是给你两个数乘起来就成, 但是, 长度五万哦

想都不用想, 肯定是是字符串做啦

如果直接暴力模拟乘法和进位那肯定是要超时的, 但是究其原理很像多项式相乘, 而多项式相乘可以用快速傅里叶变换优化.

举个例子

计算162*25

先写成这个样子

(1*x^2+6*x+2)*(2*x+5) , 其中x=10带进去就是原式了, 这不难理解, 有一点生成函数的味道(但是并不是啊!),

好,换个写法 :  (  2* x^0  +  6* x^1  +  1* x^2  ) * (  5* x^0  +  2* x^1  ) 里面的顺序就是个 十 百

拆开 化简 那就是 10*x^0 + 34*x^1 + 17*x^2 + 2*x^3  顺序依然是个 十 百 千

但是我们是十进制啊, 哪会有说哪一位是 10 34 17的呢?

还记不记得把x=10带进去就是结果这件事情? 我们把系数上的十位进上去, 然后把个位留下来,试试看?

10* x^0 + 34* x^1 + 17* x^2 + 2* x^3

0* x^0 + (34+10/10)* x^1 + 17* x^2 + 2* x^3

0* x^0 + 35* x^1 + 17* x^2 + 2* x^3

0* x^0 + 5* x^1 + (17+30/10)* x^2 + 2* x^3

0* x^0 + 5* x^1 + 20* x^2 + 2* x^3

0* x^0 + 5* x^1 + 0* x^2 + (2+20/10)* x^3

0* x^0 + 5* x^1 + 0* x^2 + 4* x^3

根据个十百千的顺序, 得到结果就是4050, 笔算或者按按计算器发现162*25=4050.

所以我们发现, 我们可以用多项式相乘来模拟大位数乘法, 乘完之后可能并不是我们想要的结果, 但是取余进位就可以处理完成, 并且这个处理取余的时间并不是很长, 所以关键是如何更加快速地进行这个多项式的相乘.

这里先介绍FFT即快速傅里叶变换, 使用方法大致如下.

先把多项式中的系数放入复数数组的实部之中, 虚部置0

像母函数一样, 比如162放进去就是(.a指实部, .b指虚部)

a[0].a=2; a[0].b=0

a[1].a=6; a[1].b=0;

a[2].a=1; a[2].b=0;

a[3].a=0; a[3].b=0; 行了 后面不管实部虚部都是0

所以读入字符串后放进复数数组之后记得倒置.

第二个数字b同理.

接下来就是寻找fft函数的另一个参数了----len( 不一定是len, 随便取了个变量名而已)

怎么找?

一个多项式和另一个多项式分别取最高次数然后求和, 然后找最小的但是大于这个和的2次幂, 参见代码那个while循环和注释.

然后分别输入那两个复数数组, len, 进行正向的变换, 转换成(我自称的)"fft"形式

有趣的事情发生了!fft形式的多项式可以直接线性相乘!

所以 c[ i ] = a[ i ] * b[ i ] ( 复数数组c是存a*b的结果的, 乘完也是fft形式, 并且这里的乘是重载的复数相乘哦 )

得到了结果多项式的fft形式!

所以逆变换一下, 参数仍然是转换的复数数组, len, 但是后面的on参数为-1, 就是逆变换的意思

逆变换完了!但是由于一些原因, c[i].a并不是整数, 解决方法就是四舍五入.

最后就是刚刚详细举例的取余进位啦, 最后输出的时候注意一下, 比如数字8和字符'8'的区别, 不过这个是C语言基础知识了, 不必多说.

代码如下:

#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
using namespace std;
const double PI=acos(-1);

struct complex
{
    double a,b;//a表示实部,b表示虚部
    complex(double aa=0.0,double bb=0.0)
    {
        a=aa;
        b=bb;
    }
    complex operator +(const complex &e)
    {
        return complex(a+e.a,b+e.b);
    }
    complex operator -(const complex &e)
    {
        return complex(a-e.a,b-e.b);
    }
    complex operator *(const complex &e)
    {
        return complex(a*e.a-b*e.b,a*e.b+b*e.a);
    }
};

void change(complex *y,int len)
{
    int i,j,k;
    for(i=1,j=len/2;i<len-1; i++)
    {
        if(i<j)
            swap(y[i],y[j]);
        k=len/2;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        if(j<k)
        {
            j+=k;
        }
    }
}
void fft(complex *y,int len,int on)
{
    change(y, len);
    int i,j,k;
    for(i=2;i<=len;i<<=1)
    {
        complex wn(cos(-on*2*PI/i),sin(-on*2*PI/i));
        for(j=0;j<len;j+=i)
        {
            complex w(1,0);
            for(k=j;k<j+i/2;k++)
            {
                complex u=y[k],t=w*y[k+i/2];
                y[k]=u+t;
                y[k+i/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)
        for(i=0;i<len;i++)
            y[i].a/=len;
}

char A[200005],B[200005];
complex a[200005],b[200005],c[200005];
long long int res[200005];
int len_a,len_b,len_c,len;
//怕函数体里面塞不下所以.....

int main()
{
    int i;
    while(scanf("%s%s",A,B)!=EOF)
    {
        len=1;
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        memset(c,0,sizeof(c));//置零
        memset(res,0,sizeof(res));
        len_a=(int)strlen(A);
        len_b=(int)strlen(B);
        for(i=0;i<len_a;i++)
            a[len_a-1-i].a=A[i]-'0';//倒置,fft优化的多项式存储时,下标为k则代表多项式中对应的项次数为k,所以倒置
        for(i=0;i<len_b;i++)
            b[len_b-1-i].a=B[i]-'0';//倒置
        while(len<(len_a+len_b))
            len=len<<1;//len为大于结果多项式次数的2次幂,比如最后结果最高次项为10,那么取len为2^4=16
        fft(a,len,1);//注意后面on的参数,和多项式长度参数传的是什么
        fft(b,len,1);//on=1,表示从原始的多项式转化成便于处理的"fft"形式
        for(i=0;i<len;i++)//接下来开始线性处理
            c[i]=a[i]*b[i];//这里是虚数的乘法哦,在结构体声明里面重载过
        fft(c,len,-1);//on==-1,就是把"fft"形式转化成普遍的多项式形式
        for(i=0;i<len;i++)
            res[i]=(int)(c[i].a+0.5);//转化过来的普遍的多项式形式是小数,用四舍五入的方法约成整数
        for(i=0;i<=len;i++)
        {
            res[i+1]+=res[i]*0.1;
            res[i]=res[i]%10;//向高位进位
        }
        len_c=0;
        for(i=len-1;i>=0;i--)
        {
            if(res[i]!=0)
            {
                len_c=i;
                break;
            }
        }
        for(i=len_c;i>=0;i--)
        {
            putchar((int)res[i]+'0');
        }
        puts("");
    }
    return 0;
}