JZOJ【NOI2017模拟3.30】原谅

Description

终其一生，我们在寻找一个原谅。
犯下了太多错，要原谅的那个人，永远都是自己。
Samjia在深夜中望见了没有边界的人生，他没有想到过自己犯下了这么多的错误，他想在他的一生中寻求一个原谅。
他的人生是一个没有边界的平面，平面上有n个错误，每个错误是一个点，每个点i有一定的坐标(x[i],y[i])，有一个参数p 表示每个点有p的概率出现在平面上，注意两个不同的点的出现互相没有影响，Samjia可以在两个点之间连一条线段，两条线段不能在除了端点以外的地方相交，现在Samjia想知道他最多可以连的线段数的期望。
温馨提示：请看最后面的提示：）
本题的答案在mod 100000007意义下计算

Input

第一行两个正整数n和p表示平面上有n个点以及每个点出现概率为p
接下来n行第i行两个实数x[i]和y[i]表示一个点的坐标
保证不存在三点共线

Output

一行一个正整数表示Samjia最多可以连的线段数的期望

Sample Input

3 10000001
0 1
1 0
0 0

Sample Output

39000003
样例说明：
数据中的p为0.3
1×3×(0.3×0.3×0.7)+3×0.3×0.3×0.3=0.27
贡献为1的情况有三种
贡献为3的情况有一种

Data Constraint

对于20%的数据，p=1,1<=n<=5
对于40%的数据，p=1,1<=n<=1000
对于70%的数据，1<=n<=200
对于100%的数据，1<=n<=1000
坐标的绝对值小于等于10^4
保证0<=p<100000007

Hint
欧拉公式
在一个平面图内，设点数为V，边数为E，有界面数为F
一定满足：V+F-E=1

Source:

不愧是大爷出的题目，，感觉被虐惨了23333。。。。
比赛的时候看见欧拉公式感觉很可做，想了半天，移项发现V+F-1=E。
好啊，如果知道F不就直接知道E了？那我40分秒到手啊，认真搞一搞说不定切了。。
然后，F表示有界面数，要F最大，我就想到有一种构造凸包的方法叫做中心法，然后用中心法构造了一波凸包，过了大样例，感觉十分兹瓷，然而并没有什么卵用，一拍就错，再拍再错，，只能xjb打个10分的水法跑人。。

正解果然十分高端。。
如果概率为1，明显这题就是求n个点不相交最大线段数，那其实可以证明（我自己瞎想，不是题解证法）：
给定 n 个点时，平面图最多连 3n - 6 条边，设凸包大小为 x 以后，至少要在凸包外部连 x - 3 条曲线才能达到这个上界
那么只能连直线的时候答案上界就是 3n - 6 - (x - 3) = 3n - x - 3

题解说是显然，三角剖分证一波随手的事情（%%%）。。
然而光知道了这个结论并没有什么卵用，你要知道每个点在凸包上的概率，但是这个不好统计，此时我们可以发现，把点换成边会好做很多而且对结果无影响（在这道题目中），那么我们就可以看每条边在凸包上出现的概率，如何统计呢？用给出的概率p，那么每一条边在凸包上的概率不就是两个点出现的概率和其他不符合要求的点不出现的概率吗，那就是p^2*(1-p)。。

可能我这么说有点懵比，，那我就来引用一波专业的题解（by czy）
显然对于任意一个点集，三角剖分一定是最优的答案。
根据欧拉公式，平面点集的任意一个三角剖分的边数都可以由下式决定：

3V−d−3=E

其中，V是点数，d是凸包上边数，E是三角剖分边数。
那么我们显然可以分开来求这条式子的期望。
E[3V]=3E[V]=3np。求E[d]是经典问题：直接枚举一条边的起点，然后按照极角序（这里插一句，，一开始看别人用atan2一脸懵逼，不知道和atan有什么区别，结果换上去就炸了，atan在x,y趋近很大时会出现误差，而atan2不会）枚举终点，使用指针来维护某个方向上的点数，限定某个方向不能有点（就是我说的不符合要求点。。）来统计这条边出现在凸包上的概率。
//怎么我最近天天做鬼题。。（
code:

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const int N=2e3+5;
const int inf=1e9;
const int mo=1e8+7;
struct node
{
    double x,y;
}a[N],sta[N];
ll tot,cnt,n,p,f[N];
double cp(node a,node b)
{
    return a.x*b.y-a.y*b.x;
}
node operator +(node a,node b)
{
    node c;
    c.x=a.x+b.x;
    c.y=a.y+b.y;
    return c;
}
node operator -(node a,node b)
{
    node c;
    c.x=a.x-b.x;
    c.y=a.y-b.y;
    return c;
}
bool cmp(node a,node b)
{
    return atan2(a.y,a.x)<atan2(b.y,b.x);
}
int main()
{
    freopen("forgive.in","r",stdin);
    freopen("forgive.out","w",stdout);
    scanf("%lld%lld",&n,&p);
    fo(i,1,n)scanf("%lf%lf",&a[i].x,&a[i].y);
    f[0]=1;
    fo(i,1,n)
    f[i]=f[i-1]*(1-p+mo)%mo;
    ll ans=(3*n*p-3)%mo;
    fo(i,1,n)
    {
        tot=0;
        fo(j,1,n)
        if (i!=j)sta[++tot]=a[j]-a[i];
        sort(sta+1,sta+1+tot,cmp);
        ll k=1;
        cnt=0;
        fo(j,1,tot)
        {
            while (cp(sta[j],sta[k%tot+1])>0&&k%tot+1!=j)
            k=k%tot+1,cnt++;
            ans=(ans-f[n-2-cnt]*p%mo*p%mo+mo)%mo;
            if (cnt)cnt--;
            else k++;
        }
    }
    printf("%lld\n",(ans+3*f[n])%mo);
}