2

This amazing code golf answer to Is this number Loeschian? in its entirety:

Python, 49 bytes

lambda n:0in[(n-3*i*i+0j)**.5%1for i in range(n)]

Uses the equivalent quadratic form given on OEIS of n == 3*i*i+j*j. Check whether n-3*i*i is a perfect square for any i by taking its square root and checking if it's an integer, i.e. equals 0 modulo 1. Note that Python computes square roots of perfect squares exactly, without floating point error. The +0j makes it a complex number to avoid an error on the square root of a negative.

How does Python do this? Does **.5 "detect" that a given number is a perfect square somehow? Is this only reliable for integer input or will it work on floats up to some size as well?

I've also added a parenthetical Why? to the question; is this something that programmers rely upon? Is it for speed? Does it come with a cost?

Community
  • 1
  • 1
uhoh
  • 3,713
  • 6
  • 42
  • 95
  • 2
    It's required by the IEEE 754 standard. Square root results must be the closest double to the actual square root. For perfect squares, the closest double is an integer. – Robby Cornelissen Jun 08 '20 at 04:20
  • 1
    IEEE 754 does require `sqrt` to be correctly rounded, but there are an awful lot of gaps to fill in to get from there to `z**0.5` being correctly rounded in Python for a complex number `z` (even a complex number with zero imaginary part), and in general `z**0.5` will _not_ be correctly rounded. At least, it's certainly not on my Mac laptop, where `math.sqrt(x) == x**0.5` fails for many positive floats `x`. On a typical machine, you can trace a direct path from `math.sqrt` to the underlying correctly-rounded hardware sqrt, but the libm `pow` that underlies `x**0.5` is another matter entirely. – Mark Dickinson Jul 14 '20 at 18:00
  • BTW, I'm a bit puzzled about how the solution is expected to work, since `z%1` isn't valid for a complex number `z`. `(lambda n:0in[(n-3*i*i+0j)**.5%1for i in range(n)])(35)` gives me `TypeError: can't mod complex numbers.` – Mark Dickinson Jul 15 '20 at 08:34
  • @MarkDickinson maybe you can leave a comment [there](https://codegolf.stackexchange.com/a/88796/85527)? – uhoh Jul 15 '20 at 12:24
  • @uhoh: Good idea; will look at that later when I have time. (I don't think I even have a login there at the moment.) – Mark Dickinson Jul 15 '20 at 14:41

1 Answers1

2

You can check out the source code here. They describe the algorithm they use for computing the (approximate) square root of nonnegative integers, and show that for perfect squares the algorithm gives the exact answer. The code is C, but they give a translation of the code into Python:

def isqrt(n):
    """
    Return the integer part of the square root of the input.
    """
    n = operator.index(n)
    if n < 0:
        raise ValueError("isqrt() argument must be nonnegative")
    if n == 0:
        return 0
    c = (n.bit_length() - 1) // 2
    a = 1
    d = 0
    for s in reversed(range(c.bit_length())):
        # Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2
        e = d
        d = c >> s
        a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a
    return a - (a*a > n)

I assume, but haven't yet checked, that when computing a power at runtime, Python first checks that 1. the base is a nonnegative integer, 2. the exponent is exactly 0.5, and if those both hold then it invokes the code I linked to above.

BallpointBen
  • 9,406
  • 1
  • 32
  • 62
  • I'm afraid your assumption isn't correct. For a float `x`, Python (at least CPython) passes an `x**0.5` operation straight down to the underlying C math library. It's *conceivable* that that math library could special-case an exponent of `0.5` to use the CPU's `sqrt` instruction (if such exists), which one can reasonably expect to be correctly rounded, but in my experience that doesn't usually happen. It certainly doesn't on my Mac laptop: `x**0.5` and `math.sqrt(x)` give observably different results for many `x`. – Mark Dickinson Jul 14 '20 at 17:49
  • It wouldn't make sense to substitute `isqrt(x)` for `x**0.5`, anyway: if `x` is not a perfect square, the two will usually have quite different values. For example, `isqrt(8)` is `2`, but `8**0.5` is `2.8284271247461903`. – Mark Dickinson Jul 14 '20 at 17:51