It doesn't look like roundoff to me. Nonexpert opinion.
hhr
On 2/4/2023 10:04 PM, 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
