# cubicUtils.py
# Andrew Davison, ad@coe.psu.ac.th, Nov. 2025

'''
https://en.wikipedia.org/wiki/Cubic_equation

'''


import math, cmath


def quadratic(z, a, b, c):
  return a*z**2 + b*z + c

def quadraticSolve(a, b, c):
  # Solve quadratic equation: ax^2 + bx + c = 0
  if a == 0:
    return linearSolve(b, c)
  disc = b**2 - 4*a*c
  sqrtD = cmath.sqrt(disc)
  x1 = (-b + sqrtD)/(2*a)
  x2 = (-b - sqrtD)/(2*a)
  return [x1, x2]


def linearSolve(a, b):
  # Solve linear equation: ax + b = 0
  if a == 0:
    if b == 0:
      return [math.nan]  # 0 = 0, infinite solns
    else:
      return []  # Inconsistent equ, no soln
  return [-b/a]



def cubic(z, a, b, c, d):
  return a*z**3 + b*z**2 + c*z + d


def toString(a, b, c, d):
  terms = []
  # x^3 term
  if a != 0:
    terms.append(f"{a}x^3" if a != 1 else "x^3")

  # x^2 term
  if b != 0:
    sign = " + " if b > 0 else " - "
    coef = abs(b)
    terms.append(f"{sign}{coef}x^2" if coef != 1 else f"{sign}x^2")

  # x term
  if c != 0:
    sign = " + " if c > 0 else " - "
    coef = abs(c)
    terms.append(f"{sign}{coef}x" if coef != 1 else f"{sign}x")

  # constant term
  if d != 0:
    sign = " + " if d > 0 else " - "
    terms.append(f"{sign}{abs(d)}")

  equ = "".join(terms) + " = 0"
  # Clean up leading plus sign if present
  if equ.startswith(" + "):
    equ = equ[3:]
  elif equ.startswith(" - "):
    equ = "-" + equ[3:]
  return equ



def divideCubicByLinear(cubicCoeffs, linearCoeffs):
  """Divide a cubic polynomial by a linear polynomial.
  cubicCoeffs: [a3, a2, a1, a0] for a3*x^3 + a2*x^2 + a1*x + a0
  linearCoeffs: [b1, b0] for b1*x + b0
  Returns (quotientCoeffs, remainder)
  e.g. divideCubicByLinear([1,6,11,6], [1,1]) -> [1,5,6],0
  """
  a3, a2, a1, a0 = cubicCoeffs
  b1, b0 = linearCoeffs
  
  # Initialize quotient coefficients for x^2, x, constant
  q2 = a3 / b1
  q1 = (a2 - q2 * b0) / b1
  q0 = (a1 - q1 * b0) / b1
  
  # Remainder is the leftover constant term
  remainder = a0 - q0 * b0
  
  return [q2, q1, q0], remainder


def linspace(start, stop, num):
  # Create evenly spaced numbers over a specified interval
  if num == 1:
    return [start]
  step = (stop - start) / (num - 1)
  return [start + step * i for i in range(num)]