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;
}