docs: 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)