[NOI2014模拟][数学推导/二维欧几里得]数学题

题面

给出两个向量$\vec{a}=(x_0,y_0)$, $\vec{b}=(x_1,y_1)$，求整数$x$, $y$满足不同时为$0$，使$|x\vec{a}+y\vec{b}|$最小。
$T$组数据,每组数据输出$|x\vec{a}+y\vec{b}{|}^2$。

对于20%的数据点:$T\le 100$,存在一组最优的$(x,y)$,满足$|x|,|y|\le 30$。
对于60%的数据点:$T\le 10000$,存在一组最优的$(x,y)$,满足 min ⁡ { ∣ x ∣ , ∣ y ∣ } ≤ 30 \min \{|x|,|y|\}≤30 min{∣x∣,∣y∣}≤30。
对于100%的数据点:$T\le 100000$,$|x_0|,|y_0|,|x_1|,|y_1|\le 10^4$,存在一组最优的$(x,y)$,满足$|x|,|y|\le 10^4$。

解法1

设$(k_1,k_2)$为两整数。
我们的任务就是求$(k_1{x}_1+k_2{x}_2{)}^2+(k_1{y}_1+k_2{y}_2{)}^2$的最小值。
展开得:
$$(x_1^2+y_1^2)k_1^2+2(x_1{x}_2+y_1{y}_2)k_1{k}_2+(x_2^2+y_2^2)k_2^2$$
我们可以将式子看成一个以$k_1$为自变量的函数。
令$A=x_1^2+y_1^2,B=x_1{x}_2+y_1{y}_2,C=x_2^2+y_2^2$
原式化为$$A\times k_1^2+2B\times k_2\times k_1+C\times k_2^2$$
根据二次函数顶点式得$k_1=-\frac{B}{A}{k}_2$时函数取到最小值。
那么我们就能够得到一个60分的算法了,每次枚举$k_1$或$k_2$解出与顶点最接近的另一个数即可。
考虑继续化简，由于实际的$k_1$与取到最小的$k_1$最多相差$\frac{1}{2}$，我们可以尝试将$k_1$代入先消掉。
得:$$(C-\frac{B^2}{A})k_2^2$$
代入A,B,C,得:
$$\frac{(x_1^2+y_1^2)(x_2^2+y_2^2)-(x_1{x}_2+y_1{y}_2{)}^2}{x_1^2+y_1^2}{k}_2^2$$
拆项:
$$\frac{x_1^2{y}_2^2+x_2^2{y}_1^2-2x_1{x}_2{y}_1{y}_2}{x_1^2+y_1^2}{k}_2^2$$
合并:
$$\frac{(x_1{y}_2-x_2{y}_1{)}^2}{x_1^2+y_1^2}{k}_2^2$$
可以发现这一定是一个非负数，而且仅仅只与$k_2$有关。
这样我们就可以将原式写成一个顶点式了:
$$A(k_1+\frac{B}{A}{k}_2{)}^2+\frac{(x_1{y}_2-x_2{y}_1{)}^2}{x_1^2+y_1^2}{k}_2^2$$
$A>0$，因此:
$$\frac{1}{A}(A{k}_1+B{k}_2{)}^2+\frac{1}{A}(x_1{y}_2-x_2{y}_1{)}^2{k}_2^2$$
至此，我们可以将式子划分为两部分,一部分只与$k_2$，一部分与$k_1,k_2$之间的最小距离有关。
由于$k_1,k_2$之间的最小距离在$\frac{1}{2}$范围内，因此可以计算出第一个部分的范围。
根据前文有:
$$-\frac{B}{A}{k}_2-\frac{1}{2}\le k_1\le -\frac{B}{A}{k}_2+\frac{1}{2}$$
$$-\frac{1}{2}A\le k_{1}A+B{k}_2\le \frac{1}{2}A$$
因此:$$(k_{1}A+B{k}_2{)}^2\le \frac{1}{4}{A}^2$$
我们还发现当$k_2>1,(x_1{y}_2-x_2{y}_1{)}^2>0$时,$(x_1{y}_2-x_2{y}_1{)}^2{k}_2^2$的增长速度实际上相当的快。
因此我们可以直接从$1$开始枚举$k_2$,由于后面的部分一定是递增的，前面的部分有一定的范围，因此当后面的部分与$k_2=1$取到的值已经相差$\frac{1}{4}{A}^2$时，我们就可以不再枚举了。
由于$\frac{1}{4}{A}^2\le 10^8$
因此这样做的时间复杂度为$O(T\sqrt{C})$的。
其中$C$小于$10^8$

解法2

见金斌的集训队论文《欧几里得算法的应用》-二维欧几里得算法
比第一种方法简单多了，虽然时间复杂度不明…

实现

这是解法2的实现，解法1的懒得写了。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
#define LL long long
void exgcd(LL x1,LL y1,LL x2,LL y2,LL &a,LL &b)
{
	LL x=x1*x2+y1*y2,y=x1*x1+y1*y1,z=x2*x2+y2*y2;
	if(x<0){
		exgcd(-x1,-y1,x2,y2,a,b),a=-a;
		return ;
	}
	if(4*x*x<y*z||x==0){
		if(y>z)b=1;
		else a=1;
		return ;
	}
	bool F=0;
	if(y<z)swap(x1,x2),swap(y1,y2),F=1;
	y=x1*x1+y1*y1,z=x2*x2+y2*y2;
	LL t=x/z,t1=t+1;
	if(x-t*z<=t1*z-x&&t)
		exgcd(x1-t*x2,y1-t*y2,x2,y2,a,b),b-=(t*a);
	else exgcd(x1-t1*x2,y1-t1*y2,x2,y2,a,b),b-=(t1*a);
	if(F)swap(a,b);
}
int main()
{
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
	LL x1,x2,y1,y2;
	while(~scanf("%lld%lld%lld%lld",&x1,&y1,&x2,&y2))
	{
		LL a=0,b=0;
		exgcd(x1,y1,x2,y2,a,b);
		LL res1=(x1*a+x2*b)*(x1*a+x2*b),res2=(y1*a+y2*b)*(y1*a+y2*b);
		printf("%lld\n",res1+res2);
	}
}