How to fix problem in serpentine case using Blinn/Loop algorithm

I'm trying to implement Blinn/Loop algorithm to render cubic bezier on shader. I calculate the curve type using code below (v0, v1, v2, v3 are cubic bezier segment):

CurveType CurveRenderer::DetermineType (Point    v0,
                                        Point    v1,
                                        Point    v2,
                                        Point    v3,
                                        double& d0, double& d1,
                                        double& d2, double& d3)
{

    CurveType curve_type;

    Eigen::MatrixXd M3(4,4);
    M3 << 1,  0,  0, 0,
         -3,  3,  0, 0,
          3, -6,  3, 0,
         -1,  3, -3, 1;

    Eigen::MatrixXd B(4, 4);
    B << v0.x, v0.y, 0, 1,
         v1.x, v1.y, 0, 1,
         v2.x, v2.y, 0, 1,
         v3.x, v3.y, 0, 1;

    Eigen::MatrixXd  C = M3 * B;

    Eigen::MatrixXd m0(3,3);
    m0 << C(3,0), C(3,1), C(3,3),
          C(2,0), C(2,1), C(2,3),
          C(1,0), C(1,1), C(1,3);

    Eigen::MatrixXd m1(3,3);
    m1 << C(3,0), C(3,1), C(3,3),
          C(2,0), C(2,1), C(2,3),
          C(0,0), C(0,1), C(0,3);

    Eigen::MatrixXd m2(3,3);
    m2 << C(3,0), C(3,1), C(3,3),
          C(1,0), C(1,1), C(1,3),
          C(0,0), C(0,1), C(0,3);

    Eigen::MatrixXd m3(3,3);
    m3 << C(2,0), C(2,1), C(2,3),
          C(1,0), C(1,1), C(1,3),
          C(0,0), C(0,1), C(0,3);

    d0 =  m0.determinant();
    d1 = -m1.determinant();
    d2 =  m2.determinant();
    d3 = -m3.determinant();

    double D = 3.0 * d2 * d2 - 4.0 * d1 * d3;

    curve_type = CURVE_TYPE_UNKNOWN;

    if(d1 != 0 && D > 0)                    curve_type = CURVE_TYPE_SERPENTINE;
    else if(d1 != 0 && D < 0)               curve_type = CURVE_TYPE_LOOP;
    else if(d1 != 0 && D == 0)              curve_type = CURVE_TYPE_CUSP_WITH_INFLECTION_AT_INFINITY;
    else if(d1 == 1 && d2 != 0)             curve_type = CURVE_TYPE_CUPS_WITH_CUPS_AT_INFINITY;
    else if(d1 == 0 && d2 == 0 && d3 != 0)  curve_type = CURVE_TYPE_QUADRATIC;
    else if(d1 == 0 && d2 == 0 && d3 == 0)  curve_type = CURVE_TYPE_LINE;

    return curve_type;
}

Where finally I obtain klm from matrix rMat (4 by 4), which is rMat = inverse M3 * F:

...
Eigen::MatrixXd rMat = invM3 * F;
...

To calculate Matrix F, I use this code:

...
case CURVE_TYPE_SERPENTINE:             
    tl = d2 + ((1.0 / sqrt(3.0)) * sqrt(3.0 * d2 * d2 - 4.0 * d1 * d3));
    sl = 2.0 * d1;              
    tm = d2 - ((1.0 / sqrt(3.0)) * sqrt(3.0 * d2 * d2 - 4.0 * d1 * d3));
    sm = 2.0 * d1;

    F << tl * tm,                   tl * tl * tl,           tm * tm * tm,           1,
        -(sm * tl) -(sl * tm),      -(3.0 * sl * tl * tl),  -(3.0 * sm * tm * tm),  0,
        sl * sm,                    3.0 * sl * sl * tl,     3.0 * sm * sm *  tm,    0,
        0,                          -(sl * sl * sl),        -(sm * sm * sm),        0;  

    break;

case CURVE_TYPE_LOOP:
    td = d2 + sqrt(4.0 * d1 * d3 - 3.0 * d2 *d2);
    sd = 2.0 * d1;
    te = d2 - sqrt(4.0 * d1 * d3 - 3.0 * d2 * d2);
    se = 2.0 * d1;


    F << td * te,                   td * td * te,                               td * te * te,                           1,      
        (-se * td) - (se * te), (-se * td * td) - (2.0 * sd * te * td),     (-sd * te * te) - (2.0 * se * td * te), 0,
        sd * se,                    te * sd * sd + 2.0 * se * td* sd,           td * se * se + 2 * sd * te * se,        0,
        0,                          -sd * sd * se,                              -sd * se * se,                          0;

    break;
...

I'm encountering a problem on serpentine type where the result appears wrong, From the paper, they say if(curveType == SERPENTINE && d1 < 0) i need to flip the k^3 - lm <= 0 equation to k^3 - lm > 0. It doesn't work. Strangely, the problem also occurs even when d1 > 0.