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.