حل المعادلة التكعيبية في C++ وإيجاد الجذور الحقيقية بالتنصيف

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

std::vector<long double> cubicRealRoots(long double a, long double b,
                                        long double c, long double d) {
    for (auto value : {a, b, c, d})
        if (!std::isfinite(value)) throw std::invalid_argument("Use finite coefficients.");
    if (a == 0) throw std::invalid_argument("a must be nonzero for a cubic equation.");
    long double A = b / a, B = c / a, C = d / a;
    long double bound = 1 + std::max({std::abs(A), std::abs(B), std::abs(C)});
    if (!std::isfinite(bound) || bound > 1e6L)
        throw std::invalid_argument("Normalized coefficients are too large for this demo.");
    auto f = [=](long double x) { return ((x + A) * x + B) * x + C; };
    auto nearZero = [=](long double x) {
        long double scale = 1 + std::abs(x*x*x) + std::abs(A*x*x)
                              + std::abs(B*x) + std::abs(C);
        return std::abs(f(x)) <= 1e-14L * scale;
    };
    std::vector<long double> points{-bound};
    long double criticalDiscriminant = A * A - 3 * B;
    if (criticalDiscriminant >= 0) {
        long double root = std::sqrt(criticalDiscriminant);
        for (auto x : {(-A - root) / 3, (-A + root) / 3})
            if (x > -bound && x < bound) points.push_back(x);
    }
    points.push_back(bound);
    std::sort(points.begin(), points.end());
    std::vector<long double> roots;
    for (auto x : points) if (nearZero(x)) roots.push_back(x);
    for (std::size_t i = 1; i < points.size(); ++i) {
        long double left = points[i - 1], right = points[i];
        if (nearZero(left) || nearZero(right)) continue;
        if (std::signbit(f(left)) == std::signbit(f(right))) continue;
        for (int step = 0; step < 160; ++step) {
            long double middle = left + (right - left) / 2;
            long double value = f(middle);
            if (value == 0) { left = right = middle; break; }
            if (std::signbit(value) == std::signbit(f(left))) left = middle;
            else right = middle;
        }
        roots.push_back(left + (right - left) / 2);
    }
    std::sort(roots.begin(), roots.end());
    roots.erase(std::unique(roots.begin(), roots.end(), [](auto x, auto y) {
        return std::abs(x-y) <= 1e-7L * std::max({1.0L, std::abs(x), std::abs(y)});
    }), roots.end());
    return roots;
}

int main() {
    long double a{}, b{}, c{}, d{};
    std::cout << "Coefficients a b c d: ";
    if (!(std::cin >> a >> b >> c >> d)) {
        std::cout << "Invalid input.\n";
        return 1;
    }
    try {
        auto roots = cubicRealRoots(a, b, c, d);
        std::cout << std::setprecision(12) << "Distinct real roots:\n";
        for (auto x : roots) std::cout << x << '\n';
    } catch (const std::exception& error) {
        std::cout << error.what() << '\n';
        return 1;
    }
}