China-Final2016K Justice Rains From Above

这个题题意就是给你一个圆锥的顶点和母线与轴的夹角，再给一堆二位平面上的点，问最多能用圆锥覆盖多少个点。

这个题和给一堆二维点，问用半径为R的圆最多能覆盖多少个点思路差不多。在二维情况下，我们总能让某个点在圆上，然后来卡其他的点，再用极角扫描线判最多能覆盖多少个点。所以对于圆锥，我们也可以枚举所有点，让这个点在圆锥上，然后再用类似的方法来卡其他的点。

先说一下下文的一些约定，用 O 表示圆锥的顶点，Alpha 表示圆锥母线与轴的夹角，P[i] (1 <= i <= N) 表示二位平面上的点，V[i] (1 <= i <= N) 表示向量 O->P[i] 的单位向量。

首先枚举 P[i] 作为圆锥上的点，我们以 O 为原点，V[i] 为Z 轴重新定义一个坐标系，再枚举其余的点 P[j]，如果 V[i] 和 V[j] 的夹角已经大于了2·Alpha，那么肯定不能用圆锥同时覆盖这两个点，直接跳过。接着再计算圆锥的轴在什么范围内能覆盖 P[j]，范围用轴的向量在重定系下在XOY 平面内的投影与X 轴的有向夹角的范围来表示。有了所有覆盖区间，我们只要再扫一遍看看被覆盖次数最多的区间被覆盖了多少次更新答案即可。

这个题时间卡的比较近，开始我懒得推公式，直接二分的范围就T了，虽然复杂度都是 O(N^2·logN) 的。后来推了一下好像也不是很麻烦，在CF上看了别人AC的代码，感觉大家的公式都不太一样，而且都看不懂……

其他的就是一些细节处理和卡精的问题了，下面给我我AC的代码。

#include <bits/stdc++.h>
#define EPS 1e-8
using namespace std;
const int MAXN = 1005;
const double PI = acos(-1.0);
inline int sgn(double x)
{
    if (x < -EPS) return -1;
    else if (x > EPS) return 1;
    else return 0;
}
inline double sqr(double x)
{
    return x * x;
}
inline double to_deg(double rad)
{
    return rad / PI * 180;
}
inline double to_rad(double deg)
{
    return deg / 180 * PI;
}
inline double to_std(double rad)
{
    while (rad < 0) rad += 2 * PI;
    while (rad >= 2 * PI) rad -= 2 * PI;
    return rad;
}
struct Point
{
    double x, y, z;
    Point() {}
    Point(double x, double y, double z) : x(x), y(y), z(z) {}
    inline double length() const
    {
        return sqrt(sqr(x) + sqr(y) + sqr(z));
    }
};
typedef Point Vector;
Vector operator + (Vector a, Vector b)
{
    return Vector(a.x + b.x, a.y + b.y, a.z + b.z);
}
Vector operator - (Vector a, Vector b)
{
    return Vector(a.x - b.x, a.y - b.y, a.z - b.z);
}
Vector operator * (Vector v, double mul)
{
    return Vector(v.x * mul, v.y * mul, v.z * mul);
}
Vector operator / (Vector v, double div)
{
    return Vector(v.x / div, v.y / div, v.z / div);
}
double dot(Vector a, Vector b)
{//点积
    return a.x * b.x + a.y * b.y + a.z * b.z;
}
Vector cross(Vector a, Vector b)
{//叉积
    return Vector(a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
}
double angle(Vector a, Vector b)
{//向量夹角
    return acos(dot(a, b) / a.length() / b.length());
}
Vector unit(Vector v)
{//单位向量
    return v / v.length();
}
typedef pair<double, double> Interval;
typedef pair<double, int> Event;
int n;
double alpha;
Point p[MAXN];
Vector v[MAXN];
int main()
{
    int t;
    scanf("%d", &t);
    for (int kase = 1; kase <= t; kase++)
    {
        scanf("%d%lf", &n, &alpha);
        alpha = to_rad(alpha);
        scanf("%lf%lf%lf", &p[0].x, &p[0].y, &p[0].z);
        for (int i = 1; i <= n; i ++)
        {
            scanf("%lf%lf", &p[i].x, &p[i].y);
            p[i].z = 0;
            v[i] = unit(p[i] - p[0]);
        }
        int ans = 1;
        for (int i = 1; i <= n; i ++)
        {
            v[0] = unit(Vector(rand(), rand(), rand()));
            Vector vx = unit(cross(v[i], v[0]));
            vector<Interval> itv;
            for (int j = 1; j <= n; j ++) if (j != i)
            {
                double theta = angle(v[i], v[j]) / 2;
                if (sgn(theta -  alpha) > 0) continue;
                double h = sqrt(sqr(cos(theta) / cos(alpha)) - 1) * cos(theta);
                Vector vm = unit(v[i] + v[j]) * cos(theta),
                    vn = unit(cross(v[i], v[j]));
                vm = vm - v[i] * dot(vm, v[i]);
                Vector vl = unit(vm - vn * h), vr = unit(vm + vn * h);
                double l = to_std(sgn(dot(v[i], cross(vx, vl))) * angle(vx, vl)) - EPS,
                    r = to_std(sgn(dot(v[i], cross(vx, vr))) * angle(vx, vr)) + EPS;
                if (l > r)
                {
                    itv.emplace_back(l, 2 * PI);
                    itv.emplace_back(0, r);
                }
                else itv.emplace_back(l, r);
            }
            vector<Event> evt;
            for (auto& j : itv)
            {
                evt.emplace_back(j.first, 1);
                evt.emplace_back(j.second, -1);
            }
            sort(evt.begin(), evt.end(),
            [](const Event& a, const Event& b)->bool
            {
                if (sgn(a.first - b.first)) return a.first < b.first;
                else return a.second > b.second;
            });
            int temp = 1;
            for (auto& j : evt)
            {
                temp += j.second;
                ans = max(ans, temp);
            }
        }
        printf("Case #%d: %d\n", kase, ans);
    }
    return 0;
}

欢迎各位dalao指点批评。