uvalive6886 快速傅里叶变换模板

来源:互联网 发布:政务数据开放流程 编辑:程序博客网 时间:2024/06/01 08:23
#include <bits/stdc++.h>using namespace std;const int INF  = 0x3f3f3f3f;const int maxn = 200060;const int Mod  = 1e9 + 7;#define ll       long long#define mem(x,y) memset(x,y,sizeof(x))#define IO       ios_base::sync_with_stdio(0), cin.tie(0);inline ll gcd(ll a, ll b) {return a % b == 0 ? b : gcd(b, a % b);}inline ll lcm(ll a, ll b) {return a / gcd(a, b) * b;}inline ll quick_pow(ll x, int k) {ll ans = 1; while (k) { if (k & 1) ans = (ans * x) % Mod; x = x * x % Mod;  k >>= 1; } return ans;}const double eps(1e-8);typedef long long lint;const double PI = acos(-1.0);struct Complex{    double real, image;    Complex(double _real, double _image)    {        real = _real;        image = _image;    }    Complex() {    }};Complex operator + (const Complex &c1, const Complex &c2){    return Complex(c1.real + c2.real, c1.image + c2.image);}Complex operator - (const Complex &c1, const Complex &c2){    return Complex(c1.real - c2.real, c1.image - c2.image);}Complex operator * (const Complex &c1, const Complex &c2){    return Complex(c1.real * c2.real - c1.image * c2.image, c1.real * c2.image + c1.image * c2.real);}int rev(int id, int len){    int ret = 0;    for (int i = 0; (1 << i) < len; i++)    {        ret <<= 1;        if (id & (1 << i)) ret |= 1;    }    return ret;}Complex A[maxn * 3] ;void FFT(Complex* a, int len, int DFT)//对a进行DFT或者逆DFT, 结果存在a当中{    //Complex* A = new Complex[len]; 这么写会爆栈    for (int i = 0; i < len; i++)        A[rev(i, len)] = a[i];    for (int s = 1; (1 << s) <= len; s++)    {        int m = (1 << s);        Complex wm = Complex(cos(DFT * 2 * PI / m), sin(DFT * 2 * PI / m));        for (int k = 0; k < len; k += m)        {            Complex w = Complex(1, 0);            for (int j = 0; j < (m >> 1); j++)            {                Complex t = w * A[k + j + (m >> 1)];                Complex u = A[k + j];                A[k + j] = u + t;                A[k + j + (m >> 1)] = u - t;                w = w * wm;            }        }    }    if (DFT == -1) for (int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len;    for (int i = 0; i < len; i++) a[i] = A[i];    return;}int numA[maxn], numB[maxn];//以每一位为系数, 那么多项式长度不超过Complex a[maxn * 3], b[maxn * 3]; //对应的乘积的长度不会超过int ans[maxn * 3];//  多项式 t =  a0*x^0 + a1*x^1 + .... + an*x^n//  numA[i] 保存多项式A的x^i 的系数, numB[i] 保存多项式B的x^i 的系数,C=A*B//  ans[i] 保存多项式C的x^i 的系数,int main(){    int N, M;    while (scanf("%d", &N) != EOF)    {        int tmp, L = 0, answer = 0;        mem(numA, 0);        mem(numB, 0);        mem(ans, 0);        for (int i = 1; i <= N; i++) {            scanf("%d", &tmp);            numA[tmp] = 1, numB[tmp] = 1;            L = max(L, tmp);        }        int lenA = L + 1;        int sa = 0;        while ((1 << sa) < lenA) sa++;        int lenB = L + 1;        int sb = 0;        while ((1 << sb) < lenB) sb++;        //那么乘积多项式的次数不会超过(1 << (max(sa, sb) + 1))        int len = (1 << (max(sa, sb) + 1));        //cout << len << endl;        for (int i = 0; i < len; i++)        {            if (i < lenA) a[i] = Complex(numA[i], 0);            else a[i] = Complex(0, 0);            if (i < lenB) b[i] = Complex(numB[i], 0);            else b[i] = Complex(0, 0);        }        FFT(a, len, 1);        FFT(b, len, 1);//把A和B换成点值表达        for (int i = 0; i < len; i++) //做点值表达的成乘法            a[i] = a[i] * b[i];        FFT(a, len, -1);//逆DFT换回原来的系数, 虚部一定是0        for (int i = 0; i < len; i++)            ans[i] = (int)(a[i].real + 0.5); //取整误差的处理        cin >> M;        for (int i = 1; i <= M; i++) {            scanf("%d", &tmp);            if (ans[tmp] || numA[tmp])                answer++;        }        cout << answer << endl;    }    return 0;}
原创粉丝点击