I think I have managed to solve it. The results are near to perfect (sometimes inverted, but that's probably a different problem).
This is where I went wrong, and I hope I can help other people to not waste all the time I have wasted searching this.
I have based my code on the blinn-phong document. I had coordinates b0, b1, b2, b3. I used to view them as 2D coordinates with a w, but I have changed this view, and this solved the problem. By viewing them as 3D coordinates with z = 0, and making them homogenous 4D coordinates for transformation (w = 1), the solution arrived.
By calculating the C matrix: C = M3 * B, I got these new coordinates. When calculating the determinants d0, d1, d2, d3, I used to take the x, y coordinates from columns 0 and 1 in the C matrix, and the w factor from column 2. WRONG! When you think of it, the coordinates are actually 3D coordinates, so, for the w-factors, one should take column 3 and ignore column 2.
This gave me correct determinants, resulting in a discriminant that was able to sort out what type of curve I was handling.
But beware, what made my search even longer was the fact that I assumed that when it is visibly a serpentine, the result of the discriminant should always be > 0 (serpentine). But this is not always the case, when you have a mathematically perfect sepentine (coordinates are so that the mean is exact middle), the determinant will say it's a cusp (determinant = 0). I used to think that this result was wrong, but it isn't. So don't be fooled by this.