cyCubicRoots yields incorrect results if coef[3] equals zero.

[paulhoux]

Hi Cem,

first of all: thank you for sharing this amazing piece of code with us.

I noticed that cy::CubicRoots contains a division by zero if coef[3] equals zero. One of the roots then becomes infinity and won't count towards the correct number of roots. It would be better to check for this situation and call cy::QuadraticRoots instead.

Steps to reproduce:

1. Construct a cubic curve using points P0 = Vec2(10, 10), P1 = Vec2(70, 10), P2 = Vec2(100, 40) and P3 = Vec2(100, 100).
2. Calculate the coefficients for a horizontal line through P0:

   float coef[4];
   coef[3] = p3.y - 3 * ( p2.y - p1.y ) - p0.y; // x^3
   coef[2] = 3 * ( p2.y - 2 * p1.y + p0.y );    // x^2
   coef[1] = 3 * ( p1.y - p0.y );               // x
   coef[0] = 0;                                 // 1

3. This yields the following coefficients: [ 0, 0, 90, 0 ]
4. In cy::CubicRoots, the value of delta_4 becomes 8100, but rv0 now contains a division by zero ftype rv0 = q / a.
5. The correct result would be a single root at 0.0.
6. If we call cy::QuadraticRoots with the same coefficients instead, it yields the proper result.

It's easy enough for me to do something like:

int numRoots =                                                        
!approxZero( coef[3] ) ? cy::CubicRoots( roots, coef, 0.0f, 1.0f ) : 
!approxZero( coef[2] ) ? cy::QuadraticRoots( roots, coef, 0.0f, 1.0f ) :
             cy::LinearRoot( roots[0], coef, 0.0f, 1.0f );

, but it would be better if cyCodeBase does this check for us.

Is this something you could change, or are there reasons not to?

[paulhoux]

Also note: when using intrinsics, cy::QuadraticRoots returns 2 (!) roots, both of which are 0. This is also incorrect.