METODO DE FERRARI EN R Y PYTHON

resolver_cuartica <- function(a, b, c, d, e) {
  if (a == 0) stop("El coeficiente 'a' debe ser distinto de cero")

  # Normalizar los coeficientes
  b <- b / a
  c <- c / a
  d <- d / a
  e <- e / a
  
  # Resolver la ecuación cúbica auxiliar
  cubic <- function(p, q, r) {
    Q <- (3 * q - p^2) / 9
    R <- (9 * p * q - 27 * r - 2 * p^3) / 54
    D <- Q^3 + R^2
    
    if (D >= 0) {
      S <- sign(R + sqrt(D)) * abs(R + sqrt(D))^(1/3)
      T <- sign(R - sqrt(D)) * abs(R - sqrt(D))^(1/3)
      roots <- c(S + T - p / 3)
    } else {
      theta <- acos(R / sqrt(-Q^3))
      roots <- c(2 * sqrt(-Q) * cos(theta / 3) - p / 3,
                 2 * sqrt(-Q) * cos((theta + 2 * pi) / 3) - p / 3,
                 2 * sqrt(-Q) * cos((theta + 4 * pi) / 3) - p / 3)
    }
    return(roots)
  }
  
  p <- -b^2 / 12 - e
  q <- -b^3 / 108 + b * e / 3 - d^2 / 8
  r <- complex(real = -q / 2, imaginary = sqrt(abs(q^2 / 4 + p^3 / 27)))
  u <- sign(Re(r)) * abs(r)^(1/3)
  
  y <- u - p / (3 * u)
  
  w <- sqrt(complex(real = b^2 / 4 - 2 * Re(y), imaginary = sqrt(abs((b^2 / 4 - Re(y))^2 - e))))
  z <- sqrt(complex(real = b^2 / 4 - 2 * Re(y), imaginary = -sqrt(abs((b^2 / 4 - Re(y))^2 - e))))
  
  x1 <- -b / 4 + (w + z) / 2
  x2 <- -b / 4 + (w - z) / 2
  x3 <- -b / 4 - (w + z) / 2
  x4 <- -b / 4 - (w - z) / 2
  
  return(c(x1, x2, x3, x4))
}

# Ejemplo de uso
a <- 1
b <- -8
c <- 14
d <- -8
e <- 1
raices <- resolver_cuartica(a, b, c, d, e)
print(raices)
[1]  5.657097+0.000000i  2.000000+1.784999i -1.657097+0.000000i
[4]  2.000000-1.784999i
import cmath

def resolver_cuartica(a, b, c, d, e):
    if a == 0:
        raise ValueError("El coeficiente 'a' debe ser distinto de cero")
    
    # Normalizar los coeficientes
    b /= a
    c /= a
    d /= a
    e /= a
    
    # Resolver la ecuación cúbica auxiliar
    def cubic(p, q, r):
        Q = (3 * q - p**2) / 9
        R = (9 * p * q - 27 * r - 2 * p**3) / 54
        D = Q**3 + R**2
        
        if D >= 0:
            S = cmath.exp(cmath.log(R + cmath.sqrt(D)) / 3)
            T = cmath.exp(cmath.log(R - cmath.sqrt(D)) / 3)
            roots = [S + T - p / 3]
        else:
            theta = cmath.acos(R / cmath.sqrt(-Q**3))
            roots = [
                2 * cmath.sqrt(-Q) * cmath.cos(theta / 3) - p / 3,
                2 * cmath.sqrt(-Q) * cmath.cos((theta + 2 * cmath.pi) / 3) - p / 3,
                2 * cmath.sqrt(-Q) * cmath.cos((theta + 4 * cmath.pi) / 3) - p / 3
            ]
        return roots
    
    p = -b**2 / 12 - e
    q = -b**3 / 108 + b * e / 3 - d**2 / 8
    r = cmath.sqrt(q**2 / 4 + p**3 / 27)
    u = cmath.exp(cmath.log(-q / 2 + r) / 3)
    
    y = u - p / (3 * u)
    
    w = cmath.sqrt(b**2 / 4 - 2 * y + cmath.sqrt((b**2 / 4 - y)**2 - e))
    z = cmath.sqrt(b**2 / 4 - 2 * y - cmath.sqrt((b**2 / 4 - y)**2 - e))
    
    x1 = -b / 4 + (w + z) / 2
    x2 = -b / 4 + (w - z) / 2
    x3 = -b / 4 - (w + z) / 2
    x4 = -b / 4 - (w - z) / 2
    
    return [x1, x2, x3, x4]

# Ejemplo de uso
a = 1
b = -8
c = 14
d = -8
e = 1
raices = resolver_cuartica(a, b, c, d, e)
print(raices)
[(4.412319051064921+0.8450854100461852j), (4.412319051064921-0.8450854100461852j), (-0.4123190510649217-0.8450854100461852j), (-0.4123190510649217+0.8450854100461852j)]