In this case you could try to substitute the values in the exact solution. An example of this sort is given in pages 46-50 of Sage for power users (1).
Here is a minimal working example:
var('x, y, x0, y0, x1, y1, d0, d1')
# two circles
eq1 = (x-x0)^2 +(y-y0)^2 == d0^2
eq2 = (x-x1)^2 + (y-y1)^2 == d1^2
Let's ask Sage to solve exactly this system of equations:
%time solns = solve([eq1, eq2], x, y)
CPU times: user 7.28 s, sys: 8.01 ms, total: 7.28 s
Wall time: 7.28 s
It takes a couple of seconds, but it does a good job in finding the exact formulas,
sage: solns[0][0].rhs().variables()
(d0, d1, x0, x1, y0, y1)
Of course, the advantage is that the solutions obtained with this method are exact,
sage: solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1)
-1/8*sqrt(7) + 19/8
sage: _.n()
2.04428108611693
Let's see how much it takes to evaluate it numerically:
%timeit _ = solns[0][0].rhs().subs(x0=1, y0=2, x1=3, y1=4, d0=2, d1=1).n()
1000 loops, best of 3: 704 µs per loop
Now let's wrap the solution into a Cython function (1), (2):
%%cython
from math import sqrt
def sol1_x(double x0, double y0, double x1, double y1, double d0, double d1):
return -0.5*(d0**2*x0 - d1**2*x0 - x0**3 + x0*x1**2 - x1**3 - x0*y0**2 - (x0 + x1)*y1**2 - (d0**2 - d1**2 - x0**2 + y0**2)*x1 + (2*x0*y0 + 2*x1*y0 + sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1))*y1 - sqrt(-d0**4 - d1**4 + 2*d0**2*x0**2 - x0**4 + 4*x0*x1**3 - x1**4 - y0**4 + 4*y0*y1**3 - y1**4 + 2*(d0**2 + x0**2 + y0**2)*d1**2 + 2*(d0**2 + d1**2 - 3*x0**2 - y0**2)*x1**2 + 2*(d0**2 - x0**2)*y0**2 + 2*(d0**2 + d1**2 - x0**2 + 2*x0*x1 - x1**2 - 3*y0**2)*y1**2 - 4*(d0**2*x0 + d1**2*x0 - x0**3 - x0*y0**2)*x1 - 4*(d1**2*y0 + 2*x0*x1*y0 - x1**2*y0 - y0**3 + (d0**2 - x0**2)*y0)*y1)*y0)/(x0**2 - 2*x0*x1 + x1**2 + y0**2 - 2*y0*y1 + y1**2)
Then,
%timeit _ = sol1_x(1.,2.,3.,4.,2.,1.)
100000 loops, best of 3: 7.03 µs per loop
sol1_x(1.,2.,3.,4.,2.,1.)
2.0442810861169263
That's a factor 100x of faster!
(1) if you are in the Jupyter notebook, create a new Code cell and in the first line use %%cython (and not %cython).
(2) Two observations:
- replace the operator
^ with **, for example with the command str(solns[0][0].rhs()).replace('^', '**') - use 0.5 instead of 1/2
Finally, Sage does provide a function to do this automatically, (eg. operating with Sage objects defined in the global scope), it is cython_lambda. I wasn't able to apply it here.. strange.
If you have to evaluate thousands of values, then it will be even more efficient to work with arrays, see the Cython documentation for more info or ask here :)
Without passing through the exact solution obtained with solve, here is the same answer with SciPy's fsolve interface:
from scipy.optimize import fsolve
def circles(variables, *args):
x, y = variables
x0, y0, x1, y1, d0, d1 = args
eq1 = (x-x0)^2 +(y-y0)^2 - d0^2
eq2 = (x-x1)^2 + (y-y1)^2 - d1^2
return (eq1, eq2)
%timeit _ = fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
1000 loops, best of 3: 1.24 ms per loop
sage: fsolve(circles, (-2.0, 2.0), args=(1.,2.,3.,4.,2.,1.))
array([ 2.04428109, 3.70571891])
