蝴蝶操作

蝴蝶操作

 

FFT code:

#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<complex>
using namespace std;

typedef complex<double> cd;//复数类的定义
const int maxl=2094153;//nlogn的最大长度(来自leo学长的博客)
const double PI=3.14159265358979;//圆周率,不解释

cd a[maxl],b[maxl];//用于储存变换的中间结果
int rev[maxl];//用于储存二进制反转的结果
void getrev(int bit){
    for(int i=0;i<(1<<bit);i++){//高位决定二进制数的大小
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }//能保证(x>>1)<x,满足递推性质
}

void fft(cd* a,int n,int dft){//变换主要过程
    for(int i=0;i<n;i++){//按照二进制反转
        if(i<rev[i])//保证只把前面的数和后面的数交换,(否则数组会被翻回来)
            swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<n;step<<=1){//枚举步长的一半
        cd wn=exp(cd(0,dft*PI/step));//计算单位复根
        for(int j=0;j<n;j+=step<<1){//对于每一块
            cd wnk(1,0);//!!每一块都是一个独立序列,都是以零次方位为起始的
            for(int k=j;k<j+step;k++){//蝴蝶操作处理这一块
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;
                a[k+step]=x-y;
                wnk*=wn;//计算下一次的复根
            }
        }
    }
    if(dft==-1){//如果是反变换,则要将序列除以n
        for(int i=0;i<n;i++)
            a[i]/=n;
    }
}

int output[maxl];
char s1[maxl],s2[maxl];
int main(){
    scanf("%s%s",s1,s2);//读入两个数
    int l1=strlen(s1),l2=strlen(s2);//就算"次数界"
    int bit=1,s=2;//s表示分割之前整块的长度
    for(bit=1;(1<<bit)<l1+l2-1;bit++){
        s<<=1;//找到第一个二的整数次幂使得其可以容纳这两个数的乘积
    }
    for(int i=0;i<l1;i++){//第一个数装入a
        a[i]=(double)(s1[l1-i-1]-'0');
    }
    for(int i=0;i<l2;i++){//第二个数装入b
        b[i]=(double)(s2[l2-i-1]-'0');
    }
    getrev(bit);fft(a,s,1);fft(b,s,1);//dft
    for(int i=0;i<s;i++)a[i]*=b[i];//对应相乘
    fft(a,s,-1);//idft
    for(int i=0;i<s;i++){//还原成十进制数
        output[i]+=(int)(a[i].real()+0.5);//注意精度误差
        output[i+1]+=output[i]/10;
        output[i]%=10;
    }
    int i;  
    for(i=l1+l2;!output[i]&&i>=0;i--);//去掉前导零
    if(i==-1)printf("0");//特判长度为0的情况
    for(;i>=0;i--){//输出这个十进制数
        printf("%d",output[i]);
    }
    putchar('\n');
    return 0;
}

 

 NTT（快速数论变换）的多项式乘法的代码:

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
//#include<complex>
using namespace std;
//typedef complex<double> cd;

typedef long long LL;

void exgcd(int a,int b,int& x,int& y){
    if(b==0){
        x=1;
        y=0;
        return;
    }
    int x0,y0;
    exgcd(b,a%b,x0,y0);
    x=y0;y=x0-int(a/b)*y0;
}

int Inv(int a,int p){
    int x,y;
    exgcd(a,p,x,y);
    x%=p;
    while(x<0)x+=p;
    return x;
}

int qpow(int a,int b,int p){
    if(b<0){
        b=-b;
        a=Inv(a,p);
    }
    LL ans=1,mul=a%p;
    while(b){
        if(b&1)ans=ans*mul%p;
        mul=mul*mul%p;
        b>>=1;
    }
    return ans;
}

#define maxn (65537*2)
const int MOD=479*(1<<21)+1,G=3;

int rev[maxn];
void get_rev(int bit){
    for(int i=0;i<(1<<bit);i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
}

//from internet
//for(int i=0; i<NUM; i++)  
//{  
//    int t = 1 << i;  
//    wn[i] = quick_mod(G, (P - 1) / t, P);  
//}  

LL a[maxn],b[maxn];
void ntt(LL* a,int n,int dft){
    for(int i=0;i<n;i++){
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<n;step<<=1){
        LL wn;
        wn=qpow(G,dft*(MOD-1)/(step*2),MOD);
        for(int j=0;j<n;j+=step<<1){
            LL wnk=1;//这里一定要用long long不然会迷之溢出
            for(int k=j;k<j+step;k++){
                LL x=a[k]%MOD,y=(wnk*a[k+step])%MOD;//这里也要用long long
                a[k]=(x+y)%MOD;
                a[k+step]=((x-y)%MOD+MOD)%MOD;
                wnk=(wnk*wn)%MOD;
            }
        }
    }
    if(dft==-1){
        int nI=Inv(n,MOD);
        for(int i=0;i<n;i++)
            a[i]=a[i]*nI%MOD;
    }
}

#include<cstring>
char s1[maxn],s2[maxn];

int main(){
    //scanf("%*d");
    scanf("%s%s",s1,s2);
    int l1=strlen(s1),l2=strlen(s2);
    for(int i=0;i<l1;i++)a[i]=s1[l1-i-1]-'0';
    for(int i=0;i<l2;i++)b[i]=s2[l2-i-1]-'0';
    int bit,s=2;
    for(bit=1;(1<<bit)<(l1+l2-1);bit++)s<<=1;
    get_rev(bit);ntt(a,s,1);ntt(b,s,1);
    for(int i=0;i<s;i++)
        a[i]=a[i]*b[i]%MOD;
    ntt(a,s,-1);
    for(int i=0;i<s;i++){
        a[i+1]+=a[i]/10;
        a[i]%=10;
    }
    int cnt=s;
    while(cnt>=0 && a[cnt]==0)cnt--;
    if(cnt==-1)printf("0");
    for(int i=cnt;i>=0;i--){
        printf("%d",a[i]);
    }
    putchar('\n');
    return 0;
}

 

FFT详解

转载于:https://www.cnblogs.com/guxuanqing/p/9502631.html