【trick】myyFFT

71 阅读 0 评论 47 点赞

我是靠谱客的博主糟糕花瓣，最近开发中收集的这篇文章主要介绍【trick】myyFFT，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

说在前面

发现 $C l a r i s$ 的代码里 $F F T$ 卷积的姿势不太一样，于是学习了一发 $m y y F F T$ 。

三次 $F F T$ 可以通过 $m y y F F T$ 优化成两次，卡常数必备 $t r i c k$ 。orzmyy

共轭复数

两个实部相等，虚部互为相反数的复数互为共轭复数( $c o n j u g a t e c o m p l e x n u m b e r$ )。将复数 $z$ 的共轭复数记做 $c o n j (z)$ 。

则有运算:

$x = c o n j (c o n j (x))$
$c o n j (x) c o n j (y) = c o n j (x y)$
$c o n j (x + y) = c o n j (x) + c o n j (y)$

myyFFT

设当前卷积的两个多项式分别为 $A, B$ ，设进行 $F F T$ 的多项式为 $C$ 。
复数 $z$ 的实部记做 $z . r$ ，虚部记做 $z . i$ 。

则 $C_k=A_k+B_k·i$ 。另设多项式 $D,D_k=A_k-B_k·i$ 。

由 $conj(omega_n^k)=omega_n^{-k}$ ，得：

$C(omega_n^k)=sum limits_{j=0}^{n-1}(A_j+B_j·i)omega_n^{jk}$

$D(omega_n^k)=sum limits_{j=0}^{n-1}(A_j-B_j·i)omega_n^{jk}=sumlimits_{j=0}^{n-1}conj((A_j+B_j·i)conj(omega_{n}^{jk}))=conj(sum_{j=0}^{n-1}(A_j+B_j·i)omega_{n}^{-jk})=conj(C(omega_n^{-k}))=conj(C(omega_n^{n-k}))$

在将 $C$ 转化成点积表示后，由：
$A = C + D 2 , B = C − D 2 i A=dfrac{C+D}{2},B=dfrac{C-D}{2i}$
得到 $A * B$ 的卷积表示为： $dfrac{C^2-D^2}{4i}=dfrac{C^2-D^2}{-4dfrac{1}{i}}=(D^2-C^2)dfrac 14i$

代码

uoj34:多项式乘法

#include<bits/stdc++.h>
#define RI register
#define db double
using namespace std;
const int N=5e5+10;
const db pi=acos(-1);

int n,m,L,len,rv[N];

struct cc{
	db r,i;
	cc(){r=i=0;};
	cc(db r_,db i_){r=r_;i=i_;}
	cc operator +(const cc&ky){return cc(r+ky.r,i+ky.i);}
	cc operator -(const cc&ky){return cc(r-ky.r,i-ky.i);}
	cc operator *(const cc&ky){return cc(r*ky.r-i*ky.i,r*ky.i+i*ky.r);}
	cc operator /(const int&ky){return cc(r/(db)ky,i/(db)ky);}
	inline cc conj(){return cc(r,-i);}
}a[N],b[N];

inline void FFT(cc *e,int ptr)
{
	RI int i,j,k,t;cc ix,iy,bs,ori;
	for(i=1;i<len;++i) if(i<rv[i]) swap(e[i],e[rv[i]]);
    for(i=1;i<len;i<<=1){
    	ori=cc(cos(pi/i),ptr*sin(pi/i));
    	for(j=0;j<len;j+=(i<<1)){
    		bs=cc(1,0);
    		for(k=0;k<i;++k,bs=bs*ori){
    			ix=e[j+k];iy=bs*e[i+j+k];
    			e[j+k]=ix+iy;e[i+j+k]=ix-iy;
    		}
    	}
    }
    if(ptr==1) return;
    for(i=0;i<len;++i) e[i]=e[i]/len;
}

int main(){
	RI int i,j;
	scanf("%d%d",&n,&m);
	for(i=0;i<=n;++i) scanf("%lf",&a[i].r);
	for(i=0;i<=m;++i) scanf("%lf",&a[i].i);
	for(len=1,L=0;len<=n+m;len<<=1) L++;
	for(i=1;i<len;++i) rv[i]=(rv[i>>1]>>1)|((i&1)<<(L-1)); 
	FFT(a,1);
	for(i=0;i<len;++i){
		j=(len-i)&(len-1);
		b[i]=(a[i]*a[i]-(a[j]*a[j]).conj())*cc(0,-0.25);
	}
	FFT(b,-1);
	for(i=0;i<=n+m;++i) 
	 printf("%d ",int(b[i].r+0.5));
	return 0;
}