Finding intersection points of two ellipses

I'm writing a basic 2D shape library in Python (primarily for manipulating SVG drawings), and I'm at a loss for how to efficiently calculate the intersection points of two ellipses.

Each ellipse is defined by the following variables (all floats):

c: center point (x, y)
hradius: "horizontal" radius
vradius: "vertical" radius
phi: rotation from coordinate system's x-axis to ellipse's horizontal axis

Ignoring when the ellipses are identical, there could be 0 through 4 intersection points (no intersection, tangent, partially overlapping, partially overlapping and internally tangent, and fully overlapping).

I've found a few potential solutions:

• SymPy geometry module - This basically just plugs the ellipse equations into SymPy's solver. I'm not sure whether this makes sense without already having the solver. (Incidentally, I would have used SymPy instead of rolling my own, but it performs horribly when dealing with crazy floats)
• How to detect if an ellipse intersects(collides with) a circle - This could probably be adapted for two ellipses, but I'm a little fuzzy on how to turn it into sensible code.
• How Ellipse to Ellipse intersection? - The library the answer references (CADEMIA) might have a good algorithm, but I can't even figure out if it's open source.
• Wikipedia: Intersecting Two Conics - I don't have enough of a grasp of linear algebra to understand this solution.

Any suggestions on how I should go about calculating the intersections? Speed (it might have to calculate a lot of intersections) and elegance are the primary criteria. Code would be fantastic, but even a good direction to go in would be helpful.

In math, you need to express the ellipses as bivariate quadratic equations, and solve it. I found a doucument. All the calculations are in the document, but it may take a while to implement it in Python.

So another method is to approximate the ellipses as polylines, and use shapely to find the intersections, here is the code:

import numpy as np
from shapely.geometry.polygon import LinearRing

def ellipse_polyline(ellipses, n=100):
    t = linspace(0, 2*np.pi, n, endpoint=False)
    st = np.sin(t)
    ct = np.cos(t)
    result = []
    for x0, y0, a, b, angle in ellipses:
        angle = np.deg2rad(angle)
        sa = np.sin(angle)
        ca = np.cos(angle)
        p = np.empty((n, 2))
        p[:, 0] = x0 + a * ca * ct - b * sa * st
        p[:, 1] = y0 + a * sa * ct + b * ca * st
        result.append(p)
    return result

def intersections(a, b):
    ea = LinearRing(a)
    eb = LinearRing(b)
    mp = ea.intersection(eb)

    x = [p.x for p in mp]
    y = [p.y for p in mp]
    return x, y

ellipses = [(1, 1, 2, 1, 45), (2, 0.5, 5, 1.5, -30)]
a, b = ellipse_polyline(ellipses)
x, y = intersections(a, b)
plot(x, y, "o")
plot(a[:,0], a[:,1])
plot(b[:,0], b[:,1])

and the output:

It takes about 1.5ms on my PC.