On Feb 12, 6:41 am, Steven D'Aprano <steve +comp.lang.pyt...@pearwood.info> 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! 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. 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... so it looks as though the values you're getting are at least alternating around the exact value. -- Mark -- http://mail.python.org/mailman/listinfo/python-list
