[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

On 8/4/2017 11:44 AM, Chris Angelico wrote: > On Sat, Aug 5, 2017 at 12:51 AM, Steve D'Aprano > <steve+python at pearwood.info> wrote: >> def isqrt_float(n): >> """Integer square root using floating point sqrt.""" >> return int(math.sqrt(n)) The operations in the integer version are well-defined and so the answer should be the same on any Python 3 and indeed in any language that does proper integer arithmetic (but most do not). In CPython, Math.sqrt wraps the C compiler's double sqrt function. To make it deterministic, we must specify the float representation, rounding rule, and accuracy of the sqrt function. Let's assume that the float answer is IEEE binary64, with 53 binary digits, with rounding rule 'round to nearest, ties to even' and that the algorithm is as exact as possible, which is to say, the exact mathematical answer rounded as specified. With 53 binary digits, all counts from 0 to 2**53 - 1 are exactly represented with a exponent of 0, 2**53 = 2**52 * 2, so it is exactly represented with an exponent of 1. Many other higher counts can be exactly represented with various exponents, but 2**53 + 1 requires 54 digits. So it is the first count that does not have an exact binary64 representation. Under these assumptions, how can isqrt_float(n) fail? Suppose 0 <= m <= 2**26.5, so that m*m <= 2**53. Under the assumptions that m.sqrt is as exact as possible, sqrt(m*m) == float(m) and int(float(m)) == m. The true square root of m*m -1 will be m - e (e < 1) and the integer square root is m-1. However, if e is less than (and I believe equal) to 1/2 of the interval between n and the next lowest float, m - e will be rounded up to float(m) and ifloat(m*m-1) will incorrectly return m instead of m-1. As m increases, e decreases while the interval between floats increases, so we should expect an eventual error In this range of m, sqrt(m*m +1) will always be greater than m, so rounding down cannot lead to an error. For larger m, where float(m*m) is inaccurate enough, it might. As m increases, e decreases while the interval between floats increases, so we should expect an eventual error. Whether it occurs for m <= 2* or not might be determined by careful analysis, but Chris Angelico already gave us the brute force answer. >> We know that: >> >> - for n <= 2**53, isqrt_float(n) is exact; but it is not ;-). In this range, the only possibility is math.sqrt(m*m-1) being m instead of m - delta and that is the case for Chris' answer 4503599761588224 = 67108865**2 - 1. int(math.float(4503599761588224) is 67108865 instead of 67108864. >> - for n >= 2**1024, isqrt_float(n) will raise OverflowError >> - between those two values, 2**53 < n < 2**1024, isqrt_float(n) >> will sometimes be exact, and sometimes not exact; >> - there is some value, let's call it M, which is the smallest >> integer where isqrt_float is not exact. >> Your mission, should you choose to accept it, is to find M. > Hmm. Thinking aloud a bit here. We know that isqrt_float(n) is not > technically *exact* for any n that is not a square. So what I'd do is > assume that for (n*n+1) to (n+1)*(n+1)-1, it's going to return the > same value. I would then probe every perfect square, and the values > one above and one below it. Now, that's still probably too many to > check, but it's going to be a lot faster; for starters, we don't > actually need to use isqrt_newton to compare against it. > > for root in range(2**26, 2**53): > square = root * root > if isqrt_float(square - 1) != root - 1: > print(square - 1) > break > if isqrt_float(square) != root: > print(square) > break > if isqrt_float(square + 1) != root: > print(square + 1) > break > > That gave me this result almost instantaneously: > > 4503599761588224 Your limit was barely low enough ;-). >>> math.log2(math.sqrt(4503599761588224)) 26.000000021497833 If you had started with, say int(2**26.5)-1, you would have missed this. I reran with range(2, 2**26) and there is nothing lower. >>> math.log2(4503599761588224) 52.00000004299566 >>> import decimal as d >>> I = d.Decimal(4503599761588224) >>> I.sqrt() Decimal('67108864.99999999254941951410') >>> float(I.sqrt()) 67108865.0 > which has been rounded up instead of down. I don't know if that counts > as sufficiently wrong? The isqrt_float answer is definitely wrong, and so is the claim for n < 2**53. >>>> isqrt_float(4503599761588224) > 67108865 >>>> isqrt_newton(4503599761588224) > 67108864 >>>> 67108865 ** 2 > 4503599761588225 -- Terry Jan Reedy

- Prev by Date:
**Challenge: find the first value where two functions differ** - Next by Date:
**PyYaml not using Yaml 1.2?** - Previous by thread:
**Challenge: find the first value where two functions differ** - Next by thread:
**Challenge: find the first value where two functions differ** - Index(es):