Numerically Stable Method for Solving Quadratic Equations

from numpy import allclose
from cmath import sqrt

def solve_quad(b, c):
    """Решает квадратное уравнение, x**2 + bx + c = 0.

    Parameteres
    ----------
    b, c : float
       Коэффициенты

    Returns
    -------
    x1, x2 : float or complex
       Корни.
    """
    if b >= 0:
        x1 = (-b - sqrt(b ** 2 - 4 * c)) / 2
        x2 = 2 * c / (-b - sqrt(b ** 2 - 4 * c))
    else:
        x1 = 2 * c / (-b + sqrt(b ** 2 - 4 * c))
        x2 = (- b + sqrt(b ** 2 - 4 * c)) / 2
    return x1, x2

variants = [{'b': 4.0, 'c': 3.0},
            {'b': 2.0, 'c': 1.0},
            {'b': 0.5, 'c': 4.0},
            {'b': 1e10, 'c': 3.0},
            {'b': -1e10, 'c': 4.0},]

for var in variants:
    x1, x2 = solve_quad(**var)
    print(allclose(x1*x2, var['c']))