By default, isqrt(n) gives the floor of the exact square of n.  It would be nice to have a flag to give a rounded result:

    y = isqrt(n, round=True)

Alternatively, set a mode argument to one of {'floor', 'round', 'ceil'}:

    y = isqrt(n, mode='round')

I would like something better than this:

    def risqrt(x):
        'Big integer version of: round(sqrt(x)).'
        y = isqrt(x)
        s = y ** 2
        return y if x <= s + y else y + 1

    def cisqrt(x):
        'Big integer version of: ceil(sqrt(x)).'
        return isqrt(x - 1) + 1

My use case arose when building a table of square roots incorporated in arbitrary precision functions implemented with scaled integer arithmetic:

    def get_root_table(base, steps, scale):
        s = []
        x = round(base * scale)
        for i in range(steps):
            x = risqrt(x * scale)
            s.append(x)
        return s

FWIW, when this need has turned up for me (which it has, a couple of times), I've used this:

def risqrt(n):
    return (isqrt(n<<2) + 1) >> 1

But I'll admit that that's a bit non-obvious.

[Mark]
> def risqrt(n):
>    return (isqrt(n<<2) + 1) >> 1

Sweet! New one on me - did you invent this? It's slick :-)

I'd be happy to see recipes added to the docs for rounded and ceiling flavors of isqrt, but am dubious about the value of building them in.

If they are added via an optional rounding argument, can we use the decimal module's names for the supported rounding modes?

> Sweet! New one on me 

Tim already knows this but for the record the derivation is isn't tricky.

With y=isqrt(x), then next root is at y+1 and the half way point is y+1/2 which isn't an integer.  The corresponding squares are y**2, (y+1/2)**2, and (y+1)**2.  With a multiple of four, those become 4*y**2, 4*(y+1/2)**2, and 4*y+1)**2 which are equivalent to the squares of consecutive integers: (2y)**2, (2y+1)**2, and (2y+2)**2.  Adding one gives the roots 2y+1, 2y+2, and 2y+3.  The floor division eliminates the constant term giving the desired roots y, y+1, and y+1.

> can we use the decimal module's names for the supported rounding modes?

I'm not sure those make sense because we never get to exactly half.  There is only floor, ceil, and round, not half_up, half_even, etc.

>> can we use the decimal module's names for the supported
>> rounding modes?

> I'm not sure those make sense because we never get to
> exactly half.  There is only floor, ceil, and round,
> not half_up, half_even, etc.

So use decimal's ROUND_CEILING, ROUND_FLOOR, and ROUND_HALF_EVEN. It's irrelevant that the halfway case can't occur: it's still following the ROUND_HALF_EVEN rules, it's just that one of those rules never happens to apply in this context. So what? Your _intent_ is to supply "best possible rounding", and that's what ROUND_HALF_EVEN means. It's not doing any favors to make up a new name for a rounding mode that only applies to values that can never tie. Tail. dog, wag ;-) If someone knows what half/even does, they already know what _this_ rounding mode does. Why complicate it with a useless distinction?

FYI, I had a simpler derivation in mind. Say

sqrt(n) = r + f

where r = isqrt(n) and 0 <= f < 1. Then

sqrt(4n) = 2 * sqrt(n) = 2*(r + f) = 2r + 2f, with 0 <= 2f < 2.

If f < 0.5, 2f < 1, so isqrt(4n) = 2r, and we shouldn't round r up either.

If f > 0.5, 2f > 1, so sqrt(4n) = 2r + 1 + (2f - 1), with 0 <= 2f - 1 < 1, so isqrt(4n) = 2*r + 1. In this case (f > 0.5) we need to round r up.

f = 0.5 can't happen.

Regardless, I don't believe I would have thought of this myself! It was an unexpected delight :-

> So use decimal's ROUND_CEILING, ROUND_FLOOR, and ROUND_HALF_EVEN

It think is would suck to have to type those out.  Sorry, I think you're headed down the path of foolish consistency with an unrelated module and a more complicated topic.  What's wrong with keeping consistent the name of the existing functions: round, floor, and ceil which are self-explanatory and much better known than the decimal module constants.  Also, the intent for these to be substitutable for existing code which uses round(sqrt(x)), floor(sqrt(x)), and ceil(sqrt(x)).  If those are the starting point, then round, floor, and ceil are a logical progression from the existing code.  It is also consistent with how the numeric tower approaches various ways of rounding.

All cool. Since I doubt these new rounding modes will get much use anyway, it's not worth arguing about. If I were to do this for myself, the rounding argument would be one of the three values {-1, 0, +1}, which are already more than good enough as mnemonics for "go down, get close, go up". Passing strings of any kind is already pedantic overkill ;-)

> did you invent this?

The idea is no more than: "compute an extra bit, then use that extra bit to determine which way to round". More generally, for any real number x, the nearest integer to x (rounding ties towards +infinity) is `⌊(⌊2x⌋ + 1) / 2⌋`. Now put x = √n, then ⌊2x⌋ is ⌊√(4n)⌋, and the rest follows.

> I'd be happy to see recipes added to the docs for rounded and ceiling flavors of isqrt, but am dubious about the value of building them in.

I'd similarly prefer to see recipes in the docs. We already have such a recipe for ceil(√n).

[Mark]
> The idea is no more than: "compute an extra bit, then use
> that extra bit to determine which way to round".

Thanks! Despite that this is what the derivation I sketched boiled down to, I was still missing how general the underlying principle was. Duh! Indeed, nothing special about sqrt here. So I'm delightfully surprised again :-)

> I'd similarly prefer to see recipes in the docs.

Can you elaborate on why you prefer having this in the docs rather than as built-in functionality?

For me, the rounded case would be the most common.  I don't think I'm better-off writing a wrapper with a semi-opaque trick, having to give it a different name and having it run more slowly than built in functionality. Also as builtin, I can be assured that it is maintained, documented, and well-tested.

Writing "y = isqrt(x, 'rounded')" is clear, fast, and convenient. I don't see the advantage is putting that functionality a little more out of reach.

The string options have been successful elsewhere in Python. For example, str.encode() has the error modes: "strict", "ignore", and "replace" which as clear, fast, and much better than users having to write their own wrappers.

Pretend this were the itertools module and you were responding to someone suggesting a new optional argument to one of its functions. You have tons of experience rejecting those, and much the same arguments apply, from "not every one-line function needs to be built in" to "in seventeen centuries this is the first time someone asked for it" ;-)

Seriously, the only complaint I've seen about isqrt() before is that there isn't a version that raises an exception if the argument isn't a perfect square. That's so niche I wouldn't support adding that to the core either. Sure, sometimes that's what someone may want, and that's fine - they can do it on their own.

So, are you going to reject my feature request?  I don't understand why.  Both Mark and I have had valid use cases.  The implementation is straight-forward and simple.  The workaround is slow and non-obvious.  The default remains the same so that usability isn't impaired in any way.  There is basically zero downside other that feeling icky about mode flags in general.

With respect to itertools, I did add *initial* to accumulate() solely because you wanted it ;-)

I'm not rejecting anything. Mark wrote the function, and I'll defer to him.

Yes, you added `initial` to `accumulate()`, and I spent hours in all explaining to you why it was an utterly natural addition, conspicuous by absence, including writing up several real-life use cases that just happened to pop up while that discussion dragged on. Thanks! I've used it many more times since it was added. Adding it brought accumulate() into line with other languages' workalikes. Python was the oddball there.

OTOH, not even gmp's mpz supports rounded isqrt directly, just truncated, and a _different_ function that returns not only isqrt(n), but also "the remainder" (n - isqrt(n)**2). The latter shows they're not averse to adding cruft that can easily be done - which makes it all the more suggestive that they _didn't_ also add other rounding modes.

Note: the primary intended "use case" for the version with the remainder appears to be a fast way to detect if the argument was a perfect square. As mentioned earlier, the lack of a builtin way to do that is the only complaint I've seen before about Python's isqrt.

Mark didn't mention his use case for rounded isqrt, I don't have any myself, and your only named use case was "building a table of square roots incorporated in arbitrary precision functions implemented with scaled integer arithmetic". It's very hard to believe that minor speed boosts make any difference at all when building static tables, so "it's faster built in" falls flat for me.

To the contrary, adding _any_ optional argument (be it a string or of any other type) complicates the calling sequence and requires burning cycles on every invocation just to work out which arguments _were_ passed. So, despite that I have no expectation of ever asking for a rounded isqrt, adding the possibility would cost _me_ cycles on every isqrt() call I already have. None of which are in one-shot table-builders, but mostly in loops.

If there's a _compelling_ use case for rounded isqrt, I haven't seen one, but I'd be much happier to see it added it as a new function of its own. Like, e.g., gmp added mpz_sqrtrem() rather than add an annoying "mode" argument to the widely used mpz_sqrt().

I don't really object to string arguments per se, but in the context of what one hopes to be blazing fast integer functions, throwing in string comparisons just to set a flag seems ludicrously expensive. _PyUnicode_EqualToASCIIId() is not my idea of "cheap" ;-)

Suppose we added isqrt_rem(n), returning the integer pair (i, rem) such that

    n == i**2 + rem
    0 <= rem <= 2*i

Then:

- Want the floor of sqrt(n)? i.
- The ceiling? i + (rem != 0).
- Rounded? i + (rem > i).
- Is n a perfect square? not rem.

That's how mpz addresses these, although it has a different function to compute the floor without returning the remainder too.

I wouldn't object to that - just +0, though. Depending on implementation details, which I haven't investigated, it may or may not be materially faster than doing:

def isqrt_rem(n):
    return (i := isqrt(n)), n - i*i

myself.

A new function isqrt_rem seems like a reasonably natural addition. (Though I'd call it "isqrtrem", partly by analogy with "divmod", and partly because the math module isn't very good at doing underscores.)

> Mark didn't mention his use case for rounded isqrt

Mainly for emulation of floating-point sqrt. But the number of times I've needed rounded integer square root is small compared with the number of times I've needed rounded integer division.

The name "isqrtrem" works fine for me too, and, indeed, is even more reminiscent of mpz's workalike function.

> the number of times I've needed rounded integer square root
> is small compared with the number of times I've needed rounded
> integer division

Speaking of which, the analogy to divmod extends just beyond the lack of an underscore in the name ;-) divmod() allows easy emulation of any division rounding mode, and an instant answer to "was the division exact?". isqrtrem() would support the same for square roots.

FWIW, gmpy2 uses isqrt_rem. Of course, gmpy2 uses underscores a lot. And uses next_toward instead of nextafter, and copy_sign instead of copysign, and is_nan instead of isnan... :(

I would keep the math module consistent and use isqrtrem.

I'll look at adding aliases to remain consistent with the math module.

> divmod() allows easy emulation of any division rounding mode

It could be used that way, but generally isn't.

Please consider my original request.  Adding a keyword argument is easy, clear, and has almost no mental overhead.   

It reads very well in code, `y = isqrt(x, 'ceil')` or `y = isqrt(x, 'round')`.  The equivalents with isqrt_rem are awkward and don't read well (reminding me of my Fortran days).

I've made several good-faith efforts to find any hint of demand for rounded isqrt on the web; I've found no direct support for it in other languages/environments(*); and the one use case you presented isn't compelling. Building static tables to help implement extended precision functions via scaled integer arithmetic is quite niche, and anyone numerically knowledgeable enough to pull that off will have no trouble figuring out how to roll their own isqrt rounding. For them, isqrtrem would make it dead easy.

Close to no demonstrated demand, and no demonstrated prior art, coupled with slowing down _every_ isqrt() call to cater to an overwhelmingly unused possibility (and, to boot, via the slowest possible way of specifying an option, on what should be a fast function), leaves me -1. No, I'm not a fan of the "big"/"little" mandatory arguments to int.{from,to}_bytes() either.

+0 on isqrtrem, though, because it follows established practice in a related context (mpz), and also offers efficient support for the related use case I actually found some (albeit small) demand for ("and was it a perfect square?").

(*) Other languages with an integer square root include Common Lisp, Julia, and Maple's programming language. None have an option for rounding. Maple's doesn't define the rounding used, other than to promise the result is less than 1 away from the infinitely precise result. Common Lisp and Julia return the floor. A number of environments essentially inherit isqrt from their Lisp base (e.g., wxMaxima).

Is

   i, rem = isqrt_rem(n)
   i + (rem != 0)

better than

   (isqrt(n<<2) + 1) >> 1

or

   n and isqrt(n-1) + 1

?

As for use cases, there were few cases in my experience when I needed the ceiling square root, mostly in simple experiments in REPL, but it was so rary, and workarounds satisfied me. So it would be a nice to have feature which I would use perhaps once in a year or two years, but can live without it. And I do not want to pay a cost of significantly complicating API or slowing down isqrt().

> Is
>
>  i, rem = isqrt_rem(n)
>  i + (rem != 0)
>
> better than
>
>   (isqrt(n<<2) + 1) >> 1
>
> or
>
>   n and isqrt(n-1) + 1
>
> ?

Define "better"? The first way is by far the most obvious of the three, and the second way the least obvious. The first way also "wins" on being  a variation of a uniform pattern that can deliver the floor, ceiling, or rounded result, depending on which simple comparison result is added. It's not "a trick" - it's the opposite of clever ;-)

The first way is also unique in being the only one of the three that does _not_ do any Python-level arithmetic on integers as wide as `n`. `i` and `rem` are no more than about half the bit length of `n`. `n << 2` and `n - 1` in the others have to create new int objects at least as wide as `n`.

I've been keeping my eyes open. The only mariginally relevant thing I've noticed is Hart's "one line factoring" method:

http://wrap.warwick.ac.uk/54707/1/WRAP_Hart_S1446788712000146a.pdf

That wants the _ceiling_ of the square root. And in another place it wants to know "is it a perfect square?". But the part that wants the ceiling doesn't care whether it's a perfect square, while the part asking "is it a perfect square?" doesn't care what the ceiling may be.

Thanks, Tim; very interesting. I hadn't seen this factoring algorithm before.

> That wants the _ceiling_ of the square root.

Looks like what it actually wants is the ceiling analog of isqrtrem: that is, it needs both the ceiling of the square root *and* the difference between the square of that ceiling and the original number.

The description of the algorithm in section 2 is a bit odd: they define m := s*s % n, using an expensive modulo operation, when all they need is a subtraction: m := s*s - n*i. This is noted in section 3 ("to reduce modulo Mn at step 3, one may simply subtract Mni from s2"), but they fail to note that the two things aren't equivalent for large enough i, possibly because that large an i won't be used in practice. And in the case that the two quantities differ, it's the subtraction that's needed to make the algorithm work, not the mod result.

Here's a Python version of Hart's algorithm:


from itertools import count
from math import gcd, isqrt

def isqrtrem(n):
    """ For n >= 0, return s, r satisfying s*s + r == n, 0 <= r <= 2*s. """
    s = isqrt(n)
    return s, n - s*s

def csqrtrem(n):
    """ For n > 0, return s, r satisfying n + s*s == r, 0 <= r <= 2*(s-1). """
    s = 1 + isqrt(n-1)
    return s, s*s - n

def factor(n):
    """ Attempt to use Hart's algorithm to find a factor of n. """
    for i in count(start=1):
        s, m = csqrtrem(n*i)
        t, r = isqrtrem(m)
        if not r:
            return gcd(n, s-t)

Ah, https://math.mit.edu/research/highschool/primes/materials/2019/Gopalakrishna.pdf is interesting - they conjecture a bound on the number of iterations required, and note that under that conjecture the mod can be replaced by a subtraction.