Classifying cubic bezier curves according to Loop & Blinn 2005

Question

I am trying flatten a bezier path (remove all intersections and replace with end points) as part of the implementation of the rendering algorithm described here and found the algorithm described by Loop and Blinn in GPU Gems 3 Ch. 25, with correction for $a3$'s cross product to detect curve self-intersection (loops).

The algorithm calls for the evaluation of the expression discr=d₁²(3d₂²-4d₁d₃) such that the curve is self-intersecting iff the discr < 0.

where

d₁=a₁-2a₂+3a₃, d₂=-a₂+3a₃, d₃=3a₃

a₁=b₀·(b₃⨯b₂), a₂=b₁·(b₀⨯b₃), a₃=b₂·(b₁⨯b₀)

I have implemented the algorithm in the sample code below, as well as that of a few other algorithms I have been able to find. None of them agree, and all of them are wrong. Specifically, I have implemented

1. The algorithm described in GPU Gems 3
2. The same algorithm as described in the original paper.
3. The same algorithm as implemented in paper.js (issue)(sketch).
4. The curve canonicalization algorithm described in Pomax' Bezier Primer.

That all three algorithms disagree suggests a major misunderstanding. What is it, and (more importantly) how might I identify such an issue myself?

I referred to the questions here and here, as well as the paper in which Loop & Blinn describe their approach, to no effect.

#[derive(Clone, Copy, Debug, PartialEq)]
enum Kind {
    Serpentine,
    Cusp,
    Loop,
    Quadratic,
    Line,
}

#[derive(Clone, Copy)]
pub struct Point {
    x: f32,
    y: f32,
}

#[derive(Clone, Copy)]
struct Vec3 {
    x: f32,
    y: f32,
    z: f32,
}

impl Vec3 {
    fn from_point(p: Point) -> Vec3 {
        Vec3 {
            x: p.x,
            y: p.y,
            z: 1.0,
        }
    }

    fn dot(&self, other: Vec3) -> f32 {
        self.x * other.x + self.y * other.y + self.z * other.z
    }

    fn cross(&self, other: Vec3) -> Vec3 {
        Vec3 {
            x: self.y * other.z - self.z * other.y,
            y: self.z * other.x - self.x * other.z,
            z: self.x * other.y - self.y * other.x,
        }
    }
}

fn is_zero(f: f32) -> bool {
    f.abs() < 0.000001
}

fn approx_eq(a: f32, b: f32) -> bool {
    is_zero((a - b).abs())
}

fn normalize(a: f32, b: f32, c: f32) -> (f32, f32, f32) {
    let len = (a * a + b * b + c * c).sqrt();
    if is_zero(len) {
        (0.0, 0.0, 0.0)
    } else {
        (a / len, b / len, c / len)
    }
}

pub fn classify(curve: &[Point; 4]) {
    let p0 = Vec3::from_point(curve[0]);
    let p1 = Vec3::from_point(curve[1]);
    let p2 = Vec3::from_point(curve[2]);
    let p3 = Vec3::from_point(curve[3]);

    let loop_blinn = {
        let det1 = -p3.dot(p2.cross(p0));
        let det2 = p3.dot(p2.cross(p0));
        let det3 = -p2.dot(p1.cross(p0));

        let (det1, det2, det3) = normalize(det1, det2, det3);

        let discr = det2 * det2 * (3.0 * det2 * det2 - 4.0 * det3 * det1);

        if is_zero(det1) {
            if !is_zero(det2) {
                Kind::Cusp
            } else if is_zero(det3) {
                Kind::Line
            } else {
                Kind::Quadratic
            }
        } else if is_zero(discr) {
            Kind::Cusp
        } else if discr > 0.0 {
            Kind::Serpentine
        } else {
            Kind::Loop
        }
    };

    let gpu_gems = {
        let a1 = p0.dot(p3.cross(p2));
        let a2 = p1.dot(p0.cross(p3));
        let a3 = p2.dot(p1.cross(p0));
        let d1 = a1 - 2.0 * a2 + 3.0 * a3;
        let d2 = -a2 + 3.0 * a3;
        let d3 = 3.0 * a3;
        let discr = d2 * d2 * (3.0 * d2 * d2 - 4.0 * d1 * d3);

        if is_zero(d1) {
            if is_zero(d2) && is_zero(d3) {
                Kind::Line
            } else {
                Kind::Quadratic
            }
        } else if is_zero(discr) {
            Kind::Cusp
        } else if discr < 0.0 {
            Kind::Serpentine
        } else {
            Kind::Loop
        }
    };

    let paper_js = {
        let x1 = p0.x;
        let y1 = p0.y;
        let x2 = p1.x;
        let y2 = p1.y;
        let x3 = p2.x;
        let y3 = p2.y;
        let x4 = p3.x;
        let y4 = p3.y;

        let a1 = x1 * (y4 - y3) + y1 * (x3 - x4) + x4 * y3 - x3 * y4;
        let a2 = x2 * (y1 - y4) + y2 * (x4 - x1) + x1 * y4 - y1 * x4;
        let a3 = x3 * (y2 - y1) + y3 * (x1 - x2) + x2 * y1 - y2 * x1;
        let d3 = 3.0 * a3;
        let d2 = d3 - a2;
        let d1 = d2 - a2 + a1;
        let (d1, d2, d3) = normalize(d1, d2, d3);

        let d = 3.0 * d2 * d2 - 4.0 * d1 * d3;

        if is_zero(d1) {
            if is_zero(d2) {
                if is_zero(d3) {
                    Kind::Line
                } else {
                    Kind::Quadratic
                }
            } else {
                Kind::Serpentine
            }
        } else if is_zero(d) {
            Kind::Cusp
        } else if d > 0.0 {
            Kind::Serpentine
        } else {
            Kind::Loop
        }
    };

    let pomax_primer = {
        let y31 = p3.y / p1.y;
        let y21 = p2.y / p1.y;
        let x32 = (p3.x - p2.x * y31) / (p2.x - p1.x * y21);

        let x = x32;
        let y = y31 + x32 * (1.0 - y21);

        let cusp_line = (-(x * x) + 2.0 * x + 3.0) / 4.0;
        let loop_at_1 = ((3.0 * 4.0 * x - x * x).sqrt() - x) / 2.0;
        let loop_at_0 = (-(x * x) + 3.0 * x) / 3.0;

        if x > 1.0 || y > 1.0 {
            Kind::Serpentine
        } else if (0.0..1.0).contains(&x) {
            if approx_eq(loop_at_1, y) {
                Kind::Loop
            } else if approx_eq(cusp_line, y) {
                Kind::Cusp
            } else if y < cusp_line || y > loop_at_1 {
                Kind::Serpentine
            } else {
                Kind::Loop
            }
        } else if approx_eq(loop_at_0, y) {
            Kind::Loop
        } else if approx_eq(cusp_line, y) {
            Kind::Cusp
        } else if y < cusp_line || y > loop_at_0 {
            Kind::Loop
        } else {
            Kind::Serpentine
        }
    };

    println!("\tMethod 1 (Loop Blinn 2005): {:?}", loop_blinn);
    println!("\tMethod 2 (GPU Gems 3.25):   {:?}", gpu_gems);
    println!("\tMethod 3 (Paper.js):        {:?}", paper_js);
    println!("\tMethod 4 (Pomax Primer):    {:?}", pomax_primer);
}

pub fn main() {
    let point = [
        Point { x: 1.0, y: 1.0 },
        Point { x: 1.0, y: 1.0 },
        Point { x: 1.0, y: 1.0 },
        Point { x: 1.0, y: 1.0 },
    ];
    println!("Expecting: Kind::Line");
    classify(&point);

    let line = [
        Point { x: 1.0, y: 1.0 },
        Point { x: 2.0, y: 1.0 },
        Point { x: 3.0, y: 1.0 },
        Point { x: 4.0, y: 1.0 },
    ];
    println!("Expecting: Kind::Line");
    classify(&line);

    // loop
    let normal_loop = [
        Point { x: 75.0, y: 98.0 },
        Point { x: 195.0, y: 201.0 },
        Point { x: 63.0, y: 198.0 },
        Point { x: 135.0, y: 103.0 },
    ];
    println!("Expecting: Kind::Loop");
    classify(&normal_loop);

    let end_loop = [
        Point { x: 120.0, y: 331.0 },
        Point { x: 246.0, y: 261.0 },
        Point { x: 187.0, y: 242.0 },
        Point { x: 148.0, y: 314.0 },
    ];
    println!("Expecting: Kind::Loop");
    classify(&end_loop);

    // cusp
    let cusp = [
        Point { x: 272.0, y: 305.0 },
        Point { x: 93.0, y: 223.0 },
        Point { x: 221.0, y: 223.0 },
        Point { x: 148.0, y: 314.0 },
    ];
    println!("Expecting: Kind::Cusp");
    classify(&cusp);

    // serpentine
    let serpentine = [
        Point { x: 148.0, y: 314.0 },
        Point { x: 187.0, y: 242.0 },
        Point { x: 246.0, y: 261.0 },
        Point { x: 272.0, y: 305.0 },
    ];
    println!("Expecting: Kind::Serpentine");
    classify(&serpentine);

    let serpentine2 = [
        Point { x: 110.0, y: 150.0 },
        Point { x: 25.0, y: 190.0 },
        Point { x: 210.0, y: 250.0 },
        Point { x: 210.0, y: 30.0 },
    ];
    println!("Expecting: Kind::Serpentine");
    classify(&serpentine2);
}