Codeforces 559D Pick 定理 + 数学期望

题解

传送门 Codeforces 559D Randomizer

题解

已知 $n$ 个点构成的严格凸多边形，求其中可构成非退化的凸多边形的点集，所构成的多边形内部格点（整数坐标点）数量的数学期望。

满足条件的点集规模至少为 $3$。这样的点集构成的凸多边形可以看作是原始凸多边形切割数个凸多边形后得到。那么可以枚举 $O(n^2)$ 的被切割的凸多边形，计算不被包含在多边形内部格点数量的数学期望。假设所有原始凸多边形内部的格点都被点的子集构成的凸多边形包含，最后减去上述数学期望，即为所求。

被切割的凸多边形除了一条边以外，其余的边都是原始凸多边形上的位置连续边。那么可以固定端点 $i$，根据 Pick 定理，不断计算 $i,i+1,\cdots ,i+k$ 所构成的凸多边形上，属于原始凸多边形内部的格点。被切割的凸多边形的点集为 $i,i+1,\cdots ,i+k$ 时，满足条件的凸多边形数量，占总的子集构成的凸多边形数量为
$$\frac{2^{n-k-1}-\left(\genfrac{}{}{0ex}{}{n}{0}\right)}{2^n-\left(\genfrac{}{}{0ex}{}{n}{0}\right)-\left(\genfrac{}{}{0ex}{}{n}{1}\right)-\left(\genfrac{}{}{0ex}{}{n}{2}\right)}$$ 实际上，$2^{10^5}$ 超出了 long double 型的表示范围，而且 $O(n^2)$ 时间复杂度过大。考虑误差，取 $k$ 至多为 $64$ 即可；当 $n$ 过大时，$2^n$ 远大于 $\left(\genfrac{}{}{0ex}{}{n}{k}\right),0\le k\le 2$，那么将上述数量的比值取 $\frac{1}{2^{k+1}}$ 即可。

#include <bits/stdc++.h>
using namespace std;
typedef long double db;
typedef long long ll;
const int MAXN = 1E5 + 5;
struct P
{
    ll x, y;
    ll det(P p) { return x * p.y - p.x * y; }
    P operator-(P p) { return {x - p.x, y - p.y}; }
} ps[MAXN];
int N;
db pw[105], f[105];

ll gcd(ll a, ll b) { return b == 0 ? a : gcd(b, a % b); }

int add(int x, int d)
{
    x += d;
    return x >= N ? x - N : x;
}

db get(int k) { return N >= 100 ? pw[k + 1] : f[k]; }

int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> N;
    for (int i = 0; i < N; ++i)
        cin >> ps[i].x >> ps[i].y;
    db sum = 2;
    for (int i = 0; i < N; ++i)
    {
        int j = add(i, 1);
        ll x = ps[i].x - ps[j].x, y = ps[i].y - ps[j].y;
        sum += abs(ps[i].det(ps[j])) - gcd(abs(x), abs(y));
    }
    pw[0] = 1;
    for (int i = 0; i < 100; ++i)
        pw[i + 1] = pw[i] / 2;
    if (N < 100)
    {
        db down = pow((db)2, N) - 1 - N - (db)N * (N - 1) / 2;
        for (int i = 1; i < min(N, 64); ++i)
            f[i] = (pow((db)2, N - i - 1) - 1) / down;
    }
    db del = 0;
    for (int i = 0; i < N; ++i)
    {
        db t = 0;
        for (int k = 1; k < min(N, 64); ++k)
        {
            int u = add(i, k), v = add(i, k - 1);
            ll x = ps[u].x - ps[v].x, y = ps[u].y - ps[v].y;
            t += abs((ps[u] - ps[i]).det(ps[v] - ps[i])) - gcd(abs(x), abs(y));
            if (k > 1)
            {
                ll x2 = ps[u].x - ps[i].x, y2 = ps[u].y - ps[i].y;
                del += get(k) * (t + gcd(abs(x2), abs(y2)));
            }
        }
    }
    cout << fixed << setprecision(9) << (sum - del) / 2 << '\n';
    return 0;
}