I have no idea how robust the root finding algorithm in j602, but it
indeed sometimes gave wrong results in
I recommend using the gsl library or the lapack library, or spawning
octave which comes with a roots function that uses lapack. Gsl's root
finding is not recommended if there are roots with multiplicity or
very close to each other.
Here's an example that uses gsl.
NB. Example for finding roots of a polynomial with the gsl library
NB. Both functions take a list of _real_ coefficients, and return
NB. a list of roots. Apart from numerical errors, these should
NB. be functionally identical.
NB. First, the built-in J function.
rootsj =: /:~@:>@:{:@:p.
NB. Then the library call from gsl.
rootsgsl =: 3 :0
y =. ,y
s =. #y
pad =. {.'libgsl.so.0.9.0 gsl_poly_complex_workspace_alloc * x' 15!:0 ,<s
all =. 'libgsl.so.0.9.0 gsl_poly_complex_solve i *d x * *d' (15!:0)
y;s;pad;(0#~+:<:s)
'libgsl.so.0.9.0 gsl_poly_complex_workspace_free n *' (15!:0) ,<pad
if. 0~:>{.all do. %'root finding failed' end.
/:~_2+.^:_1\>4{all
)
NB. Example of usage.
p =: _1 1 1
echo rootsj p
echo rootsgsl p
NB. End
Ambrus
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
