func SolveCubic(a, b, c, d float64) []float64 {
	// Handle degenerate case where a is zero
	aRecip := 1.0 / a
	const oneThird = 1.0 / 3.0

	scaledB := b * (oneThird * aRecip)
	scaledC := c * (oneThird * aRecip)
	scaledD := d * aRecip

	// Check if scaling resulted in non-finite values (a is zero or too small)
	if !isFinite(scaledB) || !isFinite(scaledC) || !isFinite(scaledD) {
		// Cubic coefficient is zero or nearly so - solve as quadratic
		return SolveQuadratic(b, c, d)
	}

	// Use scaled coefficients
	c0, c1, c2 := scaledD, scaledC, scaledB

	// (d0, d1, d2) is called "Delta" in the article
	d0 := (-c2)*c2 + c1
	d1 := (-c1)*c2 + c0
	d2 := c2*c0 - c1*c1

	// d is called "Discriminant"
	disc := 4.0*d0*d2 - d1*d1

	// de is called "Depressed.x", Depressed.y = d0
	de := (-2.0*c2)*d0 + d1

	if disc < 0.0 {
		// One real root
		sq := math.Sqrt(-0.25 * disc)
		r := -0.5 * de
		t1 := math.Cbrt(r+sq) + math.Cbrt(r-sq)
		return []float64{t1 - c2}
	} else if disc == 0.0 {
		// Two real roots (one is a double root)
		t1 := math.Copysign(math.Sqrt(-d0), de)
		return []float64{t1 - c2, -2.0*t1 - c2}
	}

	// Three distinct real roots
	th := math.Atan2(math.Sqrt(disc), -de) * oneThird
	thSin, thCos := math.Sincos(th)

	r0 := thCos
	ss3 := thSin * math.Sqrt(3.0)
	r1 := 0.5 * (-thCos + ss3)
	r2 := 0.5 * (-thCos - ss3)
	t := 2.0 * math.Sqrt(-d0)

	return []float64{
		t*r0 - c2,
		t*r1 - c2,
		t*r2 - c2,
	}
}

func SolveCubicInUnitInterval(a, b, c, d float64) []float64 {
	roots := SolveCubic(a, b, c, d)
	return filterRootsToUnitInterval(roots)
}

func SolveQuadratic(a, b, c float64) []float64 {
	// Scale coefficients to avoid overflow in discriminant calculation
	sc0 := c / a
	sc1 := b / a

	// Check if coefficients are valid (not Inf/NaN)
	if !isFinite(sc0) || !isFinite(sc1) {
		return solveQuadraticLinear(b, c)
	}

	// Normal case: valid quadratic
	return solveQuadraticNormal(sc0, sc1)
}

func SolveQuadraticInUnitInterval(a, b, c float64) []float64 {
	roots := SolveQuadratic(a, b, c)
	return filterRootsToUnitInterval(roots)
}