Suppose we have two planar parametric curves f(t) and g(t), which are guaranteed to be second- or third-degree Bézier curves, which means they can be placed in "polynomial" form: $$ \textbf f(t) = \textbf a t^3 + \textbf b t^2 + \textbf c t + \textbf d\textrm{, where }\textbf a \neq (0, 0)\textrm{ or }\textbf b \neq (0, 0) $$ For a graphics application, I want to find, with certain precision, all the intersection points between f and g on their respective intervals. Those intersections are then used to subdivide the curves and place them in a modified doubly-connected edge list data structure to determine the topology of a certain path: I want to identify a path's "faces" and explicitly encode them in order to correctly implement the Loop-Blinn curve rendering technique (preferentially without using the stencil buffer).
One of the facts I can use is that, unless the curves are part of the same "matrix curve", the number of intersections is at most $\partial \textbf f \ \partial \textbf g$, since they are polynomial. That means 4 for quadratic-quadratic, 6 for cubic-quadratic and 9 for cubic-cubic curves. My initial ideas were to use bounding-box and subdivision techniques, but they fail and produce a large amount of intersections when the curves are tangent (meaning $\textbf f(t) = \textbf g(u)$ and ${\textbf f}'(t) \times {\textbf g}'(u) = 0$); trust me, I spent weeks trying to "improve" only this part of the algorithm, since it is the main one, given that missed or multiple intersections can break the rest of the code.
It wasn't until now that I decided to try more mathematical approaches, instead of more geometrical methods. We essentially have a function we want to find the roots of: $$\textbf d(u, v) = \textbf f(u) - \textbf g(v)$$ Given that function, we can use the multivariate Newton method in order to search for root: $$\textbf J(u,v) = \begin{bmatrix}f'_x(u) & -g'_x(v)\\f'_y(u) & -g'_y(v)\end{bmatrix}$$ $$\textbf J^{-1}(u,v) = -\frac 1{\textbf f'(u) \times \textbf g'(v)}\begin{bmatrix}f'_x(u) & -f'_y(u)\\ g'_x(v) & -g'_y(v)\end{bmatrix}$$
And then, choosing an iteration point $(u_0, v_0) \in [0, 1]^2$, we can use Newton method to find: $$\begin{bmatrix}u_{n+1}\\v_{n+1}\end{bmatrix} = \begin{bmatrix}u_n\\v_n\end{bmatrix} + \frac 1{\textbf f'(u_n) \times \textbf g'(v_n)}\begin{bmatrix}f'_x(u_n) & -f'_y(u_n)\\ g'_x(v_n) & -g'_y(v_n)\end{bmatrix} \left(\textbf f(u_n) - \textbf g(v_n)\right)$$
That raised some questions though: this method will also fail on tangent curves (we will have numerical instabilities because of $\textbf f'(u) \times \textbf g'(v)$ going to zero), and this has no guarantee to find but one of the roots, and then I would have to randomly test points until I find all the curves (I can fail fast if $(u_n, v_n)$ falls off the unit square). Another thought would try to find the zeros of $\left\|\textbf d(u,v)\right\|^2$, since it will only be zero when the curves intersect, but finding the zeros of a $\mathbb{R}^2$ to $\mathbb{R}$ sounds even more challenging; I could try to minimize it, but I would fall on solving a system of equations of fifth degree, which falls again to trying the roots numerically.
Another insight would be to detect tangent curves in a way; there's no way to escape those, because they are artistically pleasing and they will always appear between "consecutive" curves in a piecewise-smooth path (they are tangent at their endpoints!). So I would like some new insights on how to proceed to solve that. The results don't need to be infinitely precise, but they should at least be tractable enough so I can derive topology information from them. Thank you!
Doing this reliably is very difficult. That's why CAD companies have had hundreds of developers working for decades on problems like this. I'd say a couple of weeks to get something that works 95% of the time, and several months to get to 99%.
As you said, it's a good idea to look for tangent intersections first, and treat those as a special case. If you just treat a tangent intersection as a regular intersection, it's an ill-conditioned problem, so no algorithm is going to work reliably.
I'd guess that tangent intersections will be quite rare in the interior of curve segments, but very common at shared end-points, so you should special-case that, too. In fact, shared end-points (tangent or not) are likely to be a very common type of intersection, so you should find those first.
If you're going to resort to generic multivariate root-finding, then don't write that code yourself. Go find some code written by people whose expertise is far greater than yours or mine. The Numerical Recipes book is much maligned by true experts, but it's not a bad place to start. Certainly better than re-inventing root-finding yourself.
There's a good discussion of curve-curve intersection techniques in chapter 7 of this document, though it focuses more on performance than reliability.
This book might also be helpful. It has a pretty good discussion of polynomial root-finding and how this can be applied to intersection problems.