Ruby univariate polynomial root finder

univariate_polynomial_roots.rb

require 'matrix'

# Input: an array of the n coefficients [a_n, a_n-1,..., a_1] of 
# a univariate polynomial p of degree n ([a_n]x^n + [a_n-1]x^n-1 + ... + a_1)
#
# Output: an array of all n roots of p
#
# Exploits the fact that the eigenvalues of the companion matrix of the
# monic equivalent of p are the roots of p
#
# Currently not particularly accurate, though pretty cool...

def polynomial_roots(coeffs)
  coeffs = coeffs.map(&:to_f)
  size = coeffs.size - 1
  a = -coeffs.shift
  coeffs.map! { |c| c / a }.reverse!
  m = Matrix.build(size, size - 1) { |row, col| row - col == 1 ? 1 : 0 }
  companion_matrix = Matrix.columns(m.to_a.transpose << coeffs)
  companion_matrix.eigen.eigenvalues
end

# Examples

# solve the fifth degree Hermite polynomial

p polynomial_roots([32, 0, -160, 0, 120, 0])

# => [0.0, 0.9585724646138186, -0.9585724646138186, 2.0201828704560856, -2.0201828704560856]

# solve x^4-7x^3+17x^2-17x+6 = 0

p polynomial_roots([1, -7, 17, -17, 6])

# => [0.9999999715070083, 1.000000028492992, 1.999999999999998, 3.0000000000000018]
# actual: [1, 1, 2, 3]

# compute the 7th roots of unity

p polynomial_roots([1, 0, 0, 0, 0, 0, 0, -1])

# => [(-0.9009688679024196+0.43388373911755823i), (-0.9009688679024196-0.43388373911755823i), (-0.22252093395631456+0.9749279121818226i), (-0.22252093395631456-0.9749279121818226i), 0.9999999999999999, (0.6234898018587336+0.7818314824680304i), (0.6234898018587336-0.7818314824680304i)]