[issue27761] Private _nth_root function loses accuracy

Changes by Wolfgang Maier : -- nosy: +wolma ___ Python tracker ___

[issue27761] Private _nth_root function loses accuracy

Changes by Ned Deily : -- priority: release blocker -> normal ___ Python tracker ___ ___

[issue27761] Private _nth_root function loses accuracy

Changes by Steven D'Aprano : -- versions: +Python 3.7 -Python 3.6 ___ Python tracker ___

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > Good enough for me ;-) Me too. -- ___ Python tracker ___ ___

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Let me clarify something about the extended algorithm: the starting guess is no longer the most significant source of error. It's the `mul(D(x), pow(D2, e))` part. `D(x)` is exact, but `pow(D2, e)` may not be exactly representable with 26 decimal digits, and

[issue27761] Private _nth_root function loses accuracy

Changes by Martin Panter : -- nosy: -martin.panter ___ Python tracker ___ ___

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Let me give complete code for the last idea, also forcing the scaling multiplication to use the correct context: import decimal c = decimal.DefaultContext.copy() c.prec = 25 c.Emax = decimal.MAX_EMAX c.Emin = decimal.MIN_EMIN def erootn(x,

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Oops! The `D2**e` in that code should be `pow(D2, e)`, to make it use the correct decimal context. -- ___ Python tracker

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Mark, thanks for the counterexample! I think I can fairly accuse you of thinking ;-) I expect the same approach would be zippy for scaling x by 2**e, provided that the scaled value doesn't exceed the dynamic range of the decimal context. Like so: def

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > And it's still the case that I haven't found a case where its result isn't > correctly rounded. Here's one: :-) >>> rootn(1 + 2**-52, 2) 1.0002 The correctly rounded result would be 1.0. -- ___

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > Decimal pow doesn't special-case integer exponents; the solution will still > be based on exp and ln. Ah, sorry; I'm wrong. That was true for the Python version of the decimal module, not for the C implementation. --

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > It does use the Decimal pow(), but with an integer exponent, so this specific > use of pow() doesn't invoke the Decimal exp() or ln() either. Decimal pow doesn't special-case integer exponents; the solution will still be based on exp and ln. --

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > Mark, the code I showed in roots.py is somewhat more accurate and highly > significantly faster than the code you just posted. Okay, fair enough. In that case, we still need a solution for computing rootn(x * 2**e) in the case where x*2**e itself overflows

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Mark, the code I showed in roots.py is somewhat more accurate and highly significantly faster than the code you just posted. It's not complicated at all: it just uses Decimal to do a single Newton correction with extended precision. Since it doesn't use the

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: Here's the error analysis, for the record. It's crude, but effective. Assume our decimal context precision is p (so p = 20 for the above code). 1. d is represented exactly with value in [1.0, 2.0). (Conversions from float to Decimal are exact.) 2. c.ln(d)

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: I think this whole nth root discussion has become way more complicated than it needs to be, and there's a simple and obvious solution. Two observations: 1. What we actually need is not nth_root(x), but nth_root(x*2**e) for a float x and integer e. That's

[issue27761] Private _nth_root function loses accuracy

Steven D'Aprano added the comment: As discussed with Ned via email, this issue shouldn't affect the public geometric_function and the failing tests (currently skipped) are too strict. The existing tests demand an unrealistic precision which should be loosened, e.g. from assertEqual to

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: BTW, add this other way of writing a native-precision Newton step to see that it's much worse (numerically) than writing it in the "guess + small_correction" form used in roots.py. Mathematically they're identical, but numerically they behave differently: def

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Attched file "roots.py" you can run to get a guess as to how bad pow(x, 1/n) typically is on your box. Note that it's usually "pretty darned good" the larger `n` is. There's a reason for that. For example, when n=1000, all x satisfying 1 <= x < 2**1000 have a

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Let's spell one of these out, to better understand why sticking to native precision is inadequate. Here's one way to write the Newton step in "guess + relatively_small_correction" form: def plain(x, n): g = x**(1.0/n) return g - (g -

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: As I said, the last code I posted is "fast enough" - I can't imagine a real application can't live with being able to do "only" tens of thousands of roots per second. A geometric mean is typically an output summary statistic, not a transformation applied

[issue27761] Private _nth_root function loses accuracy

STINNER Victor added the comment: Tim Peters: "Victor, happy to add comments, but only if there's sufficient interest in actually using this." It looks like tests fail on some platforms. I guess that it depends on the different factors: quality of the C math library, quality of the IEEE 754

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Steven, you certainly _can_ ;-) check first whether `r**n == x`, but can you prove `r` is the best possible result when it's true? Offhand, I can't. I question it because it rarely seems to _be_ true (in well less than 1% of the random-ish test cases I tried).

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: That's clever, Serhiy! Where did it come from? It's not Newton's method, but it also appears to enjoy quadratic convergence. As to speed, why are you asking? You should be able to time it, yes? On my box, it's about 6 times slower than the last code I posted.

[issue27761] Private _nth_root function loses accuracy

Serhiy Storchaka added the comment: Wouldn't following implementation be faster? import decimal c = decimal.DefaultContext.copy() c.prec = 25 def rootn(x, n, D=decimal.Decimal, sub=c.subtract, mul=c.multiply, log=c.ln):

[issue27761] Private _nth_root function loses accuracy

Steven D'Aprano added the comment: This has been really eye-opening, and I just wanted to drop you a note that I am watching this thread carefully. My first priority is to get the tests all passing before beta 1 on 2016-09-12, even if (as seems likely) that means weakening the tests, and then

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Victor, happy to add comments, but only if there's sufficient interest in actually using this. In the context of this issue report, it's really only important that Mark understands it, and he already does ;-) For example, it starts with float `**` because that's

[issue27761] Private _nth_root function loses accuracy

STINNER Victor added the comment: Hum, "g = D(x**(1.0/n))" computes a float and then convert it to a decimal, right? Can we please add comments to explain why you start with float, compute decimal and finish with float? And maybe also document the algorithm? --

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Adding one more version of the last code, faster by cutting the number of extra digits used, and by playing "the usual" low-level CPython speed tricks. I don't claim it's always correctly rounded - although I haven't found a specific case where it isn't. I do

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Serhiy, I don't know what you're thinking there, and the code doesn't make much sense to me. For example, consider n=2. Then m == n, so you accept the initial `g = x**(1.0/n)` guess. But, as I said, there are cases where that doesn't give the best result,

[issue27761] Private _nth_root function loses accuracy

Serhiy Storchaka added the comment: What if use pow() with exactly represented degree in approximating step? def rootn(x, n): g = x**(1.0/n) m = 1 << (n-1).bit_length() if n != m: g = (x*g**(m-n))**(1.0/m) return g Maybe it needs several

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: I don't care about correct rounding here, but it is, e.g., a bit embarrassing that >>> 64**(1/3) 3.9996 Which you may or may not see on your box, depending on your platform pow(), but which you "should" see: 1/3 is not a third, it's a binary

[issue27761] Private _nth_root function loses accuracy

STINNER Victor added the comment: > floor_nroot does arithmetic with integers of bit-length approximately 54*n. Oh, I see. But maybe the Decimal is fast enough for such giant numbers, since we can control the precision? -- ___ Python tracker

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: Victor: by "way too slow", I really *do* mean way too slow. :-) floor_nroot does arithmetic with integers of bit-length approximately 54*n. For small n, that's fine, but if someone tried to take the geometric mean of a list of 5 values (which it seems to

[issue27761] Private _nth_root function loses accuracy

STINNER Victor added the comment: Just to share my little experience with rounding numbers. Last years, I worked on a API to convert timestamps between float and integers, the private "PyTime C API": https://haypo.github.io/pytime.html At the beginning, I used various floatting point

[issue27761] Private _nth_root function loses accuracy

STINNER Victor added the comment: About tradeoff, would it be possible to add an option to choose the quality of the accuracy? For example, a flag to choose between "fast nth root" or "accurate nth root". Python already has two kind of numbers: Decimal and float. Maybe the "flag" should be

[issue27761] Private _nth_root function loses accuracy

STINNER Victor added the comment: "Just for fun, here's a recipe for a correctly-rounded nth root operation for positive finite floats. I'm not suggesting using this in the business logic: it's likely way too slow (especially for large n), but it may have a use in the tests." I don't know

[issue27761] Private _nth_root function loses accuracy

Changes by STINNER Victor : -- nosy: +haypo ___ Python tracker ___ ___

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Noting that `floor_nroot` can be sped a lot by giving it a better starting guess. In the context of `nroot`, the latter _could_ pass `int(x**(1/n))` as an excellent starting guess. In the absence of any help, this version figures that out on its own; an

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: Thanks, Mark! I had worked out the `floor_nroot` algorithm many years ago, but missed the connection to the AM-GM inequality. As a result, instead of being easy, proving correctness was a pain that stretched over pages. Delighted to see how obvious it _can_

[issue27761] Private _nth_root function loses accuracy

Ned Deily added the comment: Since test_NTh_Root is failing on so many different buildbot platforms and in different ways, I'm marking this as a "release blocker". It would be very helpful to get this resolved prior to 3.6.0b1 next month. Suggest checking the test results outputs of the

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > for positive finite floats ... assuming IEEE 754 binary64 format, of course. -- ___ Python tracker ___

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: Just for fun, here's a recipe for a correctly-rounded nth root operation for positive finite floats. I'm not suggesting using this in the business logic: it's likely way too slow (especially for large n), but it may have a use in the tests. The logic for the

[issue27761] Private _nth_root function loses accuracy

Tim Peters added the comment: A meta-note: one iteration of Newton's method generally, roughly speaking, doubles the number of "good bits" in the initial approximation. For floating n'th root, it would an astonishingly bad libm pow() that didn't get more than half the leading bits in pow(x,

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > I thought IEEE 754 was supposed to put an end to these sorts of > cross-platform differences. Maybe one day. (Though I'm not holding my breath.) IEEE 754-2008 does indeed *recommend* (but not actually *require*) correctly rounded implementations of all the

[issue27761] Private _nth_root function loses accuracy

Steven D'Aprano added the comment: On Sun, Aug 14, 2016 at 12:05:37PM +, Mark Dickinson wrote: > But I don't think there's a real problem here so long as you don't > have an expectation of getting super-accurate (e.g., correctly rounded > or faithfully rounded) results; testing that

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > What else can I do? Since I'm only dealing with integer powers, should I > try using my own ipow(y, n) for testing? I'd expect that a square-and-multiply-style pow would be less accurate than math.pow, in general, simply because of the number of

[issue27761] Private _nth_root function loses accuracy

Steven D'Aprano added the comment: On Sun, Aug 14, 2016 at 10:52:42AM +, Mark Dickinson wrote: > Same deal here: those aren't the actual errors; they're approximations > to the errors, since the computations of the epsilons depends on (a) > the usual floating-point rounding, and more

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > - finally, compare the epsilons abs(y**n - x) for the initial guess > and improved version, returning whichever gives the smaller epsilon. > > So I'm confident that nth_root() should never be worse than pow(). Same deal here: those aren't the actual

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: > - if y**n == x, return y; Here's the problem, though: you're using the libm pow again in this test, so the result of this being True doesn't give tight information about how close y is to the nth root of x (and the test result may well differ from platform

[issue27761] Private _nth_root function loses accuracy

Steven D'Aprano added the comment: On Sun, Aug 14, 2016 at 08:29:39AM +, Mark Dickinson wrote: > Steven: can you explain why you think your code *should* be giving > exact results for exact powers? Do you have an error analysis that > says that should be the case? No error analysis, only

[issue27761] Private _nth_root function loses accuracy

Mark Dickinson added the comment: Steven: can you explain why you think your code *should* be giving exact results for exact powers? Do you have an error analysis that says that should be the case? One issue here is that libm pow functions vary hugely in quality, so any algorithm that

[issue27761] Private _nth_root function loses accuracy

Martin Panter added the comment: The problem is more widespread than just Power PC. The same failures (+/-29) are also seen on the three s390x buildbots. There are multiple failures of testExactPowers(Negatives) on

[issue27761] Private _nth_root function loses accuracy on PowerPC

Raymond Hettinger added the comment: You can do no better than to consult with Uncle Timmy. -- nosy: +rhettinger, tim.peters ___ Python tracker ___

[issue27761] Private _nth_root function loses accuracy on PowerPC

Steven D'Aprano added the comment: Tests fail on a Power PC buildbot: http://buildbot.python.org/all/builders/PPC64LE%20Fedora%203.x/builds/1476/steps/test/logs/stdio == FAIL: testExactPowers

[issue27761] Private _nth_root function loses accuracy on PowerPC

New submission from Steven D'Aprano: First reported by Martin Panter here: http://bugs.python.org/issue27181#msg272488 I'm afraid I don't know enough about PowerPC to suggest a fix other than weakening the test on that platform. -- assignee: steven.daprano messages: 272638 nosy: