FFT - A * B Problem Plus - HDU 1402

FFT - A * B Problem Plus - HDU 1402

Problem Description

Calculate A * B.

Input

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

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

Output

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

Sample Input

1
2
1000
2

Sample Output

2
2000

---

分析：

$FFT介绍可参考刘汝佳白书以及OIWIKI。$

简单总结一些规律和性质：

$①、如何理解卷积：多项式乘法实际上是多项式系数向量的卷积。$

$$\begin{gathered} ②、FFT可以加速多项式在系数表示和点值表示之间的转换。 \\ 多项式的系数表示法对应于时域，而点值表示法对应于频域。 \end{gathered}$$

$③、点值表示非常适合做乘法：只要两个多项式的点集${$x_i$}$相同，则只需要把对应的值相乘即可。$

$比如，C(x)=A(x)\cdot B(x)，则C(x_i)=A(x_i)\cdot B(x_i)$

FFT具体步骤：

$$\begin{gathered} ①、补0：在两个多项式的最前面补0，补到两个多项式的最高次项为2n次，n=max(n1,n2)。 \\ 设系数向量分别为v_1和v_2。 \end{gathered}$$

$②、求值：用FFT计算f_1=DFT(v_1)，f_2=DFT(v_2)。f_1和f_2分别是两个多项式在各单位根处的取值。$

$③、乘法：把两个向量f_1和f_2的每一维相乘，得到向量f，它对应输入多项式乘积的点值表示。$

$④、插值：用FFT计算v=IDFT(f)，v就是乘积的系数向量。$

$代码是kuang神的模板。$

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<iostream>
#include<cmath>

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);
    }
};

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 -= k;
            k /= 2;
        }
        if(j<k) j += k;
    }
}
/*
 * 做fft
 * len必须为2^k形式，
 * 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(-on*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;
            }
        }
    }
    if(on == -1)
        for(int i=0;i<len;i++)
            y[i].x /= len;
}

const int N = 2e5+10;

Complex x1[N], x2[N];
char str1[N], str2[N];
int sum[N];

int main()
{
    while(~scanf("%s%s",str1,str2))
    {
        int len1 = strlen(str1), 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];
        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--;  //前导0
        for(int i=len;i>=0;i--)             //倒序输出
            printf("%c",sum[i]+'0');
        puts("");
    }
    return 0;
}