发布于

lp3803 【模板】多项式乘法（FFT）（非递归版本）

快速傅里叶变换（FFT）的一种较高效的实现

*由于快速傅里叶变换（FFT）常常被用在通信领域，故而FFT的实现效率成为了一个需要被考察的对象。

事实上，观察FFT的实现过程，我们发现它没有非递归不可的理由。

具体来说，我们发现，如果我们想要把FFT改造为递推模式，最大的障碍就是每一次合并操作。
在快速离散傅里叶变换中，每一次合并操作都会导致两个距离为\(\frac{n}{2}\)的点的值合并在一起。因此，如果将原序列传入，那么得到的就是一个乱序的队列。
我们考虑将序列以什么顺序传入。
我们发现，每一次这样的合并操作本质上都是将两个部分「翻转」了。
以0~7为例：
0 1 2 3 4 5 6 7
0 4 2 6 1 5 3 7
将它们写成二进制：
000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111
发现变换之后每个数的二进制位都倒转了。
故而我们可以通过倒转二进制位的方式来预处理原序列，使得FFT之后的序列顺序是正常的。
但是不知道是什么原因，这个代码跑出来的速度竟然也是4500ms…这就很让人挠头了。
于是向PinkRabbit巨神请教了几个经验：
1.不用中间数组，直接在原数组上操作。
2.调整赋值操作，使得修改尽量少。
3.将蝴蝶变换预处理为一个数组。
综合几点，就成功地卡到了2500ms
当然还是有待卡常，毕竟这还是一个比较危险的时间。

#include<iostream>
#include<cstdio>
#include<cmath>
#define PI 3.141592653589793238462

struct Complex{
	double r;
	double i;
	Complex(double ri=0.0,double ii=0.0):r(ri),i(ii){}
	inline Complex operator+(const Complex &B)const{
		return (Complex){r+B.r,i+B.i}; 
	}
	inline Complex operator-(const Complex &B)const{
		return (Complex){r-B.r,i-B.i};
	}
	inline Complex operator*(const Complex &B)const{
		return (Complex){r*B.r-i*B.i,r*B.i+i*B.r};
	}
}a[1<<21|1],b[1<<21|1],s[1<<21|1];
int r[1<<21|1];
inline void FFT(Complex *A,int N,int typ){
	//ButterflyTransform
	for(int i=0;i<N;++i){
		(r[i]>i)?std::swap(A[i],A[r[i]]),0:0;
	}
	Complex bs,nw,x,y;
	for(int i=2,mid;i<=N;i<<=1){
		bs=Complex(cos(2.0*PI/i),typ*sin(2.0*PI/i));
		mid=i>>1;
		for(int j=0;j<N;j+=i){
			nw=Complex(1.0,0.0);
			for(int k=0;k<mid;++k,nw=nw*bs){
				x=A[j+k],y=(A[mid+j+k]*nw);
				A[j+k]=x+y;
				A[mid+j+k]=x-y;
			}
		}
	}
}
int n,m,cnt=1,byt=0;;
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%lf",&a[i].r);
	}
	for(int i=0;i<=m;++i){
		scanf("%lf",&b[i].r);
	}
	while(cnt<=n+m){
		cnt<<=1;
		++byt;
	}
	for(int i=0;i<cnt;++i){
		r[i]=(r[i>>1]>>1)|((i&1)<<(byt-1));
	}
	FFT(a,cnt,1);
	FFT(b,cnt,1);
	for(int i=0;i<cnt;++i){
		a[i]=a[i]*b[i];
	}
	FFT(a,cnt,-1);
	//注意这里typ应该是-1 
	for(int i=0;i<=n+m;++i){
		printf("%d ",(int)(a[i].r/cnt+0.5));
		//注意输出应该强转格式。 
	}
	puts("");
}

int main(){
	init();
	return 0; 
}
发布于

lp3195 HNOI2008 玩具装箱TOY

一道斜率优化的经典例题。
对于这道题，DP式是容易得到的：
$$f_i=min\{f_j+(\sum_{j+1}^ia_{i}+i-j-1-L)^2\}$$
前缀和优化：
$$f_i=min\{f_j+(sm_i-sm_j+i-j-1-L)^2\}$$
考虑到这是一个与取值平方相关的式子，尝试用斜率优化来化简。
令：
$$u_i=sm_i+i$$
$$f_i=min{f_j+(u_i-u_j-1-L)^2}$$
对于右式，我们假定\(1\le j<k<i\)且从\(k\)转移较优。
那么我们有：
$$f_{k}+(u_{i}-u_{k}-1-L)^2\le f_{j}+(u_{i}-u_{j}-1-L)^2$$
展开并移项化简：
$$2u_{i}(u_{j}-u_{k})\le f_{j}-f_{k}+(u_{j}+1+L)^2-(u_{k}+1+L)^2$$
令：
$$g_{i}=(u_{i}+1+L)^2$$
并不妨假定：
$$u_{j}!=u_{k}$$
可得：
$$2u_{i}\le\frac{f_{j}-f_{k}+(u_{j}+1+L)^2-(u_{k}+1+L)^2}{(u_{j}-u_{k})}$$
故而，若要使得答案更大，则上式必然成立。
考虑到这个式子符合单调队列的性质，我们可以上一个类似单调队列的数据结构来维护右侧的式子，也就是俗称的“斜率”。
这就完成了斜率优化。

#include<iostream>
#include<cstdio>
#include<queue>


int n,L;
double sm[50005],f[50005];
int q[50005],l,r;
inline double u(int X){
	return sm[X]+X;
}
inline double g(int X){
	return u(X)+L+1;
}
inline double x(int X){
	return g(X);
}
inline double y(int X){
	return f[X]+g(X)*g(X);
}
inline double k(int A,int B){
	return (y(A)-y(B))/(x(A)-x(B));
}

void init(){
	scanf("%d%d",&n,&L);
	for(int i=1;i<=n;++i){
		scanf("%lf",&sm[i]);
		sm[i]+=sm[i-1];
	}
	l=r=1;
	for(int i=1;i<=n;++i){
		while(l<r&&k(q[l],q[l+1])<2*u(i)){
			++l; 
		} 
		f[i]=f[q[l]]+(u(i)-g(q[l]))*((u(i)-g(q[l])));
		while(l<r&&k(i,q[r-1])<k(q[r-1],q[r])){
			--r;
		}
		q[++r]=i;
	}
	printf("%lld\n",(long long)f[n]);
}

int main(){
	init();
	return 0;
}