A * B Problem Plus(fft)

题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=1402

hdu_1402:A * B Problem Plus

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 15419    Accepted Submission(s): 3047

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

 

Author

DOOM III

 

题解： 练习用fft实现大数的乘法达到O(nlog(n))的算法,两个数相乘看成是两个多项式的乘法，这样多项式中的x=10,fft套用模板即可

给出代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int MAX=300005;
const double PI=acos(-1.0),eps=1e-8;;
double cof1[MAX], cof2[MAX];
int n, k, permutation[MAX];
char s1[MAX],s2[MAX];
int ans[MAX];
struct complex {//复数
    double r, v;
    complex operator + (complex& obj) {
        complex temp;
        temp.r = r + obj.r;
        temp.v = v + obj.v;
        return temp;
    }
    complex operator - (complex& obj) {
        complex temp;
        temp.r = r - obj.r;
        temp.v= v - obj.v;
        return temp;
    }
    complex operator * ( complex& obj) {
        complex temp;
        temp.r = r*obj.r - v*obj.v;
        temp.v = r*obj.v + v*obj.r;
        return temp;
    }
} p1[MAX], p2[MAX], omiga[MAX], result1[MAX], result2[MAX];
void caculate_permutation(int s, int interval, int w, int next) {
    if(interval==n) {
        permutation[w] = s;
        return ;
    }
    caculate_permutation(s,interval*2, w, next/2);
    caculate_permutation(s+interval, interval*2, w+next, next/2);
}
void fft(complex transform[], complex p[]) {
    int i, j, l, num, m;
    complex temp1, temp2;
    for(i=0; i<n; i++)transform[i] = p[ permutation[i] ] ;
    num = 1, m = n;
    for(i=1; i<=k; i++) {
        for(j=0; j<n; j+=num*2)
            for(l=0; l<num; l++)
                temp2 = omiga[m*l]*transform[j+l+num],
                        temp1 = transform[j+l],
                                transform[j+l] = temp1 + temp2,
                                                 transform[j+l+num] = temp1 - temp2;
        num*=2,m/=2;
    }
}
void polynomial_by(int n1,int n2) {//多项式乘法，cof1、cof2保存的是a[0],a[1]..a[n-1]的值（a[i]*x^i）
    int i;
    double angle;
    k = 0, n = 1;
    while(n<n1+n2-1)k++,n*=2;
    for(i=0; i<n1; i++)p1[i].r = cof1[i], p1[i].v = 0;
    while(i<n)p1[i].r = p1[i].v = 0, i++;
    for(i=0; i<n2; i++)p2[i].r = cof2[i], p2[i].v = 0;
    while(i<n)p2[i].r = p2[i].v = 0, i++;
    caculate_permutation(0,1,0,n/2);
    angle = PI/n;
    for(i=0; i<n; i++)omiga[i].r = cos(angle*i), omiga[i].v = sin(angle*i);
    fft(result1,p1);
    fft(result2,p2);
    for(i=0; i<n; i++)result1[i]= result1[i]*result2[i];
    for(i=0; i<n; i++)omiga[i].v = -omiga[i].v;
    fft(result2, result1);
    for(i=0; i<n; i++)result2[i].r/=n;
    i = n -1;
    while(i&&fabs(result2[i].r)<eps)i--;
    n = i+1;
    while(i>=0) ans[i]=(int)(result2[i].r+0.5), i--;
}
int main() {
    while(scanf("%s",s1)!=EOF){
        scanf("%s",s2);
        int n1=strlen(s1),n2=strlen(s2);
        for(int i=0;i<n1;i++){
            cof1[i]=s1[n1-1-i]-'0';
        }
        for(int i=0;i<n2;i++){
            cof2[i]=s2[n2-1-i]-'0';
        }
        memset(ans,0,sizeof(ans));
        polynomial_by(n1,n2);
        for(int i=0;i<n;i++){
            if(ans[i]>=10){
                ans[i+1]+=ans[i]/10;
                ans[i]%=10;
            }
        }
        while(ans[n]>0){
            if(ans[n]>=10){
                ans[n+1]+=ans[n]/10;
                ans[n]%=10;
            }
            n++;
        }
        for(int i=n-1;i>=0;i--){
            putchar('0'+ans[i]);
        }
        puts("");
    }
    return 0;
}