本文深入探讨了费马点问题，一种经典的数学优化问题，旨在寻找一个点到三角形三个顶点距离之和最小的位置。文章详细介绍了费马点的几何构造方法和物理解释，并提供了一种算法实现，通过求解特定角度和利用外接圆交点来精确确定费马点位置，特别关注了边界情况和精度调整。
题目：
给你三个城市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;
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
