快速傅里叶变换(FFT)

作用：加速多项式乘法

朴素高精度乘法时间O(n^2),但FFT能O(nlog2n)的时间解决

前置知识：

1.点值表示法：

f(x)={( x0,f(x0) ),( x1,f(x1) ) ,( x2, f(x2) ), ( x3, f(x3) ), ( x4, f(x4) ), ... , (xn-1, f(xn-1) )}

g(x)={( x0,g(x0) ),( x1,g(x1) ) ,( x2, g(x2) ), ( x3, g(x3) ), ( x4, g(x4) ), ... , (xn-1, g(xn-1) )}

设它们乘积为h(x)，那么

h(x)={( x0,f(x0)g(x0) ),( x1, f(x1)g(x1) ), ( x2, f(x2)g(x2) ), ... , ( xn-1, f(xn-1)g(xn-1) )}

2.复数

(a1,θ1) *（a2,θ2）为（a1a2,θ1+θ2）

快速傅里叶变换的实现：

const double PI = acos(-1);
typedef complex <double> cp;

cp omega(int n, int k) {
    return cp(cos(2 * PI * k / n), sin(2 * PI * k / n));
}

void fft(cp *a, int n, bool inv) {
    if(n == 1) return ;
    static cp buf[N];
    int m = n / 2;
    for(int i = 0; i < m ; i++){        
        buf[i] = a[2 * i];
        buf[i + m] = a[2 * i + 1];
    }
    for(int i = 0; i < n; i++)
        a[i] = buf[i];
    fft(a, m, inv);
    fft(a + m, m, inv);
    for(int i = 0; i < m; i++) {
        cp x = omega(n, i);
        if(inv) x = conj(x);
        //conj是一个自带的求共轭复数的函数，精度较高。当复数模为1时，共轭复数等于倒数
        buf[i] = a[i] + x * a[i + m];
        buf[i + m] = a[i] - x * a[i + m];
    }
    for(int i = 0; i < n; i++)
        a[i] = buf[i];
}

注1：a中i属于0-m-1是A1(i)的取值，a中i属于m到n-1是A2(i)的取值，那么递归的分析，对现在A1中的第i位置，其是由0-m/2-1的得来的，然后最后再更新a数组，a(i）就是A（omega(n, i)）的取值了。

注2：n是偶数，多项式的项数。

注3：inv表示单位根是否要取倒数，FFT的逆变换即点值表示法转化为系数表示法。

做法：把点值表示法作为系数，用取了倒数的单位根代入求个点值表示法，得到Zi再除以n就是i的系数ai（证明参考：https://www.cnblogs.com/RabbitHu/p/FFT.html）以下提到"博客"均指此篇博客。

注4：在逆FFT的时候，应该用floor(a[i].real() / n + 0.5);来得到res[i]（精度问题）

优化fft(非递归)

发现每次都往下其实都是先去递归，使得所有元素到达其应该在的地方，然后再不断往上递归对a赋值的。

规律：参考博客，可以发现，一个位置a上的数，最后所在的位置是"a二进制翻转得到的数"

据此写出非递归版本fft：先把每个数放到最后的位置上，然后不断向上，从而求出最终答案。

#include<iostream>
#include<cstdio>
#include<complex>
using namespace std;
const double PI = acos(-1);
typedef complex <double> cp;

cp a[N], b[N], omg[N], inv[N];

void init() {
    for(int i = 0; i < n; i++) {
        omg[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n));
        inv[i] = conj(omg[i]);
    }
}

void fft(cp *a, cp *omg) {
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(int i = 0; i < n; i++) {
        int t = 0;
        for(int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        if(i < t) swap(a[i], a[t]);    //i < t的限制使得每对点只被交换一次(否则交换两次相当于没交换) 
    }
    static cp buf[N];
    for(int l = 2; l <= n; l *= 2) {    //区间长度 
        int m = l / 2;
        for(int j = 0; j < n; j += l)    //区间起点 
            for(int i = 0; i < m; i++) {
                buf[j + i] = a[j + i] + omg[n / l * i] * a[j + i + m];
                buf[j + i + m] = a[j + i] - omg[n / l * i] * a[j + i + m];
            }
        for(int j = 0; j < n; j++)
            a[j] = buf[j];
    }
}

蝴蝶变换：

之前为什么需要buf数组：

如果

a[j + i] = a[j + i] + omg[n / l * i] * a[j + i + m]

a[j + i + m] = a[j + i] - omg[n / l * i] * a[j + i + m]

会对更新a[j + i + m]造成影响。

而通过蝴蝶变换

cp t = omg[n / l * i] * a[j + i + m]

a[j + i + m] = a[j + i] - t

a[j + i] = a[j + i] + t

不就顺序换了一下？？

反正就不用buf数组就是了。

FFT最终模板：

const double PI = acos(-1);
typedef complex <double> cp;

cp a[N], b[N], omg[N], inv[N];

void init() {
    for(int i = 0; i < n; i++) {
        omg[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n));
        inv[i] = conj(omg[i]);
    }
}

void fft(cp *a, cp *omg) {
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(int i = 0; i < n; i++) {
        int t = 0;
        for(int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        if(i < t) swap(a[i], a[t]);    //i < t的限制使得每对点只被交换一次(否则交换两次相当于没交换) 
    }
    for(int l = 2; l <= n; l *= 2) {    //区间长度 
        int m = l / 2;
        for(cp *p = a; p != a + n; p += l)
            for(int i = 0; i < m; i++) {
                cp t = omg[n / l * i] * p[i + m];
                p[i + m] = p[i] - t;
                p[i] += t;
            }
    }
}

题1：a*bIII http://www.acmicpc.sdnu.edu.cn/problem/show/1531

注意三点：

1.0与其他数相乘只输出一个0.

2.memset是可以对复数初始化的

3.FFT中的n需要是2的倍数

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1000005;
const double PI = acos(-1);
typedef complex <double> cp;

char sa[N], sb[N];
int n = 1, lena, lenb, ans[N];
cp a[N], b[N], omg[N], inv[N];

void init(){
    for(int i = 0; i < n; i++){
        omg[i] = cp(cos(2 * PI * i / n), sin(2 * PI * i / n));
        inv[i] = conj(omg[i]);
    }
}

void fft(cp *a, cp *omg) {
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(int i = 0; i < n; i++) {
        int t = 0;
        for(int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        if(i < t) swap(a[i], a[t]);    //i < t的限制使得每对点只被交换一次(否则交换两次相当于没交换) 
    }
    for(int l = 2; l <= n; l *= 2) {    //区间长度 
        int m = l / 2;
        for(cp *p = a; p != a + n; p += l)
            for(int i = 0; i < m; i++) {
                cp t = omg[n / l * i] * p[i + m];
                p[i + m] = p[i] - t;
                p[i] += t;
            }
    }
}


signed main()
{
    while(~scanf("%s%s",sa,sb)) {
        memset(ans, 0, sizeof(ans));
        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));
        n = 1;
        lena = strlen(sa), lenb = strlen(sb);
        if(lena == 1  && sa[0] == '0' || lenb == 1  && sb[0] == '0') {
            printf("0\n");
            continue;
        }
        while(n < lena + lenb) n *= 2;
        for(int i = 0; i < lena; i++){
            a[i].real(sa[lena - 1 - i] - '0');
        }
        for(int i = 0; i < lenb; i++){
            b[i].real(sb[lenb - 1 - i] - '0');
        }
        init();
        fft(a, omg);
        fft(b, omg);
        for(int i = 0; i < n; i++) {
            a[i] *= b[i];
        }
        fft(a, inv);
        for(int i = 0; i < n; i++) {
            ans[i] += floor(a[i].real() / n + 0.5);
            ans[i + 1] += ans[i] / 10;
            ans[i] %= 10; 
        }
        int beg;
        for(int i = n-1; i >= 0; i--) {
            if(ans[i] != 0){
                beg = i;
                break;
            }
        }
        for(int i = beg; i >= 0; i--){
            printf("%lld",ans[i]);
        }
        putchar('\n');
    }
    return 0;
}

FFT的缺点是它的复数运算double精度问题导致它实际上是k*nlongn的，会比NTT的常数大很多。

NTT

前置知识：

原根：对于g,p属于Z， 如果g^i mod p ( 1<=i<=p-1)的值互不相同，则称g为p的原根。

或者说对于任意i，j(1<=i<j <= p-1) g^i mod p /= g^j mod p，那么g为p的原根。

常见模数有：998244353,1004535809,469762049，这几个的原根都为3

在NTT中，我们拿原根来代替FFT的单位根

#define g 3
const int mod = 998244353;
const int N = 300000;

inline get_rev()
{
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1 | ((i & 1) << (lim - 1)));
    }
}

inline void ntt(int *a, int inv) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) swap(a[i], a[rev[i]]);    //i < t的限制使得每对点只被交换一次(否则交换两次相当于没交换) 
    }
    for(int l = 2; l <= n; l *= 2) {    //区间长度 
        int m = l / 2;
        int tmp = q_pow(g, (mod-1)/l);
        if(inv == -1) tmp = q_pow(tmp, mod-2);
        for(int i = 0; i < n; i += l) {
            int omega = 1;
            for(int j = 0; j < m; j++, omega = omega*tmp%mod) {
                int x = a[i + j], y = omega * a[i + j + m] % mod;
                a[i + j] = (x + y) % mod, a[i + j + m] = (x - y + mod) % mod;
            }
        }
    }
    if(inv == -1)
    {
        int nI = q_pow(n, mod-2);
        for(int i = 0; i < n; i++) {
            a[i] = a[i] * nI % mod;
        }
    }
} 

例1：http://www.acmicpc.sdnu.edu.cn/problem/show/1532

#include<bits/stdc++.h>

using namespace std;
#define int long long
#define g 3
const int mod = 998244353;
const int N = 300000;

inline int read(){
    int x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
inline void write(int x)
{
    if(x<0)x=-x,putchar('-');
    if(x>9)write(x/10);
    putchar(x%10+'0');
}

char sa[N], sb[N];
int n = 1, a[N], b[N], rev[N], lena, lenb;

inline int q_pow(int a, int b){
    int ans = 1;
    while(b > 0){
        if(b & 1){
            ans = ans * a % mod;
        }
        a = a * a % mod;
        b >>= 1; 
    } 
    return ans;
}

inline get_rev()
{
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1 | ((i & 1) << (lim - 1)));
    }
}

inline void ntt(int *a, int inv) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) swap(a[i], a[rev[i]]);    //i < t的限制使得每对点只被交换一次(否则交换两次相当于没交换) 
    }
    for(int l = 2; l <= n; l *= 2) {    //区间长度 
        int m = l / 2;
        int tmp = q_pow(g, (mod-1)/l);
        if(inv == -1) tmp = q_pow(tmp, mod-2);
        for(int i = 0; i < n; i += l) {
            int omega = 1;
            for(int j = 0; j < m; j++, omega = omega*tmp%mod) {
                int x = a[i + j], y = omega * a[i + j + m] % mod;
                a[i + j] = (x + y) % mod, a[i + j + m] = (x - y + mod) % mod;
            }
        }
    }
    if(inv == -1)
    {
        int nI = q_pow(n, mod-2);
        for(int i = 0; i < n; i++) {
            a[i] = a[i] * nI % mod;
        }
    }
} 



signed main()
{
    while(~scanf("%s%s",sa,sb))
    {
        n = 1;
        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));

        lena = strlen(sa), lenb = strlen(sb);
        for(int i = 0; i < lena; i++) {
            a[i] = sa[lena - 1 - i] - '0';
        }
        for(int i = 0; i < lenb; i++) {
            b[i] = sb[lenb - 1 - i] - '0';
        }
        while(n < lena + lenb) n <<= 1;
        get_rev();

        ntt(a, 1);
        ntt(b, 1);
        for(int i = 0; i < n; i++) {
            a[i] = a[i] * b[i] % mod;
        }
        ntt(a, -1);

        for(int i = 0; i < n; i++) {
            a[i + 1] += a[i] / 10;
            a[i] %= 10;
        }
        int cnt = n;
        while(cnt >= 0 && a[cnt] == 0) cnt--;
        if(cnt == -1) {
            printf("0");
        }
        else {
            for(int i = cnt; i >= 0; i--){
                write(a[i]);
            }
        }
        putchar('\n');
    }
    return 0;
}