vector <point> c_c(circle P, circle Q){
    vector <point> ans(2);
    double d = dist( P.c , Q.c );
    double m = (d + (( P.r*P.r - Q.r*Q.r )/d)/2);
    double n = (d - (( P.r*P.r - Q.r*Q.r )/d)/2);
    point x = ((P.c*n)+(Q.c*m))/(m+n);
    point dir = (Q.c - P.c);
    point orto = (point( -dir.y , dir.x ))/d;
    ans[0] = x + orto*( P.r*P.r - m*m );
    ans[1] = x - orto*( P.r*P.r - m*m );
}