Regarding other ways to solve polynomials, I've found that the Aberth method is fast and accurate, and it can be implemented with array operations since it works on all roots at once. For example roots from a random 20-coefficient polynomial takes about 1ms. Bini's "Numerical Computation of Polynomial Zeros by Means of Aberth's Method" seems to have all the details needed for a really solid implementation, although I didn't finish working through it.
Marshall On Sat, Feb 04, 2023 at 07:04:06PM -0800, Elijah Stone wrote: > Is it possible the culprit is floating-point roundoff? If so, perhaps the > routine could be made to detect that automatically. You can get a > conservative estimation of floating-point roundoff by using interval > arithmetic, directing the hardware to alternately round towards positive and > negative infinity for each operation; when the estimated roundoff is larger > than the refinement you made to your root estimate, the latter is likely > junk, so you should fall back to some alternate strategy. > > (When changing the rounding mode is somewhat expensive, as on x86 before > avx512, you can get a slightly looser--but still pretty good--estimate by > multiplying by 1-2^-52 or 1+2^-52. Those 'round' towards and away from > zero, so some extra fanagling is required; this can be amortised by > partitioning all terms according to sign, which is probably a good idea > anyway.) > > It could also be worth looking at other approaches. Like I said before, I > think bisection using intervals is a likely approach. I also found this > thesis https://www.math.stonybrook.edu/~scott/Papers/thesis-alt.pdf which > talks about forcing newton to converge, which seems promising. > > On Sat, 4 Feb 2023, Henry Rich wrote: > > > Somehow I lost the original post. > > > > For some reason Laguerre's method doesn't converge for this polynomial. > > There seem to be large regions of non-convergence. Changes of up to 1e_5 > > in any of the constants still fail to converge. > > > > The point seems to oscillate in a wide range but doesn't head toward any > > solution. In the implementation the initial guess is always 0j0. > > > > I don't see why this would fail. If anyone on this list can help, or > > knows someone who can help, I'd appreciate suggestions. > > > > Henry Rich > > ---------------------------------------------------------------------- > > For information about J forums see http://www.jsoftware.com/forums.htm > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
