On Sun, 12 Feb 2012 12:18:15 -0800, Mark Dickinson wrote: > On Feb 12, 6:41 am, Steven D'Aprano <steve > [email protected]> wrote: > >> err = -a/b # Estimate of the error in the current w. >> if abs(err) <= 1e-16: >> break > > If the result you're expecting is around -1.005, this exit condition is > rather optimistic: the difference between the two Python floats either > side of this value is already 2.22e-16, so you're asking for less than > half a ulp of error!
I was gradually coming to the conclusion on my own that I was being overly optimistic with my error condition, although I couldn't put it into words *why*. Thanks for this Mark, this is exactly the sort of thing I need to learn -- as is obvious, I'm no expert on numeric programming. > As to the rest; your error estimate simply doesn't have enough > precision. The main problem is in the computation of a, where you're > subtracting two almost identical values. The absolute error incurred in > computing w*exp(w) is of the same order of magnitude as the difference > 'w*exp(w) - x' itself, so err has lost essentially all of its > significant bits, and is at best only a crude indicator of the size of > the error. The solution would be to compute the quantities 'exp(w), > w*exp(w), and w*exp(w) - x' all with extended precision. Other than using Decimal, there's no way to do that in pure Python, is there? We have floats (double) and that's it. > For the other > quantities, there shouldn't be any major issues---after all, you only > need a few significant bits of 'delta' for it to be useful, but with the > subtraction that generates a, you don't even get those few significant > bits. > >> (The correct value for w is approximately -1.00572223991.) > > Are you sure? Wolfram Alpha gives me the following value for W(-1, > -0.36787344117144249455719773322925902903079986572265625): > > -1.005722239691522978... I did say *approximately*. The figure I quote comes from my HP-48GX, and seems to be accurate to the precision offered by the HP. -- Steven -- http://mail.python.org/mailman/listinfo/python-list
