Cubic Equation Solver

            from math import acos, cos, cbrt, pi, sqrt

def cubic_equation(a: float, b: float, c: float, d: float) -> float:
    if a == 0:
        raise ValueError("a cannot be zero")
    A, B, C = b / a, c / a, d / a
    p = B - A * A / 3
    q = 2 * A**3 / 27 - A * B / 3 + C
    discriminant = q * q / 4 + p**3 / 27
    if discriminant >= 0:
        return cbrt(-q / 2 + sqrt(discriminant)) + cbrt(-q / 2 - sqrt(discriminant)) - A / 3
    radius = 2 * sqrt(-p / 3)
    theta = acos((3 * q / (2 * p)) * sqrt(-3 / p))
    return min(radius * cos((theta + 2 * pi * k) / 3) - A / 3 for k in range(3))

assert abs(cubic_equation(1, -6, 11, -6) - 1) < 1e-9
          

            #include <assert.h>
#include <math.h>

double cubic_equation(double a, double b, double c, double d) {
    double A=b/a, B=c/a, C=d/a;
    double p=B-A*A/3.0, q=2*A*A*A/27.0-A*B/3.0+C;
    double disc=q*q/4.0+p*p*p/27.0;
    if (disc >= 0.0)
        return cbrt(-q/2+sqrt(disc))+cbrt(-q/2-sqrt(disc))-A/3;
    double r=2*sqrt(-p/3), theta=acos((3*q/(2*p))*sqrt(-3/p));
    double best=INFINITY;
    for(int k=0;k<3;k++) best=fmin(best,r*cos((theta+2*M_PI*k)/3)-A/3);
    return best;
}

int main(void) { assert(fabs(cubic_equation(1,-6,11,-6)-1)<1e-9); }
          

            #include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <numbers>

double cubic_equation(double a, double b, double c, double d) {
    const double A=b/a, B=c/a, C=d/a;
    const double p=B-A*A/3.0, q=2*A*A*A/27.0-A*B/3.0+C;
    const double disc=q*q/4.0+p*p*p/27.0;
    if (disc >= 0.0)
        return std::cbrt(-q/2+std::sqrt(disc))+std::cbrt(-q/2-std::sqrt(disc))-A/3;
    const double r=2*std::sqrt(-p/3), theta=std::acos((3*q/(2*p))*std::sqrt(-3/p));
    std::array<double,3> roots{};
    for(int k=0;k<3;k++) roots[k]=r*std::cos((theta+2*std::numbers::pi*k)/3)-A/3;
    return *std::min_element(roots.begin(), roots.end());
}

int main() { assert(std::abs(cubic_equation(1,-6,11,-6)-1)<1e-9); }
          

x86-64 NASM · System V ABI · Linux · double arguments and result in XMM registers

            ; double cubic_equation(double a, double b, double c, double d)
; Newton iteration returns one real root. xmm0=a, xmm1=b, xmm2=c, xmm3=d.
global cubic_equation
section .text
cubic_equation:
    movapd xmm7, xmm0       ; a
    movapd xmm6, xmm1       ; b
    movapd xmm5, xmm2       ; c
    movapd xmm4, xmm3       ; d
    pxor xmm0, xmm0         ; starting guess x=0
    mov ecx, 40
.iterate:
    movapd xmm1, xmm0
    mulsd xmm1, xmm0
    movapd xmm2, xmm1
    mulsd xmm2, xmm0
    mulsd xmm2, xmm7
    mulsd xmm1, xmm6
    addsd xmm2, xmm1
    movapd xmm1, xmm0
    mulsd xmm1, xmm5
    addsd xmm2, xmm1
    addsd xmm2, xmm4        ; f(x)
    movapd xmm1, xmm0
    mulsd xmm1, xmm0
    addsd xmm1, xmm1
    movapd xmm3, xmm0
    mulsd xmm3, xmm0
    addsd xmm1, xmm3
    mulsd xmm1, xmm7        ; 3ax²
    movapd xmm3, xmm0
    addsd xmm3, xmm3
    mulsd xmm3, xmm6
    addsd xmm1, xmm3
    addsd xmm1, xmm5        ; f'(x)
    divsd xmm2, xmm1
    subsd xmm0, xmm2
    loop .iterate
    ret
          

            function result = cubic_equation(a, b, c, x)
    values = roots([a, b, c, x]);
    realValues = real(values(abs(imag(values)) < 1e-9));
    result = min(realValues);
end
          

            ClearAll[mwCalculate];
mwCalculate[a_, b_, c_, d_] := Min[x /. NSolve[a*x^3 + b*x^2 + c*x + d == 0, x, Reals]];