Algorithm: Intersection of two conics

You can work by considering the pencil containing the two conics ($\lambda f(x,y)+(1-\lambda) g(x,y)=0$) and find its degenerate element, which is formed of two straight lines. Then intersecting the straight lines with one of the conics is straightforward.

Another approach is by eliminating $y$ to obtain a fourth degree polynomial equation.

Combine the two equations to eliminate $y^2$. Then you can express $y$ in terms of a rational function of $x$. $$ a x^2 + b y^2 + c xy + d x+e y + f = 0\\ a' x^2 + b' y^2 + c' xy + d' x+e' y + f' = 0\\ $$ $$ (ab'-a'b) x^2 + (cb'-c'b) xy + (db'-d'b) x+(eb'-e'b) y + (fb'-f'b) = 0,\\ y=-\frac{(ab'-a'b) x^2 +(db'-d'b) x+ (fb'-f'b)}{(cb'-c'b) x +(eb'-e'b)}=\frac{a''x^2+d''x+f''}{c''x+e''}. $$ Substituting in one of the original equations, you get a univariate polynomial.

$$(a x^2+ d x+ f)(c''x+e'')^2+ b (a''x^2+d''x+f'')^2 - (c x +e) (c''x+e'')(a''x^2+d''x+f'')= 0 .$$

This can be solved by closed formulas.

The book "Perspectives on Projective Geometry" by Jürgen Richter-Gebert contains a detailed description of both theory and algorithms concerning calculations with conics. The book is light of examples though.

In the algorithm for splitting a degenerate conic into two lines, there is missing a minus sign. Page 190 step 3 should be: beta = sqrt( - B_{i,i} ).

For an example see the answers Intersection of conics using matrix representation and Decomposition of a degenerate conic which use a slightly different algorithm, but is nevertheless very helpful.

I have a complete solution which works fine for me but should undergo a bit more testing. I made it because Geogebra is terrible at finding roots.

My approach will benefit only modestly from parallelization. I am not convinced matrices are an improvement, but for a graphics library it is surely worth checking. My experience is this: symbolically, they make the overall problem easier to write and organize, but numerically they make the problem sequence longer, and more prone to error. Specifically, the steps are sequential, there is no way to avoid solving a cubic equation, and the matrices are not direct transformations, but represent ordinary linear systems which have to be solved. Moreover, matrix inverse cannot in general be assumed.

Instead, I simply solve the fourth order equation directly. I focus on a solid, efficient calculation method for all roots (real and imaginary) of a third order polynomial. From there, the fourth order solution is basically free. This delimits the problem, allowing me to focus on identifying and mitigating the sources of numerical instability.

For example, consider the case that the two sections are touching (tangent). The author of Intersection of conics using matrix representation flags this as a potential source of instability. But represented as an ordinary polynomial, this problem is well defined: we have simply a double root: f(x) = 0 and f'(x) = 0.

I believe my method is stable, and I have spent time insuring that difficult cases are handled correctly, but I have not been exhaustive. We should test it more carefully for implementation in a library. If you are interested, let me know. I am happy to post it here.