FFT

FFT计算多项式乘法与大数乘法

多项式乘法

#include <iostream>
#include<stdio.h>
#include <complex> //复数模板头文件
const int maxN=1e6;
const double pi=3.1415926;
using namespace std;
void FFT(complex<double> *a,int n,int inv);
complex<double> a[maxN],b[maxN];
int n,n1,n2,cnt,rev[maxN],c[maxN],sign=1;//合成函数
int main()
{
    scanf("%d%d",&n1,&n2);
    n=1;
    while(n<n1+n2)n<<=1;
    for(int i=0;i<n1;i++){
        scanf("%lf",&a[i].real());
    }
    for(int i=0;i<n2;i++){
        scanf("%lf",&b[i].real());
    }
    FFT(a,n,1);
    FFT(b,n,1);
    for(int i=0;i<=n;i++){
        a[i]=a[i]*b[i];
    }
    FFT(a,n,-1);
    for(int i=0;i<m;i++){
        printf("%d ",(int)(a[i].real()/n+0.5));
    }
    return 0;
}
void FFT(complex<double> *a,int n,int inv){
    int bit=0;
    while((1<<bit)<n)bit++;
    for(int i=0;i<n;i++){ //为每一个值作下标的颠倒
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
        if(i<rev[i])//防止又转换回来
            swap(a[i],a[rev[i]]);
    }//下标颠倒之后a[i]代表基本值，len为2时加并一步步叠加
    for(int mid=1;mid<n;mid<<=1){ //从最底层开始叠加
        complex<double> wn(cos(pi/mid),inv*sin(pi/mid));//每一层的单位根
        for(int i=0;i<n;i+=mid<<1){//给每一层分组，i是每组起始位置
            complex<double> w(1,0);//
            for(int j=0;j<mid;j++){//组内运算
                complex<double> x=a[i+j],y=w*a[i+j+mid];
                a[i+j]=x+y;a[i+j+mid]=x-y;
                w*=wn;//运算时取实数部分
            }
        }
    }
}

大数乘法

数据
函数read将字符串s1,s2也就是大数转化为复数的实数部分；经过FFT转化化为c的系数，最后每位相加，注意进位；

#include <iostream>
#include <stdio.h>
#include <string.h>
#include <complex> //复数模板头文件
using namespace std;
const int maxN=1001000;
const double pi=acos(-1.0);
void FFT(complex<double> *a,int n,int inv);
int read(char s[],complex<double> *a);
int n,n1,n2,cnt,rev[maxN],c[maxN],sign=1;//合成函数
void trans(int c[],int n);
void trans_1(int c[],int n);
complex<double> a[maxN],b[maxN];
int main()
{
    char s1[maxN],s2[maxN];
    while(~scanf("%s%s",s1,s2)){
        sign=1;
        n=1;
        n1=read(s1,a);
        n2=read(s2,b);
        if(sign==-1)
            printf("-");
        while(n<n1+n2)n<<=1;//将n化作c的系数值,n为2的倍数;
        FFT(a,n,1);//系数化为点值表示
        FFT(b,n,1);
        for(int i=0;i<n;i++){
            a[i]=a[i]*b[i];
        }
        FFT(a,n,-1);//点值化为系数表示
        trans(c,n);
        memset(a,0,(n+10)*sizeof(complex<double>));
        memset(b,0,(n+10)*sizeof(complex<double>));
        memset(c,0,(n+10)*sizeof(int));
        //memset(rev,0,(n+10)*sizeof(int));
    }
    return 0;
}
void trans(int c[],int n){
    for(int i=0;i<n;i++){
        c[i]+=(int)(a[i].real()/n+0.5);//将c的系数存放在c中
        if(c[i]>=10)
            c[i+1]+=c[i]/10,c[i]%=10,n+=(i==n);
    }
    while(!c[n]&&n>0){
        n--;
    }
    for(int i=n;i>=0;i--){
        printf("%d",c[i]);
    }
    printf("\n");
}
void trans_1(int c[],int n){
    for(int i=0;i<n;i++){
        c[i]=a[i].real()/n+0.5;//将c的系数存放在c中
    }
    int bit_add=0,k=0;
    for(int i=0;i<n;i++){
        int temp=c[i]+bit_add;
        if(temp>=10){
            c[k++]=temp%10;
            bit_add=temp/10;
        }
        else{
            c[k++]=temp;
            bit_add=0;
        }
    }
    while(bit_add){
        c[k++]=bit_add;
        bit_add/=10;
    }
    k--;
    while(c[k]==0&&k>0)
        k--;
    for(int i=k;i>=0;i--){
        printf("%d",c[i]);
    }
    printf("\n");
}
void FFT(complex<double> *a,int n,int inv){
    int bit=0;
    while((1<<bit)<n)bit++;
    for(int i=0;i<n;i++){ //为每一个值作下标的颠倒
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
        if(i<rev[i])//防止又转换回来
            swap(a[i],a[rev[i]]);
    }//下标颠倒之后a[i]代表基本值，len为2时加并一步步叠加
    for(int mid=1;mid<n;mid<<=1){ //从最底层开始叠加
        complex<double> wn(cos(pi/mid),inv*sin(pi/mid));//每一层的单位根
        for(int i=0;i<n;i+=mid<<1){//给每一层分组，i是每组起始位置
            complex<double> w(1,0);//
            for(int j=0;j<mid;j++){//组内运算
                complex<double> x=a[i+j],y=w*a[i+j+mid];
                a[i+j]=x+y;a[i+j+mid]=x-y;
                w*=wn;//运算时取实数部分
            }
        }
    }
}
int read(char s[],complex<double> *a){
    int n=0;
    for(int i=strlen(s)-1;i>=0;i--){
        if(s[i]=='-'){
            if(sign==1)
                sign=-1;
            else
                sign=1;
        }
        else if(s[i]>='0'&&s[i]<='9'){
            complex<double> m(s[i]-'0',0);
            a[n++]=m;
        }
    }
    return n;
}