quintic.007

-- SPDX-License-Identifier: MPL-2.0
-- CONFIDENTIAL — EXPERIMENTAL — NOT FOR DISTRIBUTION
--
-- Solvable Quintic Equations
--
-- Bring-Jerrard form: x⁵ + x + a = 0
-- Solved numerically via Newton's method:
--   x_{n+1} = x_n - f(x_n) / f'(x_n)
--
-- where f(x) = x⁵ + x + a, f'(x) = 5x⁴ + 1

-- The quintic polynomial
@pure fn quintic(x: Float, a: Float) -> Float {
    return pow(x, 5.0) + x + a
}

-- Its derivative
@pure fn quintic_deriv(x: Float) -> Float {
    return 5.0 * pow(x, 4.0) + 1.0
}

-- One Newton step
@pure fn newton_step(x: Float, a: Float) -> Float {
    let fx = quintic(x, a)
    let fpx = quintic_deriv(x)
    return x - fx / fpx
}

-- Solve x⁵ + x + 1 = 0 (a = 1)
-- Starting from x₀ = -0.5
@total data a_val = 1.0
@total data x0 = neg(0.5)
@total data x1 = newton_step(x0, a_val)
@total data x2 = newton_step(x1, a_val)
@total data x3 = newton_step(x2, a_val)
@total data x4 = newton_step(x3, a_val)
@total data x5 = newton_step(x4, a_val)
@total data x6 = newton_step(x5, a_val)

-- Verify: f(x6) should be ≈ 0
@total data residual = quintic(x6, a_val)

-- Solve x⁵ + x - 2 = 0 (a = -2, obvious root x = 1)
@total data b_val = neg(2.0)
@total data y0 = 0.5
@total data y1 = newton_step(y0, b_val)
@total data y2 = newton_step(y1, b_val)
@total data y3 = newton_step(y2, b_val)
@total data y4 = newton_step(y3, b_val)
@total data residual_b = quintic(y4, b_val)