(Fix pushed. Since CSZ1 has no factors in common with 8, i&7 should visit all
elements of cyclefracs; not in any particular order, but that shouldn't
matter.)
On Sat, 4 Feb 2023, Elijah Stone wrote:
I think the j implementation does exactly the same thing as the julia and bqn
implementations:
if(!--kicktimer){kicktimer=CSZ1;
x=zplus(x,ztymes(dx,zrj0(cyclefracs[(i>>3)&8])));
Actually, now that I look at it, that's wrong. I'm not sure what (i>>3)&8
was
supposed to be, but if I replace it with i&7, the code works fine.
On Sat, 4 Feb 2023, Henry Rich wrote:
That might well be what I'm seeing.
The J implementation tries Laguerre, and if it fails perturbs the
polynomial by 1e_12 in the leading coefficient. Do you think it would
be better if instead I retried with a random starting point?
I haven't checked for what the region of convergence of the starting
point is.
hhr
On 2/4/2023 9:55 PM, Marshall Lochbaum wrote:
Tried out my own Laguerre implementation from here:
https://github.com/mlochbaum/bqn-libs/blob/master/polynomial.bqn#L80-L113
I remember I used J as a reference for this, but I think I ended up
following this one in Julia:
https://github.com/giordano/PolynomialRoots.jl/blob/master/src/PolynomialRoots.jl#L255
With the starting point at 0 it converges, but only with a few random
adjustments from this line:
{ 0=jump_i|i ? root +↩ dx × Jump i÷jump_i ;@}
Here Jump i÷jump_i just ignores the argument and produces a random
float, ?0 in J (Julia has a table to avoid computing one).
Removing this, I see a small region of non-convergence around 0, radius
between 0.1 and 0.2. I suppose this isn't what you get?
Marshall
On Sat, Feb 04, 2023 at 08:34:46PM -0500, 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
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
