NB. Find roots of polynomial by eigenvalue method.
NB. Let polynoomial be a0 + a1*x + a2*x^2 + ... + a_m*x^m
NB. Usage: poly_root a0 a1 ... a_m ==> x_1, x_2, ... x_m
NB. (This is just to show the use of LAPACK -- J has a
NB. built-in primitive, p., that gives the roots directly.)
require jpath '~addons\math\lapack\lapack.ijs'
require jpath '~addons\math\lapack\dgeev.ijs'
dgeev=: dgeev_jlapack_
NB. ---------------------------------------------
poly_root=: monad define
m=. <:n=. #coef=. y
NB. normalize coefficients by the highest:
norm=. |. }: (coef % {: coef)
NB. Form matrix: see Num. Recipes p 368.
A=. 1 |."1(1* =@i. m)
A=. (-norm) 0}A
NB. Use LAPACK to get eigenvalues, and grade up:
/:~ >1{ dgeev A
)