自动搬运

查看原文

来自洛谷，原作者为

	littlez_meow 这个小z很懒，连这个家伙很懒都没有留下

搬运于2025-08-24 22:49:26，当前版本为作者最后更新于2024-05-30 20:05:30，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

建议改为，偏导数在 OI 中的运用。

题目指路。

题意

给定 $n,m,k,E$ 和长度分别为 $n,m$ 的序列 $p_1\sim p_n,q_1\sim q_m$，记：

$$f(x)=\sum\limits_{i=1}^k\left(a_i\cos(ix)+b_i\sin(ix)\right) $$$$\bar{u}=\frac{1}{n} \sum\limits_{i=1}^{n} f\left(p_{i}\right) $$$$\bar{v}=\frac{1}{m} \sum\limits_{j=1}^{m} f\left(q_{i}\right) $$$$\operatorname{cost}(f)=\frac{1}{\bar{u}-\bar{v}} \sqrt{\sum\limits_{i=1}^{n}\left(f\left(p_{i}\right)-\bar{u}\right)^{2}+\sum\limits_{j=1}^{m}\left(f\left(q_{j}\right)-\bar{v}\right)^{2}} $$

你需要给出两个长为 $k$ 的序列 $a_1\sim a_k,b_1\sim b_k$，使得 $\bar{u}>\bar{v},\operatorname{cost}(f)$ 最小。绝对误差或相对误差小于 $10^{-E}$。

思路

下面简记 $\operatorname{cost}(f)$ 为 $c(f)$。

第一眼，$c$ 是一个函数到实数的映射。泛函？这怎么做？

再想一想，其实 $c$ 的值只和 $a_1\sim a_k,b_1\sim b_k$ 这 $2k$ 个取值有关。也就是说，$c$ 是一个 $2k$ 元函数，记为 $c(f)=c(a_1,b_1,\cdots,a_k,b_k)$。很好，至少是个正常的函数了。

$\bar{u}>\bar{v}$ 的限制太烦了，有没有什么办法能绕开呢？

有！如果 $\bar{u}<\bar{v}$，我们直接把所有 $a_i,b_i$ 全部乘 $-1$。这样，函数值乘 $-1$，方差不变，均值之差乘 $-1$，此时的 $\bar{u}>\bar{v}$。

也就是说，我们可以直接计算 $c^2(f)$。这样，根号去掉了，$\bar{u}<\bar{v}$ 时的负数取值也去掉了。根据上面的论述，让 $c(f)$ 取到最小值的 $f$ 也可以让 $c^2(f)$ 取到最小值。这容易用反证法证明。

设 $C(f)=C(a_1,b_1,\cdots,a_k,b_k)=c^2(f)=\dfrac{\sum\limits_{i=1}^n(f(p_i)-\bar{u})^2+\sum\limits_{i=1}^m(f(q_i)-\bar{v})^2}{(\bar{u}-\bar{v})^2}$，我们只需要研究 $C(f)$ 的最值就好了，且没有 $\bar{u}>\bar{v}$ 的要求。

下一步怎么办呢？直接展开 $\bar{u},\bar{v},f$？这样得到的结果会不会太复杂了呢？

最小二乘法的经验告诉我们，我们可以尝试用向量和矩阵来描述上述式子，往往能大幅简化式子。

设 $2k$ 维向量 $\vec{t}=[a_1,b_1,\cdots,a_k,b_k]^T,\vec{f}(x)=[\cos(x),\sin(x),\cdots,\cos(kx),\sin(kx)]^T$。不难验证 $f(x)=\vec t\cdot\vec f(x)$。下面如果没有特殊说明，向量都是 $2k$ 维的。

现在 $C(f)$ 就是一个关于 $\vec t$ 的函数了，记为 $C(\vec t)$。

设 $\vec p_i=\vec f(p_i),\vec q_i=\vec f(q_i)$。我们可以直接算出所有 $\vec p_i,\vec q_i$，$C(f)$ 就跟 $\cos,\sin$ 没有半点关系了。

现在来看看 $\bar{u}$ 是怎么算出来的吧。

写出原式：$\bar u=\dfrac 1 n\sum\limits_{i=1}^nf(p_i)$。

带入 $f(x)=\vec t\cdot\vec f(x)$，得到：$\bar u=\dfrac 1 n\sum\limits_{i=1}^n\vec t\cdot\vec p_i$。

使用向量点乘和加法的分配律：$\bar u=\vec t\cdot(\dfrac 1 n\sum\limits_{i=1}^n\vec p_i)$

右边就是一个常向量了，记 $\vec r=\dfrac 1 n\sum\limits_{i=1}^n\vec p_i$，则 $\bar u=\vec t\cdot\vec r$。

同理，记 $\vec s=\dfrac 1 m\sum\limits_{i=1}^m\vec q_i$，则 $\bar v=\vec t\cdot\vec s$。

现在，把这些向量带回原式：

$$C(\vec t)=\dfrac{\sum\limits_{i=1}^n(\vec t\cdot(\vec p_i-\vec r))((\vec p_i-\vec r)\cdot\vec t)+\sum\limits_{i=1}^m(\vec t\cdot(\vec q_i-\vec s))((\vec q_i-\vec s)\cdot\vec t)}{(\vec t\cdot(\vec r-\vec s))^2} $$

点积破坏了我们很多美好的性质，可不可以去掉它呢？

对于列向量 $\vec a,\vec b$，容易发现 $\vec a^T\vec b=[\vec a\cdot\vec b]$。我们只要认为 $1\times1$ 的矩阵就是一个数，就可以用矩阵乘法代替点积。

用矩阵乘法代替点积，使用结合律和分配律，得到：

$$C(\vec t)=\dfrac{\vec t^T(\sum\limits_{i=1}^n(\vec p_i-\vec r)(\vec p_i-\vec r)^T+\sum\limits_{i=1}^m(\vec q_i-\vec s)(\vec q_i-\vec s)^T)\vec t}{\vec t^T((\vec r-\vec s)(\vec r-\vec s)^T)\vec t} $$

好！中间括号里面的都变成常矩阵了！设矩阵 $A=\sum\limits_{i=1}^n(\vec p_i-\vec r)(\vec p_i-\vec r)^T+\sum\limits_{i=1}^m(\vec q_i-\vec s)(\vec q_i-\vec s)^T$，矩阵 $B=(\vec r-\vec s)(\vec r-\vec s)^T$，我们就可以化简为：

$$C(\vec t)=\dfrac{\vec t^TA\vec t}{\vec t^TB\vec t} $$

非常美的一个结果。

现在就是多元函数求最值问题了，方法是对每一元求偏导并使其等于 $0$。

设 $A_{i,j}$ 表示 $A$ 第 $i$ 行第 $j$ 列的元素，$B_{i,j}$ 同理。

为了方便，设 $x_{2i+1}=a_i,x_{2i}=b_i$。

展开矩阵并化简，我们有:

$$C(f)=\dfrac{\sum\limits_{i=1}^{2k}\sum\limits_{j=1}^{2k}x_ix_jA_{i,j}}{\sum\limits_{i=1}^{2k}\sum\limits_{j=1}^{2k}x_ix_jB_{i,j}} $$

我们发现分子分母都是关于某个自变量的二次函数，容易求偏导。

求导后，我们的结果中还存在一些高次项。加加减减消掉它们，我们发现系数又变成了矩阵的元素！只是等式右边带一个非零常数。即有以下方程：

$$A\vec t=(\vec r-\vec s)t$$

这也是合理的。因为 $a,b$ 等比例放大，标准差也等比例放大，平均数同样等比例放大，在分式上的比例最后就消掉了。样例解释一也说明了这点。高斯消元解之即可。

最后，还有高斯消元无解的情况。注意到数据范围里说 $p,q$ 随机生成，直接赌不会无解就好。

时间复杂度 $O(k^3)$，$k,n,m$ 视为同阶。

如果你的 $p,q$ 是整型，记得转成浮点型，否则你会因为溢出只有二十四分。

代码

#include<bits/stdc++.h>
#define F(i,a,b) for(int i(a),i##i##end(b);i<=i##i##end;++i)
#define R(i,a,b) for(int i(a),i##i##end(b);i>=i##i##end;--i)
#define File(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
using namespace std;
const int MAXN=601;
int n,m,k,kk,e,p[MAXN],q[MAXN];
double vecp[MAXN][MAXN],vecq[MAXN][MAXN],r[MAXN],s[MAXN];
double factor[MAXN][MAXN+1];
double a[MAXN],b[MAXN];
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0),cout.tie(0);
	cin>>n>>m>>k>>e;
	kk=k<<1;
	F(i,1,n){
		cin>>p[i];
		F(j,1,k) r[(j<<1)-1]+=(vecp[i][(j<<1)-1]=cos(j*1.0*p[i])),r[j<<1]+=(vecp[i][j<<1]=sin(j*1.0*p[i]));
	}
	F(i,1,m){
		cin>>q[i];
		F(j,1,k) s[(j<<1)-1]+=(vecq[i][(j<<1)-1]=cos(j*1.0*q[i])),s[j<<1]+=(vecq[i][j<<1]=sin(j*1.0*q[i]));
	}
	F(i,1,kk) r[i]/=n,s[i]/=m,factor[i][kk+1]=r[i]-s[i];
	F(i,1,kk) F(j,1,kk){
		F(t,1,n) factor[i][j]+=(vecp[t][i]-r[i])*(vecp[t][j]-r[j]);
		F(t,1,m) factor[i][j]+=(vecq[t][i]-s[i])*(vecq[t][j]-s[j]);
	}
	F(i,1,kk){
		int mx=i;
		F(j,i+1,kk) fabs(factor[j][i])>fabs(factor[mx][i])&&(mx=j);
		swap(factor[i],factor[mx]);
		F(j,i+1,kk){
			double delta=factor[j][i]/factor[i][i];
			F(t,i,kk+1) factor[j][t]-=factor[i][t]*delta;
		}
	}
	R(i,kk,1){
		if(i&1){
			double&qwq=a[(i+1)>>1];
			qwq=factor[i][kk+1]/factor[i][i];
			R(j,i-1,1) factor[j][kk+1]-=qwq*factor[j][i];
		}
		else{
			double&qwq=b[(i+1)>>1];
			qwq=factor[i][kk+1]/factor[i][i];
			R(j,i-1,1) factor[j][kk+1]-=qwq*factor[j][i];
		}
	}
	double u=0,v=0;
	F(i,1,k){
		u+=r[(i<<1)-1]*a[i]+r[i<<1]*b[i];
		v+=s[(i<<1)-1]*a[i]+s[i<<1]*b[i];
	}
	if(u<v) F(i,1,k) a[i]=-a[i],b[i]=-b[i];
	cout<<fixed<<setprecision(12);
	F(i,1,k) cout<<a[i]<<" "<<b[i]<<"\n";
	return 0;
}