integer nth root code

[bsquared]

I'm going back to look at my integer nth root code per this, and I'm running into a problem.

I'm trying variations like so:
nroot(3^180,150)

The actual integer root is 3, and my code's first guess is 4. So far so good. I then compute the next term in the recurrence:



and since 4^149 > 3^180, the fraction in the bracket is zero and the next term in the sequence is 3. Great! But we can't stop there because this new term is not the same as the old one (hasn't converged).

So now with , we get a huge term for the fraction in the bracket, the next term in the sequence is now something like 4 trillion, and the sequence never converges.

What do I do in situations like this, when even a correct guess of the integer nth root produces a wildly inaccurate next guess in the recurrence?

[wblipp]

This isn't modular arithemetic - you are looking for the rounded off (or perhaps down) value of the real root - right?

Section 9 of Numerical Recipes says you should always bracket a root and never let your iteration method get outside of the best known bounds at any stage. I've had very good experience with the Van Wijngaarden-Dekker-Brent method of Section 9.3 in a wide variety of numerical root finding problems.

If I understand the problem correctly, from a bracket of (3,4) you should move into an endgame of picking 3 or 4, not iterating further.

William

[paleseptember]

Originally Posted by wblipp

This isn't modular arithemetic - you are looking for the rounded off (or perhaps down) value of the real root - right?

Yes, the integer root is the floor function of the nth root.

Sorry, my programming skills are minimal, so I can't actually assist with the mechanics of the computation.

[bsquared]

Originally Posted by wblipp

This isn't modular arithemetic - you are looking for the rounded off (or perhaps down) value of the real root - right?

Yes, the floor of the real root, as paleseptember says.

Originally Posted by wblipp

Section 9 of Numerical Recipes says you should always bracket a root and never let your iteration method get outside of the best known bounds at any stage. I've had very good experience with the Van Wijngaarden-Dekker-Brent method of Section 9.3 in a wide variety of numerical root finding problems.

If I understand the problem correctly, from a bracket of (3,4) you should move into an endgame of picking 3 or 4, not iterating further.

William

I'll get out my NR and take a look, thanks. Basically, this seems to jive with what I was thinking... I know that the answer should be a certain number of bits (bits(input) / n, where n is the root, to be exact), so if the iteration goes outside that range after being inside that range, stop iterating and decide between the last couple values. I've already got code in place to detect oscillations by keeping a history of the last few iterations, so this shouldn't be hard to implement.

[jasonp]

In the specific case of needing the closest integer to the root, Crandall and Pomerance show that you can stop the iteration when the sequence of root estimates stops decreasing. That leads to a skewed loop where the convergence check is in the middle, after A/x_k^(n-1) has been computed, and x_(k+1) only gets computed if the convergence test fails.

[Gammatester]

As Jason said, you can stop the iteration if the sequence stops decreasing. But you have to choose a starting value greater than the real root. Normalyy you start with x0 = 2^ceil(bitsize(A)/n).

If n is large and the guesses are not very close to the root, convergence is only linear not quadratic. You should check the first few x_i for progress, otherwise it is better to use bisection steps until the iterates come close to the final result.

Wolfgang

[bsquared]

Originally Posted by Gammatester

As Jason said, you can stop the iteration if the sequence stops decreasing. But you have to choose a starting value greater than the real root. Normalyy you start with x0 = 2^ceil(bitsize(A)/n).

Yep. This is implemented now and seems to work fine.

Originally Posted by Gammatester

If n is large and the guesses are not very close to the root, convergence is only linear not quadratic. You should check the first few x_i for progress, otherwise it is better to use bisection steps until the iterates come close to the final result.

Wolfgang

The Van Wijngaarden-Dekker-Brent method that wblipp suggested includes allowances for that, I think. Implementing that method in YAFU is an exercise for another day. Stuff like this is very interesting and a fun diversion but I don't believe many people use YAFU for everyday bigint calculations for which they need the fastest and most bulletproof methods. I'm still concentrating on factoring improvements for now. That said, by all means continue to find and report these things as they come up, it's very much appreciated!