On Oct 2, 2015, at 2:16 PM, Louis Wasserman <wasserman.lo...@gmail.com> wrote:
> I'm pretty sure we could potentially public-domain our implementation, if > that were an issue. Personally, if it’s much better than mine (or what mine could be revised to be) I’d be happy to have the better outcome. > > The implementation I have here is effectively the one from Hacker’s Delight > > (2nd ed.). The initial guess is intended to be larger than the true result > > in order to simplify the termination condition of that algorithm as the > > monotonic convergence cannot have any oscillation between potential > > terminal values. > > This isn't a problem. The arithmetic mean of *any* two nonnegative numbers > is always greater than their geometric mean, so for *any* nonnegative a, (a + > x/a)/2 >= sqrt(a * x/a) = sqrt(x). On Oct 2, 2015, at 2:18 PM, Louis Wasserman <lowas...@google.com> wrote: > (https://en.wikipedia.org/wiki/Inequality_of_arithmetic_and_geometric_means > proves that the arithmetic mean is >= the geometric mean.) Got it. > So once you do *one* Newton iteration, the convergence is guaranteed to be > monotonically decreasing after that point. > > Newton's method doubles the number of correct digits with every iteration. > So you're paying one extra Newton iteration, but in exchange you're getting > -handwave- 50 correct bits to start out with. That *more* than pays for > itself. > > > There is certainly some room here for improvement in the range of input > > values less than Double.MAX_VALUE but this is a first stab at the > > implementation so hopefully such improvement may be brought in later if it > > is not in the first pass. > > It doesn't matter whether the input is bigger than Double.MAX_VALUE. You can > just find some even number, s, such that x.shiftRight(s) < 2^52. Then use > > doubleToBigInteger(Math.sqrt(x.shiftRight(s))).shiftLeft(s / 2) > > as your starting estimate. You're still getting ~50 correct bits to start > your Newton iteration. Excellent suggestion. I’ll look into revising it accordingly. Initially I had the thing broken into three ranges: 4 <= x <= Long.MAX_VALUE, Long.MAX_VALUE < x <= Double.MAX_VALUE, and Double.MAX_VALUE < x but found that this was lame and pointless. Thanks, Brian
