def roots_cubic(a, b, c, d):
    r'''Cubic equation solver based on a variety of sources, algorithms, and 
    numerical tools. It seems evident after some work that no analytical
    solution using floating points will ever have high-precision results
    for all cases of inputs. Some other solvers, such as NumPy's roots
    which uses eigenvalues derived using some BLAS, seem to provide bang-on 
    answers for all values coefficients. However, they can be quite slow - and 
    where possible there is still a need for analytical solutions to obtain
    15-35x speed, such as when using PyPy.
    
    A particular focus of this routine is where a=1, b is a small number in the
    range -10 to 10 - a common occurence in cubic equations of state.

    Parameters
    ----------
    a : float
        Coefficient of x^3, [-]
    b : float
        Coefficient of x^2, [-]
    c : float
        Coefficient of x, [-]
    d : float
        Added coefficient, [-]

    Returns
    -------
    roots : tuple(float)
        The evaluated roots of the polynomial, 1 value when a and b are zero,
        two values when a is zero, and three otherwise, [-]

    Notes
    -----
    For maximum speed, provide Python floats. Compare the speed with numpy via:
        
    %timeit roots_cubic(1.0, 100.0, 1000.0, 10.0)
    %timeit np.roots([1.0, 100.0, 1000.0, 10.0])
    
    %timeit roots_cubic(1.0, 2.0, 3.0, 4.0)
    %timeit np.roots([1.0, 2.0, 3.0, 4.0])
    
    The speed is ~15-35 times faster; or using PyPy, 240-370 times faster.
    
    Examples
    --------
    >>> roots_cubic(1.0, 100.0, 1000.0, 10.0)
    (-0.010010019045111562, -88.73128838313305, -11.258701597821826)
    
    References
    ----------
    .. [1] "Solving Cubic Equations." Accessed January 5, 2019.
       http://www.1728.org/cubic2.htm.

    '''    
    '''
    Notes
    -----
    Known issue is inputs that look like
    1, -0.999999999978168, 1.698247818501352e-11, -8.47396642608142e-17
    Errors grown unbound, starting when b is -.99999 and close to 1.
    '''
    if b == 0.0 and a == 0.0: 
        return (-d/c, )
    elif a == 0.0:
        D = c*c - 4.0*b*d
        b_inv_2 = 0.5/b
        if D < 0.0:
            D = (-D)**0.5
            x1 = (-c + D*1.0j)*b_inv_2
            x2 = (-c - D*1.0j)*b_inv_2
        else:
            D = D**0.5
            x1 = (D - c)*b_inv_2
            x2 = -(c + D)*b_inv_2
        return (x1, x2)
    a_inv = 1.0/a
    a_inv2 = a_inv*a_inv
    bb = b*b
    '''Herbie modifications for f:
    c*a_inv - b_a*b_a*third
    '''
    
    b_a = b*a_inv
    b_a2 = b_a*b_a
    f = c*a_inv - b_a2*third
#    f = (3.0*c*a_inv - bb*a_inv2)*third 
    g = ((2.0*(bb*b) * a_inv2*a_inv) - (9.0*b*c)*(a_inv2) + (27.0*d*a_inv))*one_27
#    g = (((2.0/(a/b))/((a/b) * (a/b)) + d*27.0/a) - (9.0/a*b)*c/a)/27.0
    
    h = (0.25*(g*g) + (f*f*f)*one_27)
#    print(f, g, h)
    '''h has no savings on precicion - 0.4 error to 0.2.
    '''
#    print(f, g, h, 'f, g, h')
    if h == 0.0 and g == 0.0 and f == 0.0:
        if d/a >= 0.0:
            x = -((d*a_inv)**(third))
        else:
            x = (-d*a_inv)**(third)
        return (x, x, x)
    elif h > 0.0:
#        print('basic')
        # 1 real root, 2 imag
        root_h = h**0.5
        R = -(0.5*g) + root_h
        if R >= 0.0:
            S = R**third
        else:
            S = -((-R)**third)
        T = -(0.5*g) - root_h
        if T >= 0.0:
            U = (T**(third))
        else:
            U = -(((-T)**(third)))

        SU = S + U
        b_3a = b*(third*a_inv)
        t1 = -0.5*SU - b_3a
        t2 = (S - U)*complex_factor
        
        x1 = SU - b_3a
        # x1 is OK actually in some tests? the issue is x2, x3?
        x2 = t1 + t2
        x3 = t1 - t2
    elif h <= 0.0:
        t2 = a*a
        t3 = d*d
        t10 = c*c
        t14 = b*b
        t15 = t14*b
        
        '''This method is inaccurate when choice_term is too small; but still
        more accurate than the other method.
        '''
        choice_term = -18.0*a*b*c*d + 4.0*a*t10*c + 4.0*t15*d - t14*t10 + 27.0*t2*t3
        if abs(choice_term) > 1e-12 or abs(b + 1.0) < 1e-7:
#            print('mine')
            t32 = 1.0/a
            t20 = csqrt(choice_term)
            t31 = (36.0*c*b*a + 12.0*root_three*t20*a - 108.0*d*t2 - 8.0*t15)**third
            t33 = t31*t32
            t32_t31 = t32/t31
            
            x1 = (t33*sixth - two_thirds*(3.0*a*c - t14)*t32_t31 - b*t32*third).real
            t40 = (3.0*a*c - t14)*t32_t31
            
            t50 = -t33*twelfth + t40*third - b*t32*third
            t51 = 0.5j*root_three*(t33*sixth + two_thirds*t40)
            x2 = (t50 + t51).real
            x3 = (t50 - t51).real
        else:
#            print('other')
            # 3 real roots
            # example is going in here
            i = (((g*g)*0.25) - h)**0.5
            j = i**third # There was a saving for j but it was very weird with if statements!
            '''Clamied nothing saved for k.
            '''
            k = acos(-(g/(2.0*i)))
            L = -j
    
            # Would be nice to be able to compute the sin and cos at the same time
            M = cos(k*third)
            N = root_three*sin(k*third)
            P = -(b*(third*a_inv))
    
            # Direct formula for x1
            x1 = 2.0*j*M - b*(third*a_inv)
            x2 = L*(M + N) + P
            x3 = L*(M - N) + P
    return (x1, x2, x3)