mersenneforum.org integer nth root code
 Register FAQ Search Today's Posts Mark Forums Read

 2009-09-28, 20:51 #1 bsquared     "Ben" Feb 2007 5·727 Posts integer nth root code 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: $x_{k+1} = \frac{1}{n}\left[{x_{k} (n-1) + \frac{A}{x_k^{n-1}}}\right]$ 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 $x_k = 3$, 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?
 2009-09-29, 04:46 #2 wblipp     "William" May 2003 New Haven 22·593 Posts 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 Last fiddled with by wblipp on 2009-09-29 at 04:47
2009-09-29, 05:06   #3
paleseptember

Jun 2008
Wollongong, .au

3×61 Posts

Quote:
 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.

2009-09-29, 12:46   #4
bsquared

"Ben"
Feb 2007

5×727 Posts

Quote:
 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.

Quote:
 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.

Last fiddled with by bsquared on 2009-09-29 at 12:46 Reason: value to range...

 2009-09-29, 13:52 #5 jasonp Tribal Bullet     Oct 2004 DD916 Posts 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.
 2009-09-29, 19:39 #6 Gammatester     Mar 2009 2·19 Posts 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
2009-09-29, 19:54   #7
bsquared

"Ben"
Feb 2007

5×727 Posts

Quote:
 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.

Quote:
 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!

 Similar Threads Thread Thread Starter Forum Replies Last Post jasong Miscellaneous Math 5 2016-04-24 03:40 Damian Math 3 2010-01-01 01:56 jasonp Software 17 2009-01-29 02:25 mfgoode Homework Help 10 2007-10-05 04:12 junky NFSNET Discussion 0 2004-03-26 02:14

All times are UTC. The time now is 21:27.

Mon Aug 15 21:27:55 UTC 2022 up 39 days, 16:15, 1 user, load averages: 1.21, 1.37, 1.37