【FFT】大整数乘法

来源：互联网发布：java二级考试题库下载编辑：程序博客网时间：2024/06/08 09:35

http://www.cnblogs.com/skyivben/archive/2008/07/23/1248413.html

整理一下模板

hdu1402

sincos需要手写，因为hdu没有

#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>#include <cmath>using namespace std;const double pi = acos(-1.0);const int maxn = 1 << 18;struct Complex {    double x, y;    Complex (double real = 0, double imag = 0) : x(real), y(imag) {}    double &real() {        return x;    }    double &imag() {        return y;    }        void print()    {        cout<<"real="<<x<<" imag=%.7lf\n"<<y<<endl;    }}Pa[500000],Pb[500000],Pc[500000];long long ans[500000],la,lb,N,A[500000],B[500000];char a[500000],b[500000];int t;Complex operator+(const Complex &a, const Complex &b) {    return Complex(a.x + b.x, a.y + b.y);}Complex operator-(const Complex &a, const Complex &b) {    return Complex(a.x - b.x, a.y - b.y);}Complex operator*(const Complex &a, const Complex &b) {    return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}inline void sincos(double theta,double &p0,double &p1){    p0=sin(theta);    p1=cos(theta);}void FFT(Complex P[], int n, int oper){    for (int i = 1, j = 0; i < n - 1; i++) {        for (int s = n; j ^= s >>= 1, ~j & s;);        if (i < j) {            swap(P[i], P[j]);        }    }    Complex unit_p0;    for (int d = 0; (1 << d) < n; d++) {        int m = 1 << d, m2 = m * 2;        double p0 = pi / m * oper;        sincos(p0, unit_p0.y, unit_p0.x);        for (int i = 0; i < n; i += m2) {            Complex unit = 1;            for (int j = 0; j < m; j++) {                Complex &P1 = P[i + j + m], &P2 = P[i + j];                Complex t = unit * P1;                P1 = P2 - t;                P2 = P2 + t;                unit = unit * unit_p0;            }        }    }}int main(){//    scanf("%d",&t);//    for (;t;t--) {  //      scanf("%s",a+1);//        scanf("%s",b+1);    for (;scanf("%s%s",a+1,b+1)==2;) {        la=strlen(a+1),lb=strlen(b+1);        for (int i=0;i<la;i++) A[i]=a[la-i]-'0';        for (int i=0;i<lb;i++) B[i]=b[lb-i]-'0';        for (N=1;N<=la+lb+1;N<<=1) ;        for (int i=0;i<la;i++) Pa[i]=Complex(A[i],0);        for (int i=0;i<N-la;i++) Pa[i+la]=Complex(0,0);        for (int i=0;i<lb;i++) Pb[i]=Complex(B[i],0);        for (int i=0;i<N-lb;i++) Pb[i+lb]=Complex(0,0);        FFT(Pa,N,1),FFT(Pb,N,1);        for (int i=0;i<N;i++) Pc[i]=Pa[i]*Pb[i];        FFT(Pc,N,-1);        for (int i=0;i<N+N;i++) ans[i]=0;        for (int i=0;i<N;i++) ans[i]=(long long)(Pc[i].x/N+0.5);        long long jw=0,len=0;        for (;len<N || jw;len++) {            long long x=(ans[len]+jw);            ans[len]=x%10;            jw=x/10;        }        for (;len>0 && !ans[len];len--) ;        for (int i=len;i>=0;i--) printf("%I64d",ans[i]);        printf("\n");    }    return 0;}

spoj MUL

#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>#include <cmath>using namespace std;const double pi = acos(-1.0);const int maxn = 1 << 18;struct Complex {    double x, y;    Complex (double real = 0, double imag = 0) : x(real), y(imag) {}    double &real() {        return x;    }    double &imag() {        return y;    }        void print()    {        cout<<"real="<<x<<" imag=%.7lf\n"<<y<<endl;    }}Pa[500000],Pb[500000],Pc[500000];int ans[500000],la,lb,N;char a[500000],b[500000];int t;Complex operator+(const Complex &a, const Complex &b) {    return Complex(a.x + b.x, a.y + b.y);}Complex operator-(const Complex &a, const Complex &b) {    return Complex(a.x - b.x, a.y - b.y);}Complex operator*(const Complex &a, const Complex &b) {    return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}void FFT(Complex P[], int n, int oper){    for (int i = 1, j = 0; i < n - 1; i++) {for (int s = n; j ^= s >>= 1, ~j & s;);if (i < j) {swap(P[i], P[j]);}}Complex unit_p0;for (int d = 0; (1 << d) < n; d++) {int m = 1 << d, m2 = m * 2;double p0 = pi / m * oper;sincos(p0, &unit_p0.y, &unit_p0.x);for (int i = 0; i < n; i += m2) {Complex unit = 1;for (int j = 0; j < m; j++) {Complex &P1 = P[i + j + m], &P2 = P[i + j];Complex t = unit * P1;P1 = P2 - t;P2 = P2 + t;unit = unit * unit_p0;}}}}int main(){    scanf("%d",&t);    for (;t;t--) {        scanf("%s",a+1);        scanf("%s",b+1);        la=strlen(a+1),lb=strlen(b+1);        for (N=1;N<=la+lb+1;N<<=1) ;        for (int i=0;i<la;i++) Pa[i]=Complex(a[la-i]-'0',0);        for (int i=0;i<N-la;i++) Pa[i+la]=Complex(0,0);        for (int i=0;i<lb;i++) Pb[i]=Complex(b[lb-i]-'0',0);        for (int i=0;i<N-lb;i++) Pb[i+lb]=Complex(0,0);        FFT(Pa,N,1),FFT(Pb,N,1);        for (int i=0;i<N;i++) Pc[i]=Pa[i]*Pb[i];        FFT(Pc,N,-1);        for (int i=0;i<N+N;i++) ans[i]=0;        for (int i=0;i<N;i++) ans[i]=(int)(Pc[i].x/N+0.5);        int jw=0,len=0;        for (;len<N || jw;len++) {            int x=(ans[len]+jw);            ans[len]=x%10;            jw=x/10;        }        for (;len>0 && !ans[len];len--) ;        for (int i=len;i>=0;i--) printf("%d",ans[i]);        printf("\n");    }    return 0;}

spoj VFMUL

需要压位来算，但是位数压多了之后就会导致精度不够，且会爆int

#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>#include <cmath>using namespace std;const long double pi = acos(-1.0);const int maxn = 1 << 18;struct Complex {    long double x, y;    inline Complex (long double real = 0, long double imag = 0) : x(real), y(imag) {}    long double &real() {        return x;    }    long double &imag() {        return y;    }        void print()    {        cout<<"real="<<x<<" imag="<<y<<endl;    }}Pa[250000],Pb[250000],Pc[250000];long long ans[250000],la,lb,N;char a[350000],b[350000];int t;inline Complex operator+(const Complex &a, const Complex &b) {    return Complex(a.x + b.x, a.y + b.y);}inline Complex operator-(const Complex &a, const Complex &b) {    return Complex(a.x - b.x, a.y - b.y);}inline Complex operator*(const Complex &a, const Complex &b) {    return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}inline void sincos(long double theta,long double &a,long double &b){a=sin(theta);b=cos(theta);}void FFT(Complex P[], int n, int oper){    for (int i = 1, j = 0; i < n - 1; i++) {for (int s = n; j ^= s >>= 1, ~j & s;);if (i < j) {swap(P[i], P[j]);}}Complex unit_p0;for (int d = 0; (1 << d) < n; d++) {int m = 1 << d, m2 = m * 2;long double p0 = pi / m * oper;//sincos(p0, &unit_p0.y, &unit_p0.x);sincos(p0,unit_p0.y,unit_p0.x);for (int i = 0; i < n; i += m2) {Complex unit = 1;for (int j = 0; j < m; j++) {Complex &P1 = P[i + j + m], &P2 = P[i + j];Complex t = unit * P1;P1 = P2 - t;P2 = P2 + t;unit = unit * unit_p0;}}}}int main(){    scanf("%d",&t);    for (;t;t--) {        scanf("%s",a+1);        scanf("%s",b+1);        la=strlen(a+1),lb=strlen(b+1);        for (N=1;N<=la/5+lb/5+3;N<<=1) ;        int cnt=0;        for (int i=0;i<la;i+=5,cnt++) {            long long x=0;            for (int j=max(1LL,la-(i+4));j<=la-i;j++) x=(x*10)+a[j]-'0';            Pa[cnt]=Complex(x,0);        }        for (int i=cnt;i<N;i++) Pa[i]=Complex(0,0);        cnt=0;        for (int i=0;i<lb;i+=5,cnt++) {            long long x=0;            for (int j=max(1LL,lb-(i+4));j<=lb-i;j++) x=(x*10)+b[j]-'0';            Pb[cnt]=Complex(x,0);        }        for (int i=cnt;i<N;i++) Pb[i]=Complex(0,0);        FFT(Pa,N,1),FFT(Pb,N,1);        for (int i=0;i<N;i++) Pc[i]=Pa[i]*Pb[i];        FFT(Pc,N,-1);        for (int i=0;i<N;i++) ans[i]=(long long)(Pc[i].x/N+0.5);        long long jw=0;int len=0;        for (;len<N;len++) {            long long x=(ans[len]+jw);            ans[len]=x%100000;            jw=x/100000;        }        for (;jw;len++) ans[len]=jw%100000,jw/=100000;        for (len--;len>0 && !ans[len];len--) ;        printf("%lld",ans[len]);          for (int i=len-1;i>=0;i--) printf("%05lld",ans[i]);        printf("\n");    }    return 0;}

0 0