require 'mathn'
module Math
  def self.roots_of_cubic_function(a,b,c,d,with_complex = false)
    f = ((3*c/a) - (b**2/a**2))/3
    g = (((2*b**3)/a**3) - (9*b*c/a**2) + (27*d/a))/27
    h = (g**2/4) + (f**3/27)
    if f == 0 && g == 0 && h == 0
      f = (3*c/a) - (b**2/a**2)
      g = (2*b**3/a**3) - (9*b*c/a**2) + (27*d/a)
      h = (g**2/4) + (f**3/27)
      [-Math.cbrt(d/a)]
    elsif h <= 0
      i = Math.sqrt((g**2/4) - h)
      j = Math.cbrt(i)
      k = Math.acos(-(g / (2*i)))
      l = j * -1
      m = Math.cos(k/3)
      n = Math.sqrt(3) * Math.sin(k/3)
      p = (b/(3*a))*-1
      [2*j*Math.cos(k/3) - (b/(3*a)), l*(m+n)+p, l*(m-n)+p]
    else
      f = (((3*c)/a) - (b**2/a**2))/3
      g = ((2*b**3/a**3) - (9*b*c/a**2) + ((27*d)/a))/27
      h = (g**2/4) + (f**3/27)
      r = -(g/2) + Math.sqrt(h)
      s = Math.cbrt(r)
      t = -(g/2) - Math.sqrt(h)
      u = Math.cbrt(t)
      roots = [(s + u) - (b/(3*a))]
      roots |= [1,-1].map do |algebraic_sign|
        Complex((-(s + u)/2) - (b/(3*a)),algebraic_sign * (s-u)*Math.sqrt(3)/2)
      end if with_complex
      roots
    end.map { |n| n.respond_to?(:round) ? n.round(10) : n }.to_set
  end

  def self.roots_of_quadratic_function(a,b,c)
    p = b/a
    q = c/a
    denom = (p/2.0)**2 - q
    return Set[] if denom < 0
    root = Math.sqrt(denom)
    Set.new([:+,:-].map{ |operator| (-(p/2.0)).send(operator, root)})
  end

  def self.roots_of_quartic_function(*args)
    a,b,c,d,e = args.map { |n| n/args.first }
    f = c - (3*b**2/8)
    g = d + (b**3/8) - (b*c/2)
    h = e - (3*b**4/256) + (b**2 * c/16) - ( b*d/4)
    roots = self.roots_of_cubic_function(1,(f/2),(((f**2)-(4*h))/16),-(g**2)/64,true)
    combs = roots.reject{ |n| n == 0}.combination(2)
    only_complex = combs.select { |comb| comb.all? { |n| n.is_a? Complex }}
    p,q = (only_complex.empty? ? combs.first : only_complex.first).map { |n| Math.sqrt(n) }
    r = (-g/(8*p*q))
    s = b/(4*a)
    roots = Set.new([p + q + r -s, p - q - r -s, -p + q - r -s, -p - q + r -s].map { |n| n.is_a?(Complex) ? (n.real if n.imaginary == 0) : n }.compact )
    roots.map { |n| n.respond_to?(:round) ? n.round(10) : n }.to_set
  end
end