题意:
给定\(n(n\le 100)\)个半径为\(1\)的圆柱体,假设太阳光是从无限远处来的平行于地面直线,圆柱会在地面上投射出无限长的矩形阴影。
宽度定义为这些矩形的宽度并,问宽度最大的太阳角度。
题解:
一开始以为是枚举两个圆心,然后垂直连线啥的,发现错了,但是网上找不到好的题解。。于是看std理解了思路,遂写此题解。
显然,我们有了太阳的角度,求宽度是很简单的,只需要逆时针旋转\(\theta\)度,然后计算\(y\)轴上的宽度并,排序+扫一遍就行。
但是如何确定角度呢。
注意到,假如圆的相交关系不变,宽度的变化跟角度的变化是连续的。
先给出这个夹角的一个特解,然后找到每一段连续相交的圆的\(y\)最小和最大的两个,那么
求个导
\[dw = \sum (A.x cos\theta - A.ysin\theta - B.x cos \theta + B.y sin \theta) d\theta \] \[dw = \sum (A.x - B.x) cos\theta - \sum (A.y - B.y)sin \theta) d\theta \]令导数为\(0\),得到\(tan\theta = \frac{\sum (A.y - B.y)}{\sum (A.x - B.x)}\),于是可以得到\(\theta\)。
那么如何得到所有的特解呢,我们发现,两个圆投影是否相交当且仅当在两圆公切线处发生改变,于是我们枚举所有公切线角度,排序后就能得到每个区间了,可以取区间角平分线作为特解。
代码比较冗长。
#include <bits/stdc++.h>
#define pt(x) cout << x << endl;
#define Mid ((l + r) / 2)
#define lson (rt << 1)
#define rson (rt << 1 | 1)
using namespace std;
const double Pi = acos(-1.0);
const double eps = 1e-11;
struct Point {
double x, y;
Point() : x(0), y(0) {}
Point(double x, double y) : x(x), y(y) {}
Point operator+(const Point &a)const{ return Point(x + a.x, y + a.y); }
Point operator-(const Point &a)const{ return Point(x - a.x, y - a.y); }
} a[109];
typedef Point Vector;
struct line{
Point a, b;
line() {}
line(const Point &a,const Point &b) : a(a), b(b) {}
};
struct circle{
Point c;
double r;
circle() : c(Point(0,0)), r(0) {}
circle(const Point &c, double r) : c(c), r(r) {}
};
Point operator*(double c, const Point &a){ return Point(c * a.x, c * a.y); }
double abs(const Point &a){ return sqrt(a.x * a.x + a.y * a.y); }
Point rot(const Point &a, double theta){ return Point(a.x * cos(theta) - a.y * sin(theta), a.x * sin(theta) + a.y * cos(theta)); }
double arg(const Point &a){
double t = atan2(a.y, a.x);
return t < 0 ? t + 2 * Pi:t;
}
int tangent(const circle &C1,const circle &C2,vector<line> &res){
double d = abs(C1.c - C2.c);
if(d < eps) return 0;
int c=0;
if(C1.r + C2.r < d - eps){
double t = acos((C1.r + C2.r) / d);
res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c), t), C2.c + rot(C2.r / d * (C1.c - C2.c), t)));
res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c),-t), C2.c + rot(C2.r / d * (C1.c - C2.c),-t)));
c += 2;
} else if(C1.r + C2.r < d + eps){
Point p = C1.c + C1.r / d * (C2.c - C1.c);
res.push_back(line(p, p + rot(C2.c - C1.c, Pi / 2)));
c++;
}
if(abs(C1.r - C2.r) < d - eps){
double t1 = acos((C1.r - C2.r) / d), t2 = Pi - t1;
res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c), t1), C2.c + rot(C2.r / d * (C1.c - C2.c),-t2)));
res.push_back(line(C1.c + rot(C1.r / d * (C2.c - C1.c),-t1), C2.c + rot(C2.r / d * (C1.c - C2.c), t2)));
c+=2;
} else if(abs(C1.r - C2.r) < d + eps){
Point p = C1.c + C1.r / d * (C2.c - C1.c);
res.push_back(line(p, p + rot(C2.c - C1.c, Pi / 2)));
c++;
}
return c;
}
int n;
double maxn, minn, sima, simi;
double calc(double theta) {
double tt[109] = {0};
for(int k = 1; k <= n; k++) {
tt[k] = rot(a[k], theta).y;
}
sort(tt + 1, tt + 1 + n);
double mr = 0, sum = 0;
sum = 2.0;
mr = tt[1] + 1.0;
for(int k = 1; k <= n; k++) {
if(tt[k] + 1.0 < mr) continue;
if(tt[k] - 1.0 < mr) {
sum += tt[k] + 1.0 - mr;
} else {
sum += 2.0;
}
mr = tt[k] + 1.0;
}
return sum;
}
double calc_peak(double phi) {
pair<double, int> p[109];
for(int i = 1; i <= n; i++) {
p[i] = make_pair(rot(a[i], phi).y, i);
}
sort(p + 1, p + 1 + n);
int pre = 1;
double x = 0, y = 0;
for(int i = 2; i <= n + 1; i++) {
if(i == n + 1 || p[i - 1].first + 1 < p[i].first - 1) {
x += (a[p[i - 1].second].x - a[p[pre].second].x);
y += (a[p[i - 1].second].y - a[p[pre].second].y);
pre = i;
}
}
double ans = arg({y, x});
if(ans > Pi) ans -= Pi;
return ans;
}
void work() {
if(n == 0) exit(0);
for(int i = 1; i <= n; i++)
scanf("%lf%lf", &a[i].x, &a[i].y);
maxn = -1e18;
minn = 1e18;
vector<line> v;
for(int i = 1; i <= n; i++) {
for(int j = i + 1; j <= n; j++) {
tangent(circle(a[i], 1.0), circle(a[j], 1.0), v);
}
}
vector<double> th;
for(int i = 0; i < v.size(); i++) {
line &l = v[i];
double tt = arg(l.b - l.a);
if(tt < Pi) tt = Pi - tt;
else tt = 2 * Pi - tt;
th.push_back(tt);
}
th.push_back(0);
th.push_back(Pi);
sort(th.begin(), th.end());
int tot = 0;
for(int i = 0; i < th.size(); i++)
if(i == 0 || th[i] > th[i - 1] + eps)
th[tot++] = th[i];
th.erase(th.begin() + tot, th.end());
for(int i = 0; i < tot; i++) {
double sum = calc(th[i]);
if(sum > maxn + eps) {
maxn = sum;
sima = th[i];
}
if(sum < minn - eps) {
minn = sum;
simi = th[i];
}
if(i - 1 >= 0) {
double phi0 = (th[i - 1] + th[i]) / 2;
double phi1 = calc_peak(phi0);
double sum = calc(phi1);
if(sum > maxn + eps) {
maxn = sum;
sima = phi1;
}
if(sum < minn - eps) {
minn = sum;
simi = phi1;
}
}
}
printf("%.15f\n%.15f\n", simi, sima);
}
signed main()
{
ios :: sync_with_stdio(0);
cin.tie(0);
while(~scanf("%d", &n)) work();
return 0;
}