@
Gym101158 J: 求圆与给定凸多边形最大面积交
传送门：点我点我
求 $10 $ 个点组成的凸多边形 $(convexquad polygon) $ ，坐标范围 $[-100,100] $ ，与一个给定半径的圆的最大面积交。圆心的位置由你确定！
模拟退火坐标。
三分套三分求最优坐标，坑点是第二个三分的上下限需要动态求。
Steepest descent
Downhill simplex (Nelder–Mead)
Quasi-Newton (BFGS)
Evolution strategy (CMA-ES)
模拟退火
抄了个这场$Gymquad rank1 $的模拟退火方法。
#include<bits/stdc++.h>
#define fi first
#define se second
#define iis std::ios::sync_with_stdio(false)
#define pb push_back
#define o2(x) (x)*(x)
#define db double
using namespace std;
typedef long long LL;
typedef pair<int, LL> pii;
#define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))
#define crossOp(p1,p2,p3) sign(cross(p1,p2,p3))
const int MXN = 1e5 + 6;
const int INF = 0x3f3f3f3f;
const double eps = 1e-9;
const double PI = acos(-1.0);
int n, m;
struct P {
db x, y;
P(){};
P(db _x, db _y):x(_x),y(_y){}
P operator+(P p){return {x+p.x,y+p.y};}
P operator-(P p){return {x-p.x,y-p.y};}
db operator *(P s) {return x * s.x + y * s.y;}
db operator ^(P s) {return x * s.y - y * s.x;}
db dot(P p) {return x*p.x+y*p.y;}
db det(P p) {return x*p.y-y*p.x;}
};
inline int sign(db a) {return a < -eps?-1:a>eps;}
inline int cmp(db a, db b) {return sign(a-b);}
db dis_pp(P a, P b) {return sqrt((b - a)*(b - a));}//两点距离
db dot_p(P a, P b, P c) {return (b - a) * (c - a);}//点乘
double disToLine(P c,P A,P B){return fabs(cross(A,B,c))/dis_pp(A,B);}
//
int judge(P a, P b, P c) { //判断c是否在ab线段上（前提是c在直线ab上）
if(c.x >= min(a.x, b.x) && c.x <= max(a.x, b.x)
&& c.y >= min(a.y, b.y) && c.y <= max(a.y, b.y))return 1;
return 0;
}
db area(P b, P c, db r) {
P a(0.0, 0.0);
if (dis_pp(b, c) < eps) return 0.0;
db h = fabs(cross(a, b, c)) / dis_pp(b, c);
if(dis_pp(a, b) > r - eps && dis_pp(a, c) > r - eps) { //两个端点都在圆的外面则分为两种情况
db angle = acos(dot_p(a, b, c) / dis_pp(a, b) / dis_pp(a, c));
if(h > r - eps) {
return 0.5 * r * r * angle;
} else if(dot_p(b, a, c) > 0 && dot_p(c, a, b) > 0) {
db angle1 = 2 * acos(h / r);
return 0.5 * r * r * fabs(angle - angle1) + 0.5 * r * r * sin(angle1);
} else {
return 0.5 * r * r * angle;
}
} else if(dis_pp(a, b) < r + eps && dis_pp(a, c) < r + eps) { //两个端点都在圆内的情况
return 0.5 * fabs(cross(a, b, c));
} else { //一个端点在圆上一个端点在圆内的情况
if(dis_pp(a, b) > dis_pp(a, c)) { //默认b在圆内
swap(b, c);
}
if(fabs(dis_pp(a, b)) < eps) { //ab距离为0直接返回0
return 0.0;
}
if(dot_p(b, a, c) < eps) {
db angle1 = acos(h / dis_pp(a, b));
db angle2 = acos(h / r) - angle1;
db angle3 = acos(h / dis_pp(a, c)) - acos(h / r);
return 0.5 * dis_pp(a, b) * r * sin(angle2) + 0.5 * r * r * angle3;
} else {
db angle1 = acos(h / dis_pp(a, b));
db angle2 = acos(h / r);
db angle3 = acos(h / dis_pp(a, c)) - angle2;
return 0.5 * r * dis_pp(a, b) * sin(angle1 + angle2) + 0.5 * r * r * angle3;
}
}
}
P myp[MXN], ar[MXN];
db get_s(int n, db R, P O) { //求圆与多边形面积并
for (int i = 1; i <= n + 1; i++) myp[i] = ar[i];//顺或逆时针顺序
for(int i = 1; i <= n + 1; i++) myp[i] = myp[i] - O;
O = P(0, 0);
db sum = 0;
for(int i = 1; i <= n; i++) {
int j = i + 1;
db s = area(myp[i], myp[j], R);
if(cross(O, myp[i], myp[j]) > 0) sum += s;
else sum -= s;
}
return fabs(sum);
}
//
bool cmp1(P a, P b) { //3 4 1 2
double p = atan2(a.y, a.x), q = atan2(b.y, b.x);
if(p != q) return p < q;
return a.x < b.x;
}
bool cmp2(P a, P b) { //2 3 4 1 2
P c(0, 0);//尽量别单独使用
double tmp = cross(a, b, c);
if(tmp == 0) return a.x < b.x;
return tmp > 0;
}
int Quadrant1(P a) { //象限排序,注意x负半轴
if(a.x>0&&a.y>=0) return 3;
if(a.x<=0&&a.y>0) return 4;
if(a.x<0&&a.y<=0) return 1;
if(a.x>=0&&a.y<0) return 2;
return -1;
}
bool cmp5(P &a, P &b) {
int qa = Quadrant1(a), qb = Quadrant1(b);
if(qa == qb) return cmp2(a, b);
return qa < qb;
}
bool intersect(db l1, db r1, db l2, db r2) {
if(l1 > r1) swap(l1, r1);if(l2 > r2) swap(l2, r2);
return !(cmp(r1,l2)==-1||cmp(r2,l1)==-1);
}
bool isSS(P p1, P p2, P q1, P q2) {//判断线段相交
return intersect(p1.x,p2.x,q1.x,q2.x)&&intersect(p1.y,p2.y,q1.y,q2.y)&&
crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<=0&&crossOp(q1,q2,p1)*crossOp(q1,q2,p2)<=0;
}
bool isSS_strict(P p1, P p2, P q1, P q2) {//严格相交
return crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<0&&crossOp(q1,q2,p1)
*crossOp(q1,q2,p2) < 0;
}
bool isMiddle(db a, db m, db b) {
return sign(a-m)==0||sign(b-m)==0||(a<m!=b<m);
}
bool isMiddle(P a, P m, P b) {
return isMiddle(a.x,m.x,b.x)&&isMiddle(a.y,m.y,b.y);
}
bool onSeg(P p1, P p2, P q) {
return crossOp(p1,p2,q) == 0 && isMiddle(p1,q,p2);
}
bool ojbk(P p1, P p2, P q1, P q2) {
P p = p2 - p1;
P q = q2 - q1;
if(q.det(p) == 0) return 1;
return 0;
}
double rad(P p1, P p2) {
return atan2(p1.det(p2),p1.dot(p2));
}
bool xielv(P p1, P p2, P q1, P q2) {
P p = p2 - p1;
P q = q2 - q1;
if(q.det(p) == 0 && cmp(rad(p, q), 0) == 0) return 1;
return 0;
}
double _R, ans;
P cur(0, 0);
mt19937 lh(std::clock());
double Rand() {return 1.0*(lh()%INF)/INF;}
void simulateAnneal() {
ans = get_s(n, _R, cur);
double step = 30;
for(int iter = 0; iter < 2048; ++ iter) {
P best = cur;
db bestArea = ans;
bool improve = 0;
db dir = Rand();
const int K = 50;
db stepAngle = PI * 2/ K;
for(int i = 0; i < K; ++i) {
P nex{best.x+cos(dir)*step, best.y+sin(dir)*step};
db area = get_s(n, _R, nex);
if(area > bestArea) {
bestArea = area;
best = nex;
improve = 1;
}
dir += stepAngle;
}
if(improve) {
ans = max(ans, bestArea);
cur = best;
step *= 1.2;
}else {
step *= 0.9;
}
}
}
int main(){
scanf("%d%lf", &n, &_R);
double L = 1000, R = -1000;
for(int i = 1; i <= n; i++) {
scanf("%lf%lf", &ar[i].x, &ar[i].y);
L = min(L, ar[i].x); R = max(R, ar[i].x);
cur = cur + ar[i];
}
cur.x /= n, cur.y /= n;;
reverse(ar+1, ar+1+n);
ar[n + 1] = ar[1];
simulateAnneal();
printf("%.8f
", ans);
return 0;
}
三分套三分
先三分 $x $ ，再三分 $y $ 。
第一个三分的上下限是 最小横坐标 和 最大横坐标。
第二个三分的上下限是 由三分(x)出的这条直线 与 凸多边形交点的 最小纵坐标 和 最大纵坐标。
#include<bits/stdc++.h>
#define fi first
#define se second
#define iis std::ios::sync_with_stdio(false)
#define pb push_back
#define o2(x) (x)*(x)
#define db double
using namespace std;
typedef long long LL;
typedef pair<int, LL> pii;
#define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))
#define crossOp(p1,p2,p3) sign(cross(p1,p2,p3))
const int MXN = 1e5 + 6;
const int INF = 0x3f3f3f3f;
const double eps = 1e-9;
const double PI = acos(-1.0);
int n, m;
struct P {
db x, y;
P(){};
P(db _x, db _y):x(_x),y(_y){}
P operator+(P p){return {x+p.x,y+p.y};}
P operator-(P p){return {x-p.x,y-p.y};}
db operator *(P s) {return x * s.x + y * s.y;}
db operator ^(P s) {return x * s.y - y * s.x;}
db dot(P p) {return x*p.x+y*p.y;}
db det(P p) {return x*p.y-y*p.x;}
};
inline int sign(db a) {return a < -eps?-1:a>eps;}
inline int cmp(db a, db b) {return sign(a-b);}
db dis_pp(P a, P b) {return sqrt((b - a)*(b - a));}//两点距离
db dot_p(P a, P b, P c) {return (b - a) * (c - a);}//点乘
double disToLine(P c,P A,P B){return fabs(cross(A,B,c))/dis_pp(A,B);}
//
int judge(P a, P b, P c) { //判断c是否在ab线段上（前提是c在直线ab上）
if(c.x >= min(a.x, b.x) && c.x <= max(a.x, b.x)
&& c.y >= min(a.y, b.y) && c.y <= max(a.y, b.y))return 1;
return 0;
}
db area(P b, P c, db r) {
P a(0.0, 0.0);
if (dis_pp(b, c) < eps) return 0.0;
db h = fabs(cross(a, b, c)) / dis_pp(b, c);
if(dis_pp(a, b) > r - eps && dis_pp(a, c) > r - eps) { //两个端点都在圆的外面则分为两种情况
db angle = acos(dot_p(a, b, c) / dis_pp(a, b) / dis_pp(a, c));
if(h > r - eps) {
return 0.5 * r * r * angle;
} else if(dot_p(b, a, c) > 0 && dot_p(c, a, b) > 0) {
db angle1 = 2 * acos(h / r);
return 0.5 * r * r * fabs(angle - angle1) + 0.5 * r * r * sin(angle1);
} else {
return 0.5 * r * r * angle;
}
} else if(dis_pp(a, b) < r + eps && dis_pp(a, c) < r + eps) { //两个端点都在圆内的情况
return 0.5 * fabs(cross(a, b, c));
} else { //一个端点在圆上一个端点在圆内的情况
if(dis_pp(a, b) > dis_pp(a, c)) { //默认b在圆内
swap(b, c);
}
if(fabs(dis_pp(a, b)) < eps) { //ab距离为0直接返回0
return 0.0;
}
if(dot_p(b, a, c) < eps) {
db angle1 = acos(h / dis_pp(a, b));
db angle2 = acos(h / r) - angle1;
db angle3 = acos(h / dis_pp(a, c)) - acos(h / r);
return 0.5 * dis_pp(a, b) * r * sin(angle2) + 0.5 * r * r * angle3;
} else {
db angle1 = acos(h / dis_pp(a, b));
db angle2 = acos(h / r);
db angle3 = acos(h / dis_pp(a, c)) - angle2;
return 0.5 * r * dis_pp(a, b) * sin(angle1 + angle2) + 0.5 * r * r * angle3;
}
}
}
P myp[MXN], ar[MXN];
db get_s(int n, db R, P O) { //求圆与多边形面积并
for (int i = 1; i <= n + 1; i++) myp[i] = ar[i];//顺或逆时针顺序
for(int i = 1; i <= n + 1; i++) myp[i] = myp[i] - O;
O = P(0, 0);
db sum = 0;
for(int i = 1; i <= n; i++) {
int j = i + 1;
db s = area(myp[i], myp[j], R);
if(cross(O, myp[i], myp[j]) > 0) sum += s;
else sum -= s;
}
return fabs(sum);
}
//
bool cmp1(P a, P b) { //3 4 1 2
double p = atan2(a.y, a.x), q = atan2(b.y, b.x);
if(p != q) return p < q;
return a.x < b.x;
}
bool cmp2(P a, P b) { //2 3 4 1 2
P c(0, 0);//尽量别单独使用
double tmp = cross(a, b, c);
if(tmp == 0) return a.x < b.x;
return tmp > 0;
}
int Quadrant1(P a) { //象限排序,注意x负半轴
if(a.x>0&&a.y>=0) return 3;
if(a.x<=0&&a.y>0) return 4;
if(a.x<0&&a.y<=0) return 1;
if(a.x>=0&&a.y<0) return 2;
return -1;
}
bool cmp5(P &a, P &b) {
int qa = Quadrant1(a), qb = Quadrant1(b);
if(qa == qb) return cmp2(a, b);
return qa < qb;
}
bool intersect(db l1, db r1, db l2, db r2) {
if(l1 > r1) swap(l1, r1);if(l2 > r2) swap(l2, r2);
return !(cmp(r1,l2)==-1||cmp(r2,l1)==-1);
}
bool isSS(P p1, P p2, P q1, P q2) {//判断线段相交
return intersect(p1.x,p2.x,q1.x,q2.x)&&intersect(p1.y,p2.y,q1.y,q2.y)&&
crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<=0&&crossOp(q1,q2,p1)*crossOp(q1,q2,p2)<=0;
}
bool isSS_strict(P p1, P p2, P q1, P q2) {//严格相交
return crossOp(p1,p2,q1)*crossOp(p1,p2,q2)<0&&crossOp(q1,q2,p1)
*crossOp(q1,q2,p2) < 0;
}
bool isMiddle(db a, db m, db b) {
return sign(a-m)==0||sign(b-m)==0||(a<m!=b<m);
}
bool isMiddle(P a, P m, P b) {
return isMiddle(a.x,m.x,b.x)&&isMiddle(a.y,m.y,b.y);
}
bool onSeg(P p1, P p2, P q) {
return crossOp(p1,p2,q) == 0 && isMiddle(p1,q,p2);
}
bool ojbk(P p1, P p2, P q1, P q2) {
P p = p2 - p1;
P q = q2 - q1;
if(q.det(p) == 0) return 1;
return 0;
}
double rad(P p1, P p2) {
return atan2(p1.det(p2),p1.dot(p2));
}
bool xielv(P p1, P p2, P q1, P q2) {
P p = p2 - p1;
P q = q2 - q1;
if(q.det(p) == 0 && cmp(rad(p, q), 0) == 0) return 1;
return 0;
}
double _R;
double check1(double x) {
double L = 10000, R = -10000, midl, midr, tmpl, tmpr, ans = 0;
for(int i = 1; i <= n; ++i) {
if(sign((ar[i].x-x)*(ar[i+1].x-x)) <= 0) {
double tmp = ar[i].y + (ar[i+1].y-ar[i].y)/(ar[i+1].x-ar[i].x)*(x-ar[i].x);
L = min(L, tmp);
R = max(R, tmp);
}
}
for(int i = 0; i < 80; ++i) {
midl = (L + L + R) / 3;
midr = (L + R + R) / 3;
tmpl = get_s(n, _R, P{x, midl}); tmpr = get_s(n, _R, P{x, midr});
if(tmpl > tmpr) {
ans = max(ans, tmpl);
R = midr;
}else {
ans = max(ans, tmpr);
L = midl;
}
}
return ans;
}
int main(){
scanf("%d%lf", &n, &_R);
double L = 1000, R = -1000;
for(int i = 1; i <= n; i++) {
scanf("%lf%lf", &ar[i].x, &ar[i].y);
L = min(L, ar[i].x); R = max(R, ar[i].x);
}
reverse(ar+1, ar+1+n);
ar[n + 1] = ar[1];
double midl, midr, tmpl, tmpr, ans = 0;
for(int i = 0; i < 80; ++i) {
midl = (L + L + R) / 3;
midr = (L + R + R) / 3;
tmpl = check1(midl); tmpr = check1(midr);
if(tmpl > tmpr) {
ans = max(ans, tmpl);
R = midr;
}else {
ans = max(ans, tmpr);
L = midl;
}
}
printf("%.8f
", ans);
return 0;
}
模拟退火套路
for(int tc = 0; tc < 100; tc++) {//模拟退火次数
pdd o;//随机起点
double d = 400;//初始温度
double area;//最优解
while (d > eps) {//终止温度
pdd nex = o + rand();//随机移动
double tot;//当前答案
if (tot > area) o = nex, area = tot;//更新最优解
d *= 0.99;//降温系数0.998,0.975
}
ans = max(ans, area);//更新答案
}
double _R, ans;
P cur(0, 0);
mt19937 lh(std::clock());
double Rand() {return 1.0*(lh()%INF)/INF;}
void simulateAnneal() {
ans = get_s(n, _R, cur);
double step = 30;//看题目范围咯
for(int iter = 0; iter < 2048; ++ iter) {//或者改成>eps
P best = cur;
db bestArea = ans;
bool improve = 0;
db dir = Rand();
const int K = 50;//试探50个方向
db stepAngle = PI * 2 / K;
for(int i = 0; i < K; ++i) {
P nex{best.x+cos(dir)*step, best.y+sin(dir)*step};
db area = get_s(n, _R, nex);
if(area > bestArea) {
bestArea = area;
best = nex;
improve = 1;
}
dir += stepAngle;
}
if(improve) {
ans = max(ans, bestArea);
cur = best;
step *= 1.2;
}else {
step *= 0.9;
}
}
}