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:
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?
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.