CodeForces 559 D.Randomizer（计算几何+期望）

本文链接：https://blog.csdn.net/V5ZSQ/article/details/79037813

Description

给出一个 $n$ 个点的凸包，等概率选则该凸包点集的大于等于三的子集形成一个新凸包，问该凸包内部整点的期望值

Input

第一行输入一整数 $n$ 表示点数，之后 $n$ 行每行输入两个整数 $x_i,y_i$ 表示第 $i$ 个点的坐标，按顺时针顺序给出这些点的坐标 $(3\le n\le 10^5,-10^9\le x_i,y_i\le 10^9)$

Output

输出子凸包内部整点的期望值，误差不超过 $10^{-9}$

Sample Input

4
0 0
2 0
2 2
0 2

Sample Output

0.2

Solution

首先要用到皮克定理：顶点为整点的凸多边形，其面积 $S$ =内部整点数 $+0.5*$ 边界上整点数 $-1$ ，故只要知道边界上点数和凸多边形面积即可得到凸多边形内部点数

问题求的是凸包内部整点数期望值，也即为原凸包内部整点数 $-$ 不被凸包包含的整点数期望值，如果子凸包中两个相邻的顶点在原凸包中分别是第 $i$ 个点和第 $j$ 个点，那么由 $i,i+1,...,j$ 这些点组成的凸包内部的整点以及第 $i$ 个点与第 $j$ 个点之间的连线上的整点（假设有 $num$ 个)都不被含有 $(i,j)$ 这条边的子凸包包含，假设 $i$ 到 $j$ 之间有 $k$ 个点，那么剩下 $n-k$ 个点中要选出至少一个点与 $i,j$ 组成子凸包，方案数 $2^{n-k}-1$ ，而总方案数为 $2^n-1-n-C_n^2$ ，那么 $(i,j)$ 这条边对不被凸包包含整点期望值的贡献为 $num\cdot \frac{2^{n-k}-1}{2^n-1-n-C_n^2}$ ，枚举 $i,j$ ，维护由第 $i,i+1,...,j$ 个点组成凸包的面积以及边界上点的数量以求出 $num$ 即可，注意到答案对精度要求不太大，而 $k$ 比较大时概率很小可以忽略不计，故只需要枚举 $k\le 35$ 的情况即可，且当 $n$ 比较大的时候，可以把概率 $\frac{2^{n-k}-1}{2^n-1-n-C_n^2}$ 近似为 $\frac{1}{2^k}$

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=100001;
struct node
{
    int x,y;
    node operator-(const node &b){return (node){x-b.x,y-b.y};}
    double operator*(const node &b){return 1.0*x*b.y-1.0*y*b.x;}
}a[maxn];
int gcd(int a,int b)
{
    return b?gcd(b,a%b):a;
}
double pick(double S,double T)
{
    return S+1.0-0.5*T;
}
double f[105];
double Get(int n,int k)
{
    if(n>100)return 1.0/f[k+1];
    return (f[n-k-1]-1)/(f[n]-1-n-0.5*n*(n-1));
}
int main()
{
    f[0]=1;
    for(int i=1;i<=100;i++)f[i]=2.0*f[i-1];
    int n;
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%d%d",&a[i].x,&a[i].y);
    double S=0,T=0;
    for(int i=0;i<n;i++)
    {
        int j=(i+1)%n;
        S+=(a[i]-a[0])*(a[j]-a[0])*0.5;
        T+=gcd(abs(a[j].x-a[i].x),abs(a[j].y-a[i].y));
    }
    double ans=0;
    for(int i=0;i<n;i++)
    {
        node pre=a[(i+1)%n]-a[i];
        double s=0,t=gcd(abs(pre.x),abs(pre.y));
        for(int j=(i+2)%n,k=2;k<=35&&(j+1)%n!=i;j=(j+1)%n,k++)
        {
            node now=a[j]-a[i];
            node temp=now-pre;
            double num=gcd(abs(now.x),abs(now.y));
            t+=num+gcd(abs(temp.x),abs(temp.y));
            s+=pre*now*0.5;
            ans+=(pick(s,t)+((j+1)%n==i?0:num-1))*Get(n,k);
            t-=num;
            pre=now;
        }
    }
    printf("%.10f\n",pick(S,T)-ans);
    return 0;
}