51nod1027 大数乘法（FFT）

题目：

Input

第1行：大数A第2行：大数B(A,B的长度 <= 1000，A,B >= 0）

Output

输出A * B

O(nlgn)的快速傅里叶算法如下：

#include <iostream>#include <stdio.h>#include <cmath>#include <cstring>using namespace std;#define N 1050*2const double PI = acos(-1.0);struct Vir //复数类型{double re, im;Vir(double _re = 0., double _im = 0.) : re(_re), im(_im) {} //初始化构造函数/*重载运算符*/Vir operator*(Vir r) { return Vir(re*r.re - im*r.im, re*r.im + im*r.re); } //复数乘法Vir operator+(Vir r) { return Vir(re + r.re, im + r.im); } //复数加法Vir operator-(Vir r) { return Vir(re - r.re, im - r.im); } //复数减法};char a[N * 2], b[N * 2];Vir pa[N * 2], pb[N * 2];int ans[N * 2];void bit_rev(Vir *a, int loglen, int len){for (int i = 0; i < len; ++i){int t = i, p = 0;for (int j = 0; j < loglen; ++j) //p是t的位反转结果{p <<= 1;p = p | (t & 1);t >>= 1;}if (p < i) //下标反转后变小，则交换a[p]与a[i]{Vir temp = a[p];a[p] = a[i];a[i] = temp;}}}void FFT(Vir *a, int loglen, int len, int on) //len = 2^loglen，且为向量a的长度{bit_rev(a, loglen, len); //位逆序置换for (int s = 1, m = 2; s <= loglen; ++s, m <<= 1) //s表示二叉树的每一层，m表示该层节点序列长度（m = 2^s）{Vir wn = Vir(cos(2 * PI*on / m), sin(2 * PI*on / m)); //wn = e ^2πi/m or e ^-2πi/mfor (int i = 0; i<len; i += m) //该层有len/m个节点序列{Vir w = Vir(1.0, 0);for (int j = 0; j<m/2; ++j) //每个节点序列有m/2次蝴蝶操作{Vir v = w*a[i + j + m / 2];Vir u = a[i + j];a[i + j] = u + v;a[i + j + m / 2] = u - v;w = w*wn;}}}if (on == -1) //逆DFT, 计算结果的每个元素除以len{for (int i = 0;i<len;++i) a[i].re /= len, a[i].im /= len;}}int main(){while (scanf("%s%s", a, b) != EOF){int lena = strlen(a);int lenb = strlen(b);int n = 1, loglen = 0;while (n < lena + lenb) //n为大于或等于lena+lenb的最小的2的幂次，幂次值为loglenn <<= 1, loglen++;/*多项式系数转换为复数（虚部为0，实部为系数），并将a,b的次数界都扩展为n*/for (int i = 0, j = lena - 1; i<n; ++i,--j)pa[i] = Vir(j >= 0 ? a[j] - '0' : 0. ,  0.);for (int i = 0, j = lenb - 1; i<n; ++i, --j)pb[i] = Vir(j >= 0 ? b[j] - '0' : 0. ,  0.);for (int i = 0;i <= n;++i)ans[i] = 0;FFT(pa, loglen, n, 1);FFT(pb, loglen, n, 1);for (int i = 0; i<n; ++i) //θ(n)完成点值对乘法pa[i] = pa[i] * pb[i];FFT(pa, loglen, n, -1);for (int i = 0; i<n; ++i) ans[i] = pa[i].re + 0.5;for (int i = 0; i<n; ++i) //从低位向高位进位ans[i + 1] += ans[i] / 10, ans[i] %= 10;int pos = lena + lenb - 1;for (; pos>0 && ans[pos] <= 0; --pos); //去掉高位的前导0for (; pos >= 0; --pos) printf("%d", ans[pos]);puts("");}return 0;}

大数乘法V2

我觉得我用的是C语言。结果提交上去就是CE...。以C++提交就可以，醉了*****

#include <stdio.h>#include <math.h>#include <string.h>#define PI acos(-1.0)#define N 100005 * 2typedef struct complex{double real;double image;}Complex;void initComplex(Complex* c, double _real, double _image) // init complex = _real + _image*i{c->real = _real, c->image = _image;}Complex AddComplex(Complex a, Complex b) //return complex s = complex a + complex b{Complex s;s.real = a.real + b.real;s.image = a.image + b.image;return s;}Complex subComplex(Complex a, Complex b) //return complex v = complex a - complex b{Complex v;v.real = a.real - b.real;v.image = a.image - b.image;return v;}Complex mulComplex(Complex a, complex b) //return complex p = complex a * complex b{complex p;p.real = a.real*b.real - a.image*b.image;p.image = a.real*b.image + a.image*b.real;return p;}void exchange(Complex* a, Complex* b) //exchange value{Complex  temp;temp = *a;*a = *b;*b = temp;}void bitReserve(Complex a[], int n, int loglen) //let a[0 .. n-1]'s index reserve, 0<=bitlen(index)<=loglen{int i, j, t, i_rev;for (i = 0; i < n - 1; i++){t = i, i_rev = 0;for (j = 0; j < loglen; j++) //each i -> i_rev{i_rev <<= 1;i_rev |= (t & 1);t >>= 1;}if (i < i_rev){exchange(&a[i], &a[i_rev]);}}}void FFT(Complex a[], int n, int loglen, int on) //get a's DFT when on  is 1, get a's DFT^-1 when on is -1{int i, j, k, m;Complex wn, w, t, u;bitReserve(a, n, loglen);m = 2; //each floor node's lengthfor (i = 1; i <= loglen; i++) //represent recursive tree's each floor{initComplex(&wn, cos(2*PI*on/m), sin(2*PI*on/m)); //e ^ 2πi/m*onfor (j = 0; j < n; j += m) //every floor has n/m nodes{initComplex(&w, 1.0, 0.0); //w = 1for (k = 0; k < m / 2; k++) //every node has m/2 butterfly operation{t = mulComplex(w, a[j + k + m / 2]); //common subexpressionu = a[j + k];a[j + k] = AddComplex(u, t);a[j + k + m / 2] = subComplex(u, t);w = mulComplex(w, wn);}}m <<= 1;}if (on == -1) //DFT^-1{for (i = 0; i < n; i++){a[i].real /= n;a[i].image /= n;}}}void convert(char str[], Complex c[], int n){int i, j, len;len = strlen(str);for (i = 0, j = len - 1; i < n; i++, j--){if (j >= 0)initComplex(&c[i], str[j] - '0', 0.0);elseinitComplex(&c[i], 0.0, 0.0);}}void multiplyBNum(char a[], char b[]) // return product of bignum a with bignum b{Complex complex_a[N * 2], complex_b[N * 2], complex_p[N * 2];int ans[N * 2];int alen, blen, n, loglen = 0, i;alen = strlen(a), blen = strlen(b);for (n = 1; n < alen + blen; n <<= 1) {//n’min >= alen+blen and n = 2^loglenloglen++;}convert(a, complex_a, n);convert(b, complex_b, n);for (i = 0; i < n; i++)ans[i] = 0;FFT(complex_a, n, loglen, 1);FFT(complex_b, n, loglen, 1);for (i = 0; i < n; i++)complex_p[i] = mulComplex(complex_a[i], complex_b[i]);FFT(complex_p, n, loglen, -1);for (i = 0; i < n; i++)ans[i] = (int)(complex_p[i].real + 0.5);for (i = 0; i < n; i++){ans[i + 1] += ans[i] / 10;ans[i] %= 10;}for (i = alen + blen - 1; i > 0 && ans[i] <= 0; i--);for (; i >= 0; i--)printf("%d", ans[i]);printf("\n");}int main(){char x[N], y[N];scanf("%s", x);scanf("%s", y);multiplyBNum(x, y);return 0;}