[BZOJ2179] FFT快速傅立叶&高精度乘法

「BZOJ2179」

Code

---

#include<bits/stdc++.h>
using namespace std;
#define pi acos(-1)
 
typedef complex<double> C;
const int N = 131100;
char s[N],t[N];
int n,m,l,r[N],c[N];
C a[N],b[N];
 
void fft(C *a, int f) {
    for(int i = 0; i < n; i++) if(r[i] > i) swap(a[i], a[r[i]]);//分配正确的顺序
    for(int i = 1; i <n; i <<= 1) {//i是当前要计算的区间的长度的一半
        C wn(cos(pi/i), f*sin(pi/i));//由于需要大区间的单位根，所以2*pi/2*i，消去2
        for(int j = 0; j < n; j += i<<1) {//j是每两个小区间的起点
            C w = 1;
            for(int k = 0; k < i; k++, w *= wn) {//在小区间内枚举复数，一次合并两个小区间
                C x = a[j+k], y = w*a[j+k+i];//蝴蝶变换
                a[j+k] = x+y, a[j+k+i] = x-y;
            }
        }
    }
}
int main() {
    scanf("%d%s%s", &m, s, t);
    for(int i = 0; i < m; i++) a[i] = s[m-i-1]-'0', b[i] = t[m-i-1]-'0';
    for(n = 1, m <<= 1; n < m; n <<= 1) l++;
    for(int i = 0; i < n; i++) r[i] = (r[i>>1]>>1)|((i&1)<<(l-1));//i到i*2要<<1,则r[i]到r[i*2]就要>>1,由于是从r[i]推到r[i*2],i*2第最低位的信息会缺失，所以要|((i&1)<<(l-1))
    fft(a, 1), fft(b, 1);//通过离散傅里叶变换（把n个复数带入多项式，得到点值表示）求出a,b的点值表示
    for(int i = 0; i < n; i++) a[i] *= b[i];//新的多项式的点至表示就是g(x){(x0,a(x0)*b(x0)));(x1,a(x1)*b(x1));...(xn-1,a(xn-1)*b(xn-1))}
    fft(a, -1);//把前面的点值表示作为新多项式的系数，并把新多项式转化成系数表示法
    for(int i = 0; i < n; i++) a[i] /= n;//以下就是输出答案了
    for(int i = 0; i < m; i++) c[i] = (int)(a[i].real()+0.1);
    for(int i = 0; i < m; i++) if(c[i] >= 10) {
        c[i+1] += c[i]/10, c[i] %= 10;//模拟每一位进位过程，画个图就明白了
    } else if(!c[i] && i == m-1) m--;//如果最高为是零，就m--
    for(int i = m-1; ~i; i--) printf("%d", c[i]);//倒序输出
    return 0;
}

• 第一道BZOJ权限题纪念