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
