两圆求交点

/*
求圆（x1,y1,R）和圆（x2,y2,r）的交点
数组p存放圆心(1~n-1)和交点(n~)，数组radius存放半径 
count 用于记录交点数（包括圆心的数目） 
*/
for(int i=0; i<n; i++)
{
    // 初始化 
    double x1, y1, R;
    x1 = p[i].x, y1 = p[i].y,R = radius[i]; 
    double x2, y2, r;
    x2 = p[j].x, y2 = p[j].y,r = radius[j];
                                    
    double d = dist(x1,y1,x2,y2);//求两圆圆心距 
    if(d > R+r || d < f(R-r)) continue;//两圆相离时 
                                      
    if(x1 == x2 && y1 == y2) continue;//同心圆 
    else if(y1 == y2 && x1 != x2)
    {
         double a = ((R*R-r*r)-(x1*x1-x2*x2))/(2*x2-2*x1);
         if(d == f(R-r) || d == R+r)//只有一个交点时 
              p[count].x = a, p[count].y = y1,count ++;
              else//两个交点 
              {
                   double t = R*R-(a-x1)*(a-x1);
                   p[count].x = a, p[count].y = y1+sqrt(t),count ++;
                   p[count].x = a, p[count].y = y1-sqrt(t),count ++;
               }                                                                                
    }
    else if(y1 != y2)
    {
        double k ,disp;
        k = (2*x1-2*x2)/(2*y2-2*y1);
        disp = ((R*R-r*r)-(x1*x1-x2*x2)-(y1*y1-y2*y2))/(2*y2-2*y1);//直线偏移量 
        double a, b, c;
        a = (k*k+1);
        b = (2*(disp-y1)*k-2*x1);
        c = (disp-y1)*(disp-y1)-R*R+x1*x1;
        double disc;
        disc = b*b-4*a*c;//一元二次方程判别式 
        if(d == f(R-r) || d == R+r)                                      
        {
            p[count].x = (-b)/(2*a),p[count].y = k*p[count].x + disp,count++;
        }
        else
        {
            p[count].x = ((-b)+sqrt(disc))/(2*a),p[count].y = k*p[count].x + disp,count++;
            p[count].x = ((-b)-sqrt(disc))/(2*a),p[count].y = k*p[count].x + disp,count++;
        }
    }
}