We are going to calculate points of intersections for 2 circels, given their centers and radiuses on a 2D plane. Say, we will implement the following C function
int intersect(struct circle_t const circles[], struct point_t intersections[]);
where point_t and circle_t are

struct point_t {
    double x;
    double y;
};

struct circle_t {
    point_t center;
    double r;
};


    The function intersect takes 2 circles as parameter, returns the number of points of intersection, and the details of the points will be put in parameter intersection. Since 2 circles may have up to 2 points of intersection, we may call that function in following way:

#include <stdio.h>

int main(void)
{
    struct circle_t circles[2]; // 2 circles
    struct point_t points[2]; // and up to 2 points

    /* here we read the information via stdin
     * in this format
     *   (x0, y0, r0)
     *   (x1, y1, r1)
     */
    scanf("%lf%lf%lf%lf%lf%lf",
          &circles[0].center.x, &circles[0].center.y, &circles[0].r,
          &circles[1].center.x, &circles[1].center.y, &circles[1].r);

    /* 2 circles are the same */
    /* NOTE: since x, y and r are all of type double,
     *       use == to compare them may cause problem.
     */
    if (circles[0].center.x == circles[1].center.x
     && circles[0].center.y == circles[1].center.y
     && circles[0].r == circles[1].r)
    {
       puts("The circles are the same.");
       return 0;
    }

    /* we call //intersect// here.
     * both //circles// and //points// are an array of 2.
     * if only one point of intersection, it will be stored in //points[0]//.
     */
    switch (intersect(circles, points)) {
        case 0:
            puts("No intersection.");
            break;
        case 1:
            printf("(%.3lf %.3lf)\n", points[0].x, points[0].y);
            break;
        case 2:
            printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",
                   points[0].x, points[0].y,
                   points[1].x, points[1].y);
    }
    return 0;
}


First, let's consider a simple situation: there is no intersection

#include <math.h> // we need //sqrt// and //fabs//

double distance_sqr(struct point_t const* a, struct point_t const* b)
{
    return (a->x - b->x) * (a->x - b->x) + (a->y - b->y) * (a->y - b->y);
}

double distance(struct point_t const* a, struct point_t const* b)
{
    return sqrt(distance_sqr(a, b));
}

int intersect(struct circle_t const circles[], struct point_t intersections[])
{
    double d = distance(&circles[0].center, &circles[1].center);
    if (d > circles[0].r + circles[1].r
     || d < fabs(circles[0].r - circles[1].r))
    {
        return 0;
    }
    /* more codes here */
}


To solute the rest of this problem, besides the usual way of representing circles like

• (x - x₀)² + (y - y₀)² = r², where (x₀, y₀) is the center and r is the radius

we will also use parametric equation for that

• x = r * cosθ + x₀
• y = r * sinθ + y₀

Now, let (x_a, y_a) be one point of intersections, then for circle₀ (that is (x₀, y₀) and r₀), there is some θ_a that meets

• x_a = r₀ * cosθ_a + x₀
• y_a = r₀ * sinθ_a + y₀

then substitute (x_a, y_a) in the equation of circle₁
(r₀ * cosθ_a + x₀ - x₁)² + (r₀ * sinθ_a + y₀ - y₁)² = r₁²

Don't panic though the equation is full of horrible calls to trigonometrical function. We do some tricks on the left and will get
  //LEFT// = ((r₀ * cosθ_a) + (x₀ - x₁))² + ((r₀ * sinθ_a) + (y₀ - y₁))²
      = (r₀ * cosθ_a)² + (x₀ - x₁)² + 2 * (r₀ * cosθ_a) * (x₀ - x₁) +
        (r₀ * sinθ_a)² + (y₀ - y₁)² + 2 * (r₀ * sinθ_a) * (y₀ - y₁)
      = (r₀ * cosθ_a)² + (r₀ * sinθ_a)² + (x₀ - x₁)² + (y₀ - y₁)² +
        2 * r₀ * (x₀ - x₁) * cosθ_a + 2 * r₀ * (y₀ - y₁) * sinθ_a
      = r₀² + (x₀ - x₁)² + (y₀ - y₁)² +
        2 * r₀ * (x₀ - x₁) * cosθ_a + 2 * r₀ * (y₀ - y₁) * sinθ_a
      = r₁² = //RIGHT//

Then let

• a = 2 * r₀ * (x₀ - x₁)
• b = 2 * r₀ * (y₀ - x₁)
• c = r₁² - (r₀² + (x₀ - x₁)² + (y₀ - y₁)²)

the equation will be transformed to
    a * cosθ_a + b * sinθ_a = c
and then represent sinθ_a as ±(1 - (cosθ_a)²)^1/2
    b * (±(1 - (cosθ_a)²))^1/2 = c - a * cosθ_a
square both sides (then the annoying ± is gone)
    b² * (1 - (cosθ_a)²) = c² + a² * (cosθ_a)² - 2 * c * a * cosθ_a
 => b² - b² * (cosθ_a)² = c² + a² * (cosθ_a)² - 2 * c * a * cosθ_a
 => (a² + b²) * (cosθ_a)² - (2 * a * c) * cosθ_a + (c² - b²) = 0

See, we are about to solute a quadratic equation of cosθ_a, which won't be very hard
    cosθ = (±(q² - 4 * p * s)^1/2 - q) / (2 * p)
where

• p = a² + b²
• q = -2 * a * c
• s = c² - b²

You may notice the sign ±, so how may we determine which it should really be? We have 2 choices here, one is to calculate for more details, until we find whether it should be positive or negative, the other is, we just //verify// the solution, by using the magical computing power of the computer. In the latter way, what we will do is substitute the following 4 solutions

• (((q² - 4 * p * s)^1/2 - q) / (2 * p), (1 - (cosθ_a)²)^1/2)
• ((-(q² - 4 * p * s)^1/2 - q) / (2 * p), (1 - (cosθ_a)²)^1/2)
• (((q² - 4 * p * s)^1/2 - q) / (2 * p), -(1 - (cosθ_a)²)^1/2)
• ((-(q² - 4 * p * s)^1/2 - q) / (2 * p), -(1 - (cosθ_a)²)^1/2)

in the origin equation, to see which will meet.
By the way, if there is only one point of intersection, (q² - 4 * p * s)^1/2 will be 0 and the sign is not important.

Here is the complete code

#include <stdio.h>
#include <math.h>

struct point_t {
    double x;
    double y;
};

struct circle_t {
    struct point_t center;
    double r;
};

double distance_sqr(struct point_t const* a, struct point_t const* b)
{
    return (a->x - b->x) * (a->x - b->x) + (a->y - b->y) * (a->y - b->y);
}

double distance(struct point_t const* a, struct point_t const* b)
{
    return sqrt(distance_sqr(a, b));
}

int intersect(struct circle_t const circles[], struct point_t intersections[])
{
    double a, b, c, p, q, s; // we have talked about those a, b, c, p, q, s
    double cos_value[2], sin_value[2]; // values for cosθa sinθa
    double d // distance between centers of the two circles
             = distance(&circles[0].center, &circles[1].center);

    if (d > circles[0].r + circles[1].r
        || d < fabs(circles[0].r - circles[1].r))
    {
        return 0;
    }

    a = 2.0 * circles[0].r * (circles[0].center.x - circles[1].center.x);
    b = 2.0 * circles[0].r * (circles[0].center.y - circles[1].center.y);
    c = circles[1].r * circles[1].r - circles[0].r * circles[0].r
        - distance_sqr(&circles[0].center, &circles[1].center);
    p = a * a + b * b;
    q = -2.0 * a * c;

    // there is only one point of intersection
    if (d == circles[0].r + circles[1].r
     || d == fabs(circles[0].r - circles[1].r))
    {
        cos_value[0] = -q / p / 2.0;
        sin_value[0] = sqrt(1 - cos_value[0] * cos_value[0]);

        intersections[0].x = circles[0].r * cos_value[0] + circles[0].center.x;
        intersections[0].y = circles[0].r * sin_value[0] + circles[0].center.y;

        // verify the solution here, if it is not, negate sinθa
        if (distance_sqr(&intersections[0], &circles[1].center)
            != circles[1].r * circles[1].r)
        {
            intersections[0].y = circles[0].center.y
                                 - circles[0].r * sin_value[0];
        }
        return 1;
    }

    s = c * c - b * b;
    cos_value[0] = (sqrt(q * q - 4.0 * p * s) - q) / p / 2.0;
    cos_value[1] = (-sqrt(q * q - 4.0 * p * s) - q) / p / 2.0;
    sin_value[0] = sqrt(1 - cos_value[0] * cos_value[0]);
    sin_value[1] = sqrt(1 - cos_value[1] * cos_value[1]);

    intersections[0].x = circles[0].r * cos_value[0] + circles[0].center.x;
    intersections[1].x = circles[0].r * cos_value[1] + circles[0].center.x;
    intersections[0].y = circles[0].r * sin_value[0] + circles[0].center.y;
    intersections[1].y = circles[0].r * sin_value[1] + circles[0].center.y;

    // verify and correct both solutions, still by negating sinθa
    if (distance_sqr(&intersections[0], &circles[1].center)
        != circles[1].r * circles[1].r)
    {
        intersections[0].y = circles[0].center.y - circles[0].r * sin_value[0];
    }
    if (distance_sqr(&intersections[1], &circles[1].center)
        != circles[1].r * circles[1].r)
    {
        intersections[1].y = circles[0].center.y - circles[0].r * sin_value[1];
    }

    /* if we got 2 same solutions, pick one and negate its Y coord
     * it is because cosθa is used to represent X coord of the solutions,
     * so when we got same solutions the problem must be on the Y coord.
     */
    if (intersections[0].y == intersections[1].y
     && intersections[0].x == intersections[1].x)
    {
        intersections[1].y = -intersections[1].y;
    }
    return 2;
}

int main(void)
{
    struct circle_t circles[2];
    struct point_t points[2];

    scanf("%lf%lf%lf%lf%lf%lf",
          &circles[0].center.x, &circles[0].center.y, &circles[0].r,
          &circles[1].center.x, &circles[1].center.y, &circles[1].r);

    if (circles[0].center.x == circles[1].center.x
     && circles[0].center.y == circles[1].center.y
     && circles[0].r == circles[1].r)
    {
       puts("The circles are the same.");
       return 0;
    }

    switch (intersect(circles, points)) {
        case 0:
            puts("No intersection.");
            break;
        case 1:
            printf("(%.3lf %.3lf)\n", points[0].x, points[0].y);
            break;
        case 2:
            printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",
                   points[0].x, points[0].y,
                   points[1].x, points[1].y);
    }
    return 0;
}