2つの円の交点を求める

/*
  （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++;
        }
    }
}