FFT NNT

算算劳资已经多久没学新算法了，又要重新开始学辣。直接扔板子，跑。。。话说FFT算法导论里讲的真不错，去看下就懂了。

//FFT
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
typedef __int64 ll;
const double pi = acos(-1.0);
const int N = 10000+5;
const double eps = 1e-6;

struct Complex {
        double x,y;
        Complex(double _x=0,double _y=0)
                        :x(_x),y(_y) {}
        Complex operator + (const Complex& tt)const { return Complex(x+tt.x,y+tt.y); }
        Complex operator - (const Complex& tt)const { return Complex(x-tt.x,y-tt.y); }
        Complex operator * (const Complex& tt)const { return Complex(x*tt.x-y*tt.y,x*tt.y+y*tt.x); }
};
void FFT(Complex *a, int n, int rev) {
        // n是大于等于相乘的两个数组长度的2的幂次
        // 从0开始表示长度，对a进行操作
        // rev==1进行DFT，==-1进行IDFT
        for (int i = 1,j = 0; i < n; ++ i) {
                for (int k = n>>1; k > (j^=k); k >>= 1);
                if (i<j) swap(a[i],a[j]);
        }
        for (int m = 2; m <= n; m <<= 1) {
                Complex wm(cos(2*pi*rev/m),sin(2*pi*rev/m));
                for (int i = 0; i < n; i += m) {
                        Complex w(1.0,0.0);
                        for (int j = i; j < i+m/2; ++ j) {
                                Complex t = w*a[j+m/2];
                                a[j+m/2] = a[j] - t;
                                a[j] = a[j] + t;
                                w = w * wm;
                        }
                }
        }
        if (rev==-1) {
                for (int i = 0; i < n; ++ i) a[i].x = (a[i].x/n+eps);
        }
}

//NNT
#include <cstdio>
#include <cmath>
#include <vector>
#include <cstring>
#include <algorithm>
using namespace std;
typedef __int64 ll;
const double pi = acos(-1.0);
const int N = 100000+5;
const double eps = 1e-6;
const int MOD = 998244353; 

int f[N<<1];

int pow_mod(int x, int n) {
    int ret = 1;
    while(n) {
        if(n&1) ret = (ll)ret*x%MOD;
        x = (ll)x*x%MOD;
        n >>= 1;
    }
    return ret;
}

int Top, m1[N<<3], m2[N<<3];

void FFT(int *a, int n, int rev) {
    for(int i = 1, j = 0;i < n; i++) {
        for(int k = n>>1;k > (j^=k); k >>= 1) ;
        if(i < j)   swap(a[i], a[j]);
    }
    for(int i = 0;i < n; i++)   m2[i] = m1[i*(Top/n)];
    for(int m = 2;m <= n; m <<= 1) {
        int wm = rev == 1 ? n-n/m : n/m;
        for(int i = 0;i < n; i += m) {
            int w = 0;
            for(int j = i;j < i+m/2; j++) {
                int t = (ll)m2[w]*a[j+m/2]%MOD;
                a[j+m/2] = a[j] - t;
                if(a[j+m/2] < 0)    a[j+m/2] += MOD;
                a[j] = a[j] + t;
                if(a[j] >= MOD) a[j] -= MOD;
                w += wm;
                if(w >= n)  w -= n;
            }
        }
    }
    if(rev == -1) {
        int inv = pow_mod(n, MOD-2);
        for(int i = 0;i < n; i++) a[i] = (ll)a[i]*inv%MOD;
    }
}

 

转载于:https://www.cnblogs.com/get-an-AC-everyday/p/5702482.html