Support negative exponents in pow() where a modulus is specified. #80208
Comments
|
Having gcd() in the math module has been nice. Here is another number theory basic that I've needed every now and then: def multinv(modulus, value):
'''Multiplicative inverse in a given modulus >>> multinv(191, 138)
18
>>> 18 * 138 % 191
1
>>> multinv(191, 38)
186
>>> 186 * 38 % 191
1
>>> multinv(120, 23)
47
>>> 47 * 23 % 120
1 |
rhettinger
added
3.8
x = (a+b) % m
x = a*b % m
x = pow(a, b, m)or in math:
pow(value, -1, modulus)which currently raises an exception. Presumably pow(value, -n, modulus)
for int n > 1 would mean the same as pow(pow(value, -1, modulus), n, modulus), if it were accepted at all. I'd be happy to stop with -1.
But I'm not sure anyone would use that _except_ to compute modular inverse. So probably not. |
Yes, that makes sense.
+1 ;-) |
|
+1 for the pow(value, -1, modulus) spelling. It should raise It would feel odd to me not to extend this to |
|
Here's an example of some code in the standard library that would have benefited from the availability of Lines 957 to 960 in e7a4bb5 if self._exp >= 0:
exp_hash = pow(10, self._exp, _PyHASH_MODULUS)
else:
exp_hash = pow(_PyHASH_10INV, -self._exp, _PyHASH_MODULUS)where: _PyHASH_10INV = pow(10, _PyHASH_MODULUS - 2, _PyHASH_MODULUS)With the proposed addition, that just becomes |
|
Agreed, extending pow(value, n, modulus) to negative n would be a great addition! To have modinv(value, modulus) next to that also makes a lot of sense to me, as this would avoid lots of confusion among users who are not so experienced with modular arithmetic. I know from working with generations of students and programmers how easy it is to make mistakes here (including lots of mistakes that I made myself;) One would implement pow() for negative n, anyway, by first computing the modular inverse and then raising it to the power -n. So, to expose the modinv() function to the outside world won't cost much effort. Modular powers, in particular, are often very confusing. Like for a prime modulus p, all of pow(a, -1,p), pow(a, p-2, p), pow(a, -p, p) are equal to eachother, but a common mistake is to take pow(a, p-1, p) instead. For a composite modulus things get much trickier still, as the exponent is then reduced in terms of the Euler phi function. And, even if you are not confused by these things, it's still a bit subtle that you have to use pow(a, -1,p) instead of pow(a, p-2, p) to let the modular inverse be computed efficiently. With modinv() available separately, one would expect --and get-- an efficient implementation with minimal overhead (e.g., not implemented via a complete extended-gcd). |
That's not 100% clear: the binary powering algorithm used to compute |
OK, I'm indeed assuming that modinv() is implemented efficiently, in CPython, like pow() is. Then, it should be considerably faster, maybe like this: >>> timeit.timeit("gmpy2.invert(1023,p)", "import gmpy2; p=2**61-1")
0.18928535383349754
>>> timeit.timeit("gmpy2.invert(1023,p)", "import gmpy2; p=2**127-1")
0.290736872836419
>>> timeit.timeit("gmpy2.invert(1023,p)", "import gmpy2; p=2**521-1")
0.33174844290715555
>>> timeit.timeit("gmpy2.powmod(1023,p-2,p)", "import gmpy2; p=2**61-1")
0.8771009990597349
>>> timeit.timeit("gmpy2.powmod(1023,p-2,p)", "import gmpy2; p=2**127-1")
3.460449585430979
>>> timeit.timeit("gmpy2.powmod(1023,p-2,p)", "import gmpy2; p=2**521-1")
84.38728888797652 |
Why would you expect that? Both algorithms involve a number of (bigint) operations that's proportional to log(p), so it's going to be down to the constants involved and the running times of the individual operations. Is there a clear reason for your expectation that the xgcd-based algorithm should be faster? Remember that Python has a subquadratic multiplication (via Karatsuba), but its division algorithm has quadratic running time. |
Yeah, good question. Maybe I'm assuming too much, like assuming that it should be faster;) It may depend a lot on the constants indeed, but ultimately the xgcd style should prevail. So the pow-based algorithm needs to do log(p) full-size muls, plus log(p) modular reductions. Karatsuba helps a bit to speed up the muls, but as far as I know it only kicks in for quite sizeable inputs. I forgot how Python is dealing with the modular reductions, but presumably that's done without divisions. The xgcd-based algorithm needs to do a division per iteration, but the numbers are getting smaller over the course of the algorithm. And, the worst-case behavior occurs for things involving Fibonacci numbers only. In any case, the overall bit complexity is quadratic, even if division is quadratic. There may be a few expensive divisions along the way, but these also reduce the numbers a lot in size, which leads to good amortized complexity for each iteration. |
I'll work up a PR using the simplest implementation. Once that's in with tests and docs, it's fair game for someone to propose algorithmic optimizations. |
|
Raymond, I doubt we can do better (faster) than straightforward egcd without heroic effort anyway. We can't even know whether a modular inverse exists without checking whether the gcd is 1, and egcd builds on what we have to do for the latter anyway. Even if we did know in advance that a modular inverse exists, using modular exponentiation to find it requires knowing the totient of the modulus, and computing the totient is believed to be no easier than factoring. The only "optimization" I'd be inclined to _try_ for Python's use is an extended binary gcd algorithm (which requires no bigint multiplies or divides, the latter of which is especially sluggish in Python): https://www.ucl.ac.uk/~ucahcjm/combopt/ext_gcd_python_programs.pdf For the rest:
|
|
Changing the title to reflect a focus on building-out pow() instead of a function in the math module. |
rhettinger
changed the title
|
In pure Python this seems to be the better option to compute inverses: def modinv(a, m): # assuming m > 0
b = m
s, s1 = 1, 0
while b:
a, (q, b) = b, divmod(a, b)
s, s1 = s1, s - q * s1
if a != 1:
raise ValueError('inverse does not exist')
return s if s >= 0 else s + mBinary xgcd algorithms coded in pure Python run much slower. |
|
I think #57475 is ready to go, but I'd appreciate a second pair of eyes on it if anyone has time. |
|
Done. Closing. |
|
PR 13266 introduced a compiler warning. Objects/longobject.c: In function ‘long_invmod’:
Objects/longobject.c:4246:25: warning: passing argument 2 of ‘long_compare’ from incompatible pointer type [-Wincompatible-pointer-types]
if (long_compare(a, _PyLong_One)) {
^~~~~~~~~~~
Objects/longobject.c:3057:1: note: expected ‘PyLongObject * {aka struct _longobject *}’ but argument is of type ‘PyObject * {aka struct _object *}’
long_compare(PyLongObject *a, PyLongObject *b)
^~~~~~~~~~~~ |
|
I will fix the compiler warning along with another one that I just introduced. |
|
@petr: Thanks for the quick fix! |
|
For tracker historians: see also bpo-457066 |
|
https://bugs.python.org/issue457066 The old is new again ;-). |
Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.
Show more details
GitHub fields:
bugs.python.org fields:
The text was updated successfully, but these errors were encountered: