URAL1815 Farm in San Andreas(费马点，圆圆相交)

题目：

给你三个城市A,B,C，要以某点为中心F，分别向三个城市建立三条道路。长度分别为FA,FB,FC.总花费为aFA+bFB+cFC

求(aFA+bFB+cFC)min

解决：

a=b=c时，那么F点即为费马点。

当不等时为费马点的拓展。

此时可以求出AFB,BFC,CFA三个角，具体如下图所示

特殊地a=b=c 三个角均为120.

百科上的介绍为如下：

对于每一个角小于120°的三角形ABC的每一条边为底边，向外做正三角形，然后做这个正三角形的外接圆。费马点即使这些圆的交点。

物理学的解释：

费马问题可以用物理的方法来解决。将平面上所给的三个给定点钻出洞来，再设桑格绳子系在同一个点上面，而绳子的末端都绑一个重量为m的重物。假设摩擦力可以忽略不计，那么绳子会被拉紧，而绳子的最后会停在平面一点上方。可以证明，这个点就是三个顶点所对应的费马点。首先由于绳子长度是固定的，而绳子竖直下垂的部分越长，重物的位置也就越低。势能也就越低。在平衡状态时，势能达到最小值，也就是绳子垂直下垂的部分长度达到最大值，而绳子水平长度为FA+FB+FC，因此FA+FB+FC最小，也就是达到费马点。有力学平衡定理，可以得到各个夹角为120°

从上面推广费马点定义，我们可以得出该题求的就是费马点。

至于公式直接解出费马，表示不知道。

我设计了如下算法：

首先我们可以求角ABF 以及 角ABC,以及角CBF.该标记对应上图。

代码对应如下：

       angAB = acos(( n*n-l*l-m*m ) / (2*l*m));

       angAC = acos(( m*m-l*l-n*n ) / (2*l*n));

       angBC = acos(( l*l-m*m-n*n ) / (2*n*m));

这些角度与l,m,n也就是上面提到的aFA+bFB+cFC三个中的a,b,c.

注意到上面画红线部分。O1,O2,O3三个圆交与同一个点 B=D=E.

该点其实就是所求费马点。

首先费马点不一定存在

当(l,m,n)三条边不能组成三角形的时候需要特判。因为此时上面的公式不能求出角度了。

另者，当有两个点重合时，也就是其实不是三角形而是一条直线时，也需要特判，用定比分点解出答案。并且，求出的费马点可能在三角形外，此时，我们只要同时算出三个顶点的花费，取最小值即可

现在我们的问题就是严格的三角形内求可存在的费马点。

O1,O2,O3是这样求出来的。

比如O1.角ABC用公式得出 r = AC / sin(角ABC)/2 O1是以两个以A，C为圆心，r为半径的两个圆的交点。注意到此时有两个交点，当角ABC为钝角时 O1与（AC的对应点F）异侧，否则为同侧。当ABC为锐角时O1与（AC的对应点F）同侧，否则为异侧。这样我们只要求A，C两圆的交点，在经过这样的判断可解除唯一的O1.同理可以求出O2,O3.

再O1,O2求交点，求出两个交点，然后再判断与O3的距离。这样就可以得到唯一的点B。

该点即为费马点。可以知道该点满足那求出来的三个角的关系。此题精度要求较高。需要在调整下精度。

 

/*
	Author : Lit
	Data : July 17 1:00
*/

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

const long double eps = 1e-10;
const double pi = acos(-1.0);
struct point 
{
	long double x,y;
} fi,se,p0,p1,p2,third,ans;
struct Circle
{
	long double x,y;
	long double r;
} O,A,B,first,second;
int Test;
struct seg
{
	long double a,b,c;
} line;
bool bo;
long double angAB,angBC,angAC,l,m,n,r,fans,a,b;

long double dis( point A, point B){
	return sqrt( (B.x-A.x)*(B.x-A.x)+(B.y-A.y)*(B.y-A.y) );
}
long double cross(point A,point B,point C){
	return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x);
}

seg get_vertical( seg line , point ret)
{
	seg lineA;
	lineA.a = line.b; lineA.b = -line.a; 
	lineA.c = -(lineA.a*ret.x+lineA.b*ret.y);
	return lineA;
}
point get_cross( seg line , seg lineA)
{
	point ret;
	ret.x = -(line.c*lineA.b - lineA.c*line.b) / (line.a*lineA.b-line.b*lineA.a);
	ret.y = -(line.a*lineA.c - lineA.a*line.c) / (line.a*lineA.b-line.b*lineA.a);
	return ret;
}
void Segment_Circle(Circle O,seg line,point &fi,point &se)
{
	point ret,rett;
	ret.x = O.x; ret.y = O.y;
	seg lineA;
	
	lineA = get_vertical( line , ret );	
	 
	if ( fabs(line.a)<eps )  
	{
		ret.x = 1; ret.y = 0;
	} else
	if ( fabs(line.b)<eps )
	{
		ret.x = 0; ret.y = 1;
	} else
	{
		ret.x = 1; ret.y = -line.a/line.b;
	}
	long double dist = fabs( line.a*O.x + line.b*O.y +line.c ) / (sqrt( line.a*line.a + line.b*line.b ) );
//	cout << O.r*O.r - dist*dist << endl;
	dist = sqrt(O.r*O.r - dist*dist+eps);
	rett.x = ret.x / sqrt( ret.x*ret.x + ret.y*ret.y);
	rett.y = ret.y / sqrt( ret.x*ret.x + ret.y*ret.y);
	ret = rett;
	rett = get_cross(line,lineA);
	fi.x = rett.x+dist*ret.x; fi.y = rett.y+dist*ret.y;
	se.x = rett.x-dist*ret.x; se.y = rett.y-dist*ret.y;
}
void Circle_Circle(Circle A,Circle B,point &fi,point &se)
{
	seg lineA;
	lineA.a = 2*(B.x-A.x);
	lineA.b = 2*(B.y-A.y);
	lineA.c = -(A.r*A.r-A.x*A.x-A.y*A.y+B.x*B.x+B.y*B.y-B.r*B.r);
	Segment_Circle(A,lineA,fi,se);
}
void out(point A){
	printf("%.6lf %.6lf\n",(double)A.x,(double)A.y);
}

void out(Circle A){
	printf("%.6lf %.6lf %.6lf\n",(double)A.x,(double)A.y,(double)A.r);
}
bool same( point A,point B){
	return (fabs(A.x-B.x)<eps) && (fabs(A.y-B.y)<eps);
}

void work(){
	cin >> p0.x >> p0.y >> p1.x >> p1.y >> p2.x >> p2.y;
	cin >> l >> m >> n;
	bo = 1;
	if ( same(p0,p1) || same(p1,p2) || same(p0,p2) ){
		bo = 0; return;
	}
	if (l>=m+n) {
		ans = p0; return;
	} else 
	if (m>=l+n) {
		ans = p1; return;
	} else
	if (n>=l+m) {
		ans = p2; return;
	}
	angAB = acos(( n*n-l*l-m*m ) / (2*l*m));
	angAC = acos(( m*m-l*l-n*n ) / (2*l*n));
	angBC = acos(( l*l-m*m-n*n ) / (2*n*m));
/*
	cout << angAB/3.1415926*180 << endl;
	cout << angBC/3.1415926*180 << endl;
	cout << angAC/3.1415926*180 << endl;
*/
	/*
		first circle
	*/

	r = dis( p0,p1) / sin( angAB ) / 2;
	A.x = p0.x; A.y = p0.y; A.r = r;
	B.x = p1.x; B.y = p1.y; B.r = r;
/*
	fi.x = A.x; fi.y = A.y;
	se.x = B.x; se.y = B.y;
	cout << dis(fi,se) << endl;
*/
	Circle_Circle(A,B,fi,se);
	
	if ( ( cross(p0,p1,fi)*cross(p0,p1,p2)<0 ) == (angAB>pi/2) ) {
		first.x = fi.x; first.y = fi.y; first.r = r;
	}  else{
		first.x = se.x; first.y = se.y; first.r = r;
	}
	/*
		second circle
	*/
//	out(A); out(B);
//	out(first);
//	cout << angBC/3.1415926*180 << endl;
	r = dis(p1,p2) / sin( angBC ) / 2;
	A.x = p1.x; A.y = p1.y; A.r = r;
	B.x = p2.x; B.y = p2.y; B.r = r;
	Circle_Circle(A,B,fi,se);
	if ( ( cross(p1,p2,fi)*cross(p1,p2,p0)<0 ) == (angBC>pi/2) ) {
		second.x = fi.x; second.y = fi.y; second.r = r;
	} else{
		second.x = se.x; second.y = se.y; second.r = r;
	}
//	out(second);
	/*
		third circle 's O
	*/
	r = dis(p2,p0) / sin( angAC ) /2;
	A.x = p2.x; A.y = p2.y; A.r = r;
	B.x = p0.x; B.y = p0.y; B.r = r;
	Circle_Circle(A,B,fi,se);
	if ( ( cross(p2,p0,fi)*cross(p2,p0,p1)<0 ) == (angAC>pi/2) ) {
		third.x = fi.x; third.y = fi.y; 
	} else{
		third.x = se.x; third.y = se.y; 
	}
//	out(third);

	Circle_Circle(first,second,fi,se);
	if ( fabs( dis(fi,third) - r )< 1e-6 ) ans = fi; else ans = se;
//	out(ans);
	
}
int main()
{
//	freopen("in.txt","r",stdin);
//	freopen("out.txt","w",stdout);
	cin >> Test;
	while (Test--) {
			work();
			
			if (!bo){
				if ( same(p0,p1) ){
					a = l+m; b = n;
					ans.x = p0.x * b/(a+b) + p1.x *a/(a+b);
					ans.y = p0.y * b/(a+b) + p1.y *a/(a+b);
				} else
				if ( same(p1,p2) ){
					a = m+n; b = a;
					ans.x = p1.x * b/(a+b) + p0.x *a/(a+b);
					ans.y = p1.y * b/(a+b) + p0.y *a/(a+b);
				} else
				if ( same(p0,p2) ){
					a = l+n; b = m;
					ans.x = p0.x * b/(a+b) + p1.x *a/(a+b);
					ans.y = p0.y * b/(a+b) + p1.y *a/(a+b);
				}
			}
			
			
			fans = l*dis(p0,ans) + m*dis(p1,ans) + n*dis(p2,ans);
			fans = min( fans, m*dis(p1,p0) + n*dis(p2,p0));
			fans = min( fans, l*dis(p0,p1) + n*dis(p2,p1));
			fans = min( fans, l*dis(p0,p2) + m*dis(p1,p2));
		out(ans);
		printf("%.10lf\n",(double) fans);
	}
	system("pause");
	return 0;
}