CodeForces 372 E.Drawing Circles is Fun（计算几何+dp）

Description
给出有n个点的点集S，要求从S中选出一个子集由一些点对(P[1],p[2]),(p[3],p[4]),…,(p[2k-1],p[2k])组成且满足下列三个条件：
1.k>=2
2.所有P[i]互不相同
3.对任意两个点对(P[2i-1],P[2i]),(p[2j-1],P[2j])，三角形OP[2i-1]P[2j-1]的外切圆与三角形OP[2i]P[2j]相切，三角形OP[2i-1]P[2j]的外切圆与三角形OP[2i]P[2j-1]相切
问满足条件的子集数
Input
第一行一整数n表示S中点数，之后路n行每行四个整数a[i],b[i],c[i],d[i]表示第i个点的横纵坐标为(a[i]/b[i],c[i]/d[i]) (1<=n<=1000,0<=|a[i]|,|c[i]|<=50,1<=|b[i]|,|d[i]|<=50,(a[i],c[i])!=(0,0))
Output
输出合法的子集个数，结果模1e9+7
Sample Input
10
-46 46 0 36
0 20 -24 48
-50 50 -49 49
-20 50 8 40
-15 30 14 28
4 10 -4 5
6 15 8 10
-20 50 -3 15
4 34 -16 34
16 34 2 17
Sample Output
2
Solution
以原点为反演中心进行反演变换，即(x,y)->(x/(x^2+y^2),y/(x^2+y^2))，那么一个以原点为圆心的圆就变成了一条不过原点的直线，两个相切的圆就变成了两条平行直线，如果OAB和OCD满足条件，那么说明AB与CD平行，AD和BC平行，即说明ABCD是一个平行四边形，我们把所有点对以（中点坐标，斜率）的形式存下来，只有中点坐标重复的点对才可以组成一个平行四边形，把这些点对排序，对于每一批中点重合的点对，设第i种有cnt[i]条（点可能重合进而点对可能重合），前i种共dp[i]种方案数，num[i]个点对，那么得到转移方程
num[i]=num[i-1]+cnt[i]
dp[i]=dp[i-1]+dp[i-1]*cnt[i]+num[i-1]*cnt[i]
dp数组转移的三项中，dp[i-1]表示不用当前这个点对的方案数，dp[i-1]*cnt[i]表示把前i-1种的所有方案后面加上当前这个点对的方案数，num[i-1]*cnt[i]表示从前面i-1种中拿出一个点对和当前这个点对组合形成的方案数
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;const int maxn=1111;const double eps=1e-8;const double PI=acos(-1.0);const ll mod=1000000007;int sign(double x){    if(abs(x)<eps)return 0;    if(x>eps)return 1;    return -1;}struct node{    double x,y,a;    bool operator <(const node &b)const    {        if(sign(x-b.x))return sign(x-b.x)<0;        if(sign(y-b.y))return sign(y-b.y)<0;        return sign(a-b.a)<0;    }}p[maxn*maxn];int n;double x[maxn],y[maxn];int main(){    while(~scanf("%d",&n))    {        for(int i=1;i<=n;i++)        {            double tx,ty;            scanf("%lf%lf",&tx,&ty);            x[i]=tx/ty;            scanf("%lf%lf",&tx,&ty);            y[i]=tx/ty;            double d=x[i]*x[i]+y[i]*y[i];            x[i]/=d,y[i]/=d;        }        int res=0;        for(int i=1;i<=n;i++)            for(int j=i+1;j<=n;j++)            {                p[res].x=0.5*(x[i]+x[j]),p[res].y=0.5*(y[i]+y[j]);                double temp=atan2(y[j]-y[i],x[j]-x[i]);                if(sign(temp)<0)temp+=PI;                if(sign(temp-PI)>=0)temp-=PI;                p[res++].a=temp;            }        sort(p,p+res);        ll ans=0,sum1,sum2;        for(int i=0;i<res;i++)        {            int j=i;            sum1=sum2=0;            while(j<res&&sign(p[j].x-p[i].x)==0&&sign(p[j].y-p[i].y)==0)            {                int k=j;                while(j<res&&!(p[k]<p[j]))j++;                sum2=(sum2+(sum1+sum2)*(j-k)%mod)%mod;                sum1=(sum1+j-k)%mod;            }            ans=(ans+sum2)%mod;            i=j-1;        }        printf("%I64d\n",ans);    }    return 0;}