A quick way to determine the intersection points P1 and P2 is to take the vector w between the center points A and B of the circles.
When now imagining a rectangle that is spanned by two vectors a and b between point A and P1, we can say
P1 = A + a + b
P2 = A + a - b
The only remaining question is, how long the vectors a and b are:
We know that |a|^2 + |b|^2 = r_A^2
and (|w| - |a|)^2 + |b|^2 = r_B^2
, setting them equal and solving for |a|
yields
|a| = (r_A^2 - r_b^2 + |w|^2) / (2|w|)
|b| = |b| = +- sqrt(r_A^2 - |a|^2)
Now the lengths of the vectors can be used to construct the actual vectors a and b by using the normalized vector w.
When implementing this idea, the solution is pretty straightforward
w = {
x: B.x - A.x,
y: B.y - A.y
}
d = hypot(w.x, w.y)
if (d <= A.r + B.r && abs(B.r - A.r) <= d) {
w.x/= d;
w.y/= d;
a = (A.r * A.r - B.r * B.r + d * d) / (2 * d);
b = Math.sqrt(A.r * A.r - a * a);
P1 = {
x: A.x + a * w.x - b * w.y,
y: A.y + a * w.y + b * w.x
}
P2 = {
x: A.x + a * w.x + b * w.y,
y: A.y + a * w.y - b * w.x
}
} else {
P1 = P2 = null
}