FFT--luoguP3723 [AH2017/HNOI2017]礼物

本文探讨了如何通过数学技巧最小化两组数据之间的误差平方和。通过设定一个增加的自然数并应用二次函数原理,文章详细解释了求解过程,并提供了一种使用FFT的方法来寻找最优解,同时注意到了精度问题。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

传送门
可以设增加的自然数为 c c c,原式就是 ∑ i = 1 n ( x i + c − y i ) 2 \sum_{i=1}^n(x_i+c-y_i)^2 i=1n(xi+cyi)2
展开以后就是 ∑ i = 1 n ( x i 2 + 2 × x i × c − 2 × x i × y i − 2 × y i × c + y i 2 ) \sum_{i=1}^n(x_i^2+2\times x_i\times c-2\times x_i\times y_i-2\times y_i\times c+y_i^2) i=1n(xi2+2×xi×c2×xi×yi2×yi×c+yi2)
把式子中不变的提出来,就变成了
∑ i = 1 n ( x i 2 + y i 2 ) + n × c 2 + 2 × ( ∑ x i − ∑ y i ) × c − ∑ ( 2 × x i × y i ) \sum_{i=1}^n (x_i^2+y_i^2) +n\times c^2 +2\times (\sum x_i-\sum y_i)\times c-\sum(2\times x_i\times y_i) i=1n(xi2+yi2)+n×c2+2×(xiyi)×c(2×xi×yi)

这样的话前面一开始输入的时候就可以求出,中间是一个二次函数,可以用二次函数求最值得方法求出,后面是一个可变的式子,要求它的最大值,就可以通过反转 x x x倍长 y y y F F T FFT FFT来求。

注意因为 c c c是整数,但二次函数顶点坐标不一定是整数,所以要 − 1 , + 1 -1,+1 1,+1都算一下取 m i n min min

还有这题卡精度!!!最后要 + e p s +eps +eps而且 e p s = 1 e − 2 eps=1e-2 eps=1e2比较好, w a wa wa的五次都贡献给精度了
最后放上代码:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define maxn 400005
#define inf 0x7fffffff
using namespace std;
const double Pi=acos(-1.0);
const double eps=1e-2;

inline int rd(){
	int x=0,f=1;char c=' ';
	while(c<'0' || c>'9') f=c=='-'?-1:1,c=getchar();
	while(c<='9' && c>='0') x=x*10+c-'0',c=getchar();
	return x*f; 
}

struct complex{
	double x,y;
	complex(double xx=0,double yy=0) {x=xx,y=yy;}
}a[maxn],b[maxn];

complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

int n,m,rev[maxn],limit=1,l,sumx,sumy,sqsum,c,ans,mx;

void FFT(complex *F,int type){
	for(int i=0;i<limit;i++)
		if(i<rev[i]) swap(F[i],F[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		complex Wn(cos(Pi/mid),type*sin(Pi/mid));
		for(int r=mid<<1,j=0;j<limit;j+=r){
			complex w(1,0);
			for(int k=0;k<mid;k++,w=w*Wn){
				complex x=F[j+k],y=F[j+mid+k]*w;
				F[j+k]=x+y,F[j+mid+k]=x-y;
			}
		}
	}
}

int main(){
	n=rd(); m=rd();
	for(int i=0;i<n;i++) scanf("%lf",&a[i].x),sumx+=a[i].x,sqsum+=a[i].x*a[i].x;
	for(int i=0;i<n;i++) scanf("%lf",&b[i].x),sumy+=b[i].x,sqsum+=b[i].x*b[i].x;
	c=(sumy-sumx)/n; ans=inf;
	for(int i=-1;i<=1;i++) {
		int tmp=c+i;
		ans=min(ans,n*tmp*tmp+2*(sumx-sumy)*tmp);
	}
	for(int i=0;2*i<n;i++) swap(a[i],a[n-i-1]);
	for(int i=0;i<n;i++) b[i+n].x=b[i].x;
	while(limit<3*n) limit<<=1,++l;
	for(int i=0;i<limit;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	FFT(a,1); FFT(b,1);
	for(int i=0;i<limit;i++) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=0;i<limit;i++) mx=max(mx,(int)((a[i].x)/(double)limit+eps));
	printf("%d\n",sqsum+ans-2*mx);
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值