题解 P1337 【[JSOI2004]平衡点 / 吊打XXX】
最近几天一直在怼计算几何……
这道题为什么要用模拟退火？明明是个：
二分！
假设我们已经确定了平衡点在下面的多边形内，那么在多边形内部选一个尽量靠近中间的点，然后确定受力方向：
然后我们发现，平衡点只能在黄色的新多边形内！
所以思路就出来了：先加上几个无穷远的点，然后不断用一个半平面去切割一个凸包得到一个新的凸包，进行二分，二分几十次后多边形会非常小，这时任意选一个顶点就是答案了，时间复杂度大概是
关键问题有两个：确定受力方向和半平面切割凸包。
每一个问题很简单，用每一个点和中心形成的单位向量乘质量，然后向量相加就行了。
对于第二个问题(convexcut)，我们可以顺次枚举原凸包上的点，如果这个点在半平面内部就加入到新的凸包，然后看这个点和下一个点形成的线段和半平面是否有交，如果有交就把交点加到新的凸包里面。
具体实现看代码：
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const double eps=1e-9;
int sign(double a){return a<-eps?-1:a>eps;}//带eps的实数比较
struct Point
{
double x,y;
Point(double xx=0,double yy=0){x=xx,y=yy;}
Point operator+(Point p){return Point(x+p.x,y+p.y);}
Point operator-(Point p){return Point(x-p.x,y-p.y);}
Point operator*(double d){return Point(x*d,y*d);}
Point operator/(double d){return Point(x/d,y/d);}
double dot(Point p){return x*p.x+y*p.y;}//点乘
double det(Point p){return x*p.y-y*p.x;}//叉乘
double abs2(){return x*x+y*y;}
double abs(){return sqrt(abs2());}
Point unit(){return sign(abs())==0?Point(0,0):(*this)/abs();}//单位向量，注意0向量的情况
Point rot90(){return Point(y,-x);}//旋转90度，用来把受力方向转变成半平面
}p[1005],last[1005],ans[1005];
int n,w[1005],top_last,top_ans;
double cross(Point p1,Point p2,Point q)
{
return (p2-p1).det(q-p1);
}
int crossOp(Point p1,Point p2,Point q){return sign(cross(p1,p2,q));}//用于判断q是否在p1p2的左侧
Point isLL(Point p1,Point p2,Point q1,Point q2)//直线交点
{
double x=(p2-p1).det(q2-p1);
double y=(q1-p1).det(p2-p1);
return q1+(q2-q1)/(x+y)*y;
}
void cut(Point q1,Point q2)//半平面切凸包
{
top_ans=0;
int n=top_last;
for(int i=1;i<=n;i++)
{
Point p1=last[i],p2=last[i%n+1];
int d1=crossOp(q1,q2,p1),d2=crossOp(q1,q2,p2);//看这两个点与半平面的位置关系
if(d1>=0)ans[++top_ans]=p1;
if(d1*d2<0)ans[++top_ans]=isLL(p1,p2,q1,q2);
}
for(int i=1;i<=top_ans;i++)last[i]=ans[i];
top_last=top_ans;
}
void work()
{
int m=top_last;
Point o(0,0);
for(int i=1;i<=m;i++)o=o+last[i];
o=o/m;//求的中心点，不知道有没有更优的取法，反正这样能AC
Point P(0,0);
for(int i=1;i<=n;i++)
P=P+(p[i]-o).unit()*w[i];//p是受力方向
if(P.abs()<=eps)//update
{
printf("%.3lf %.3lf\n",o.x,o.y);
exit(0);
}
P=o+P.rot90();
cut(o,P);
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
scanf("%lf%lf%d",&p[i].x,&p[i].y,&w[i]);
last[1]=Point(-10000,-10000);
last[2]=Point(10000,-10000);
last[3]=Point(10000,10000);
last[4]=Point(-10000,10000);
top_last=4;
for(int i=1;i<=60;i++)//二分
work();
printf("%.3lf %.3lf\n",last[1].x,last[1].y);
return 0;
}update：
一年过去了竟然又发现了一个锅（已经修正）
感谢评论区提醒。
我们注意到，每次需要求出受力方向。而如果受力刚好是
加一句特判就好了。如果受力等于