数学板块学习之拆系数FFT

拆系数FFT

比任意模数NTT要快一点的的算法(要看情况，洛谷的一个板子题我NTT跑了1071ms,而FFT跑了4927ms)
取一个$M=1<<15$
将$F[i],G[i]$变为
$$\begin{array}{rl}F[i]& =A[i]\ast M+B[i]\\ G[i]& =C[i]\ast M+D[i]\end{array}$$

所以$F[i]\ast G[i]$就变为了
$$\begin{array}{rl}F[i]\ast G[i]& =(A[i]\ast M+B[i])\ast (C[i]\ast M+D[i])\\ & =A[i]\ast C[i]\ast M^2+(B[i]\ast C[i]+A[i]\ast D[i])\ast M+B[i]\ast D[i]\end{array}$$

所以我们要求的是$A[i]\ast C[i],B[i]\ast C[i],A[i]\ast D[i],B[i]\ast D[i]$
所以最暴力的也就是八次DFT

模板代码(八次DFT):

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <cmath>
#include <string>
#include <cstring>
#include <map>
#include <set>
#include <math.h>
#include <ctime>
#include <unordered_map>
//#include <tr1/unordered_map>

using namespace std;
#define Time_begin time_t begin,end;begin=clock()
#define Time_end end=clock();cout<<"runtime:   "<<double(end-begin)/CLOCKS_PER_SEC*1000<<"ms"<<endl
#define me(x,y) memset(x,y,sizeof x)
#define MIN(x,y) x < y ? x : y
#define MAX(x,y) x > y ? x : y

typedef long long ll;
typedef unsigned long long ull;

const long double INF = 0x3f3f3f3f;
const int MOD = 1e9+7;
const double eps = 1e-06;
const long double PI = std::acos(-1);
const int M=32768;
const int MAXN = 5e5+100;


struct Complex{
	long double x,y;//实部和虚部 x+yi
	Complex(long double _x = 0.0,long double _y = 0.0){
		x = _x;
		y = _y;
	}
	Complex operator - (const Complex &b)const{
		return Complex(x - b.x,y - b.y);
	}	
	Complex operator +(const Complex &b)const{
		return Complex(x+b.x,y+b.y);
	}
	Complex operator *(const Complex &b)const{
		return Complex(x*b.x - y*b.y,x*b.y+y*b.x);
	}
};
/*
* 进行 FFT 和 IFFT 前的反转变换。
* 位置 i 和(i 二进制反转后位置)互换
* len 必须为 2 的幂
*/
void change(Complex y[],int len){
	int i,j,k;
	for(i = 1, j = len/2;i <len - 1;i++){
		if(i < j)swap(y[i],y[j]);
		//交换互为小标反转的元素,i<j 保证交换一次
		//i 做正常的 +1,j 左反转类型的 +1, 始终保持 i 和 j 是反转的
		k = len/2;
		while(j >= k){
			j -= k;
			k /= 2;
		}
		if(j < k)j += k;
	}
}
/*
* 做 FFT
* len 必须为2^k形式
* on==1 时是 DFT,on==-1 时是 IDFT
*/
void fft(Complex y[],int len,int on){
	change(y,len);
	for(int h = 2; h <= len; h <<= 1){
		Complex wn(std::cos(-on*2*PI/h),std::sin(-on*2*PI/h));
		for(int j = 0;j < len;j+=h){
			Complex w(1,0);
			for(int k = j;k < j+h/2;k++){
				Complex u = y[k];
				Complex t = w*y[k+h/2];
				y[k] = u+t;
				y[k+h/2] = u - t;
				w = w*wn;
			}
		}
	}
	if(on == -1)
		for(int i = 0;i < len;i++)
			y[i].x /= (double)len;
}
ll ans[MAXN];
int n,m,p;
int F[MAXN],G[MAXN],l=0;
Complex a[MAXN],b[MAXN],c[MAXN],d[MAXN],e[MAXN],f[MAXN],g[MAXN],h[MAXN];

/*
* F[i]=a[i]*M+b[i]
* G[i]=c[i]*M+d[i] 
*/
void mtt(){
    int len=1;
    n += m;
    while(len <= n) len<<=1;
    for(int i = 0; i < n; ++i){
        a[i].x=F[i]/M,b[i].x=F[i]%M;
        c[i].x=G[i]/M,d[i].x=G[i]%M;
    }
    fft(a,len,1),fft(b,len,1),fft(c,len,1),fft(d,len,1);
    for(int i = 0 ; i < len; ++i){
        e[i]=a[i]*c[i],f[i]=a[i]*d[i];
        g[i]=b[i]*c[i],h[i]=b[i]*d[i];
    }
    fft(e,len,-1),fft(f,len,-1),fft(g,len,-1),fft(h,len,-1);
    for(int i = 0; i < len; ++i){
        ans[i]=0;
        ans[i]=(ans[i] + (ll)(round(e[i].x))%p*M%p*M%p)%p;
        ans[i]=(ans[i] + (ll)(round(f[i].x))%p*M%p)%p;
        ans[i]=(ans[i] + (ll)(round(g[i].x))%p*M%p)%p;
        ans[i]=(ans[i] + (ll)(round(h[i].x))%p)%p;
    }

}

int main() {
    // Time_begin;
    cout<<(1<<15)<<endl;
    scanf("%d%d%d",&n,&m,&p);
    for(int i = 0; i <= n; ++i) scanf("%d",&F[i]);
    for(int i = 0; i <= m; ++i) scanf("%d",&G[i]);
    mtt();
    for(int i = 0; i <= n; ++i){
        printf("%lld%c",ans[i],i == n ? '\n' : ' ');
    }
    // Time_end;
    return 0;
}

/*

*/