What's the fastest way to find the approximate coordinates of a circle tangent to two other circles in Python?

Question

I'm trying to write a piece of code that packs circles of differing radii efficiently into a square-shaped sheet. With the algorithm I'm using, I need to find the tangent between:

• a circle with a known radius r and unknown coordinates (a, b) (on an x-y plane)
• 2 other circles ri (i = 1, 2) with known radii r1 and r2 along with coordinates (c1.x, c1.y) and (c2.x, c2.y) respectively

Equation being used (Distance between circle centers (ri and r) = ri + r):

((c1.x-a)**2 + (c1.y-b)**2)**(0.5) = r1 + r
((c2.x-a)**2 + (c2.y-b)**2)**(0.5) = r2 + r

I need to be able to solve this system of equations fast as it will need to be computed around 250 times per iteration.

Currently I'm using the solve() function from the SymPy library, since the system of equations has all non-linear elements. Currently it takes about 10 seconds to execute each solve() function call (way too long :( ).

Here's how I have the solve() function running in my program: solve([((c1.x-a)**2 + (c1.y-b)**2)**(0.5) - c1.r - circle.r, ((c2.x-a)**2 + (c2.y-b)**2)**(0.5) - c2.r - circle.r], a,b) *note that a and b are set to real numbers only (real=True)

I've been scrolling through the internet trying to find some linearization methods of circles or some kind of polynomial approximation to speed up the solving process but have found nothing and gotten nowhere.

Is there a faster and somewhat accurate way to solve a system of two circle equations in Python?

Answer 1

It sounds like you are not really using SymPy in the intended way here. The idea with using SymPy's symbolic solve function is that you would use it to find a general formula and then substitute numeric values into the formula.

First find the general formula:

In [39]: x1, x2, y1, y2, a, b, r, r1, r2 = symbols('x1, x2, y1, y2, a, b, r, r1, r2')

In [40]: eqs = [Eq((x1 - a)**2 + (y1 - b)**2, (r1 + r)**2), Eq((x2 - a)**2 + (y2 - b)**2, (r2 + r)**2)]

In [41]: eqs[0]
Out[41]: 
         2            2           2
(-a + x₁)  + (-b + y₁)  = (r + r₁) 

In [42]: eqs[1]
Out[42]: 
         2            2           2
(-a + x₂)  + (-b + y₂)  = (r + r₂) 

In [43]: sol = solve(eqs, [a, b])

In [44]: print(sol)
[((-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 + 2*y2)*(-sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)) + (-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))))/(2*(x1 - x2)), -sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)) + (-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))), ((-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 + 2*y2)*(sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)) + (-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))))/(2*(x1 - x2)), sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)) + (-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(2*(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))]

Now you can use SymPy's lambdify function to turn the solution into something that can be evaluated efficiently numerically:

In [45]: f = lambdify([x1, y1, r1, x2, y2, r2, r], sol)

In [46]: f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0)
Out[46]: [(-1.0, 5.031128874149275), (-1.0, -3.0311288741492746)]

Let's time how fast this function is:

In [47]: %timeit f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0)
27.8 µs ± 765 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

So it takes 30 microseconds to compute the answer for a given value of the parameters x1, x2 etc. Doing that 250 times would take less than 10 milliseconds.

If you install llvmlite then you can compile this expression using LLVM which is faster. We have to flatten the solution list of tuples into a single list and use SymPy's cse function to prepare the input expression:

In [69]: sol_flat = [sol[0][0],sol[0][1],sol[1][0],sol[1][1]]

In [70]: from sympy.printing.llvmjitcode import llvm_callable

In [71]: f = llvm_callable([x1, y1, r1, x2, y2, r2, r], cse(sol_flat))

In [72]: f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0)
Out[72]: (-1.0, 5.031128874149275, -1.0, -3.0311288741492746)

In [73]: %timeit f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0)
1.49 µs ± 24.5 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

Now it takes 1.5 microseconds to evaluate for a single value of the parameters so 250 evaluations would take less than 1 millisecond.

If you are restarting the program many times then it is better to save the expression for the general solution rather than calling solve again each time your program runs. In fact if you only need to evaluate the function 250 times then it is probably better not to bother with lambdify or llvm_callable because just creating the fast callable will probably take longer than 250 calls to a simple Python function. You can use pycode to print out the expression as Python code:

In [88]: print(pycode(sol))
[((1/2)*(-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 + 2*y2)*(-1/2*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))/(x1 - x2), -1/2*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)), ((1/2)*(-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 + 2*y2)*((1/2)*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))/(x1 - x2), (1/2)*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))]

Now make a function in your Python script and just paste in that code for the solution expressions:

In [91]: def f(x1,y1,r1,x2,y2,r2,r):
    ...:     return [((1/2)*(-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 +
    ...: 2*y2)*(-1/2*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y
    ...: 2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2
    ...: *y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*
    ...: y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 +
    ...: x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/
    ...: (x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))/(x1 - x2), -1/2*math.sqrt(-(-r1**2 + 2*r1*
    ...: r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1*
    ...: *2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*
    ...: x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2
    ...:  - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x
    ...: 2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y
    ...: 1*y2 + y2**2)), ((1/2)*(-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y
    ...: 1 + 2*y2)*((1/2)*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y
    ...: 2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**
    ...: 2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*
    ...: r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*
    ...: y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2
    ...: **3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))/(x1 - x2), (1/2)*math.sqrt(-(-r1**2 +
    ...:  2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2
    ...:  - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**
    ...: 2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r
    ...: *r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*
    ...: y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2
    ...:  - 2*y1*y2 + y2**2))]
    ...: 

Calling this function takes 18 microseconds:

In [93]: import math

In [94]: f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0)
Out[94]: [(-1.0, 5.031128874149275), (-1.0, -3.0311288741492746)]

In [95]: %timeit f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0)
18.8 µs ± 437 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

That is not as fast as using LLVM but avoids the overhead of compiling the function at the start of your program. It takes around 250 milliseconds on this computer to import llvm_callable and then use it to compile the function but it makes the function f about 10x faster than the pycode output. Using llvm_callable would be faster then if you needed to evaluate the function more than 10000 times per process but otherwise it is faster just to use the pycode output pasted into your Python script since it avoids import/compilation overhead at runtime.

Another option for avoiding compiler overhead is to generate an extension module so you can use ahead of time compilation but I think that this is probably overkill in your case. The simplest option here is to use pycode which will be plenty fast enough if you only need to evaluate the function 250 times. Here is complete timings for running a program that calls the function 250 times:

$ cat t.py 
import math

def f(x1,y1,r1,x2,y2,r2,r):
    return [((1/2)*(-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 + 2*y2)*(-1/2*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))/(x1 - x2), -1/2*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)), ((1/2)*(-2*r*r1 + 2*r*r2 - r1**2 + r2**2 + x1**2 - x2**2 + y1**2 - y2**2 + (-2*y1 + 2*y2)*((1/2)*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)))/(x1 - x2), (1/2)*math.sqrt(-(-r1**2 + 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2)*(-4*r**2 - 4*r*r1 - 4*r*r2 - r1**2 - 2*r1*r2 - r2**2 + x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))*(x1 - x2)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2) + (1/2)*(-2*r*r1*y1 + 2*r*r1*y2 + 2*r*r2*y1 - 2*r*r2*y2 - r1**2*y1 + r1**2*y2 + r2**2*y1 - r2**2*y2 + x1**2*y1 + x1**2*y2 - 2*x1*x2*y1 - 2*x1*x2*y2 + x2**2*y1 + x2**2*y2 + y1**3 - y1**2*y2 - y1*y2**2 + y2**3)/(x1**2 - 2*x1*x2 + x2**2 + y1**2 - 2*y1*y2 + y2**2))]

sols = [f(-3.0, 1.0, 1.5, 1.0, 1.0, 1.5, 3.0) for _ in range(250)]
$ time python t.py 

real    0m0.042s
user    0m0.036s
sys 0m0.005s

This way takes about 40 milliseconds for a complete run of the process that computes the solution 250 times. Most of this time is actually not even spent calling the function f but just importing things at startup.

Answer 2

This solution is based on Tim Roberts's comment, except that the two original circles do not have to be touching, so the distance between their centres may be more than r1+r2. (It does, however, assume no overlapping).

Use the distances between centres to find the angle between the line C1C2 and, say, C1C (where C is the new circle). Then rotate the (normalised and then scaled) line C1C2 by this angle.

Note that there are usually TWO solutions (one on either side of the original circle pair), so you will still have to choose which one you want.

import math

def getCentre( C1x, C1y, r1, C2x, C2y, r2, R ):
    dsq = ( C2x - C1x ) ** 2 + ( C2y - C1y ) ** 2
    d = math.sqrt( dsq )                                   # distance between circle centres C1 and C2

    # Use cosine formula (in reverse) to get cosine of angle between line C1C2 and C1C
    cosTheta = ( ( r1 + R ) ** 2 + dsq - ( r2 + R ) ** 2 ) / ( 2 * d * ( r1 + R ) )

    if ( abs( cosTheta ) > 1 ):
       print( "No co-tangent circle possible, as circles are too far apart" )
       print( "Just return midpoint (twice)" )
       x1, y1 = ( C1x + C2x ) / 2, ( C1y + C2y ) / 2
       x2, y2 = x1, y1

    else:
       sinTheta = math.sqrt( 1.0 - cosTheta ** 2 )

       # There are two solutions, one rotated theta from the inter-circle line, the other rotated by minus theta
       x1 = C1x + ( r1 + R ) * (  cosTheta * ( C2x - C1x ) - sinTheta * ( C2y - C1y ) ) / d
       y1 = C1y + ( r1 + R ) * (  sinTheta * ( C2x - C1x ) + cosTheta * ( C2y - C1y ) ) / d
       x2 = C1x + ( r1 + R ) * (  cosTheta * ( C2x - C1x ) + sinTheta * ( C2y - C1y ) ) / d
       y2 = C1y + ( r1 + R ) * ( -sinTheta * ( C2x - C1x ) + cosTheta * ( C2y - C1y ) ) / d

    return x1, y1, x2, y2


C1x, C1y, r1  = -3.0, 1.0, 1.5
C2x, C2y, r2  =  1.0, 1.0, 1.5
R = 3
x1, y1, x2, y2 = getCentre( C1x, C1y, r1, C2x, C2y, r2, R )

print( "Solution 1, centre: ", x1, y1 )
print( "Solution 2, centre: ", x2, y2 )

Output:

Solution 1, centre:  -1.0 5.031128874149275
Solution 2, centre:  -1.0 -3.0311288741492746