[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

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)

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

> 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`.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Serhiy Storchaka]

Serhiy Storchaka  added the comment:

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().

--
nosy: +serhiy.storchaka

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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).

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Raymond Hettinger]

Raymond Hettinger  added the comment:

> 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).

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Case Van Horsen]

Case Van Horsen  added the comment:

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.

--
nosy: +casevh

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

> 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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

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.)

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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" ;-)

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Raymond Hettinger]

Raymond Hettinger  added the comment:

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 ;-)

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Raymond Hettinger]

Raymond Hettinger  added the comment:

> 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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

[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 :-)

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

> 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).

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

> 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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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 ;-)

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Raymond Hettinger]

Raymond Hettinger  added the comment:

> 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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

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 :-

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

>> 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?

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Raymond Hettinger]

Raymond Hettinger  added the comment:

> 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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Tim Peters]

Tim Peters  added the comment:

[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?

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Mark Dickinson]

Mark Dickinson  added the comment:

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.

--

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com

[issue46187] Optionally support rounding for math.isqrt()

[Raymond Hettinger]

New submission from Raymond Hettinger :

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

--
assignee: mark.dickinson
components: Library (Lib)
messages: 409243
nosy: mark.dickinson, rhettinger, tim.peters
priority: normal
severity: normal
status: open
title: Optionally support rounding for math.isqrt()
type: enhancement
versions: Python 3.11

___
Python tracker 

___
___
Python-bugs-list mailing list
Unsubscribe: 
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com