Implement modular exponentiation for all polynomial rings

[jpflori]

sage: x = polygen(ZZ)
sage: pow(x,100,x+1)
x^100

CC: @roed314 @jdemeyer @pjbruin @sagetrac-jakobkroeker

Component: algebra

Stopgaps: wrongAnswerMarker

Author: Shashwat Singh

Branch/Commit: u/gh-shashwat1002/mod_exponentation_polynom @ c37ccc8

Issue created by migration from https://trac.sagemath.org/ticket/15867

[sagetrac-jakobkroeker]

Stopgaps: wrongAnswerMarker

[shashwat1002]

Branch: u/gh-shashwat1002/mod_exponentation_polynom

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

ff56f10

Implement modular exponentiation in polynomial rings Flint.

[sagetrac-git]

Commit: ff56f10

[slel]

comment:8

Thanks! Please add your full name in the "Authors" field.

[yyyyx4]

comment:9

Now that the modulus argument is no longer ignored, you should rename it from ignored to something else (modulus? mod?).

[shashwat1002]

Author: Shashwat Singh

[sagetrac-git]

Changed commit from ff56f10 to c2a7b16

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

c2a7b16

Implement modular exponentiation in polynomial rings Flint.

[shashwat1002]

comment:12

Replying to @slel:

Thanks! Please add your full name in the "Authors" field.

Yeah, done! :)

[shashwat1002]

comment:13

Replying to @yyyyx4:

Now that the modulus argument is no longer ignored, you should rename it from ignored to something else (modulus? mod?).

renamed to modulus and rewrote history :)

[slel]

comment:14

To compare to None, we use is and is not
rather than == and !=:

-                if modulus != None:
+                if modulus is not None:

[sagetrac-git]

Changed commit from c2a7b16 to 8293154

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

8293154

Implement modular exponentiation in polynomial rings Flint.

[shashwat1002]

comment:16

Replying to @slel:

To compare to None, we use is and is not
rather than == and !=:

-                if modulus != None:
+                if modulus is not None:

Done (and rewrote history) :)

[yyyyx4]

comment:17

This patch makes things better in that it renders the implementation of __pow__ correct, but it is still kind of misses the point: The entire reason why pow() has a "modulus" argument in the first place is that pow(x, e, m) can be exponentially faster than pow(x, e) % m! So, having an implementation of "modular exponentiation" that is simply exponentiation followed by modulo is rather pointless. It should use square-and-multiply with reductions at each step to keep things small.

The other issue is that the patch only affects modular exponentation for one particular univariate polynomial class in Sage, but there are tons of other polynomial implementations that still don't implement this correctly. Presumably, the best way to go would be to implement a generic modular-exponentiation function (do we have that already?) and simply call it in each of the polynomial classes when the modulus argument is given.

[shashwat1002]

comment:18

Replying to @yyyyx4:

This patch makes things better in that it renders the implementation of __pow__ correct, but it is still kind of misses the point: The entire reason why pow() has a "modulus" argument in the first place is that pow(x, e, m) can be exponentially faster than pow(x, e) % m! So, having an implementation of "modular exponentiation" that is simply exponentiation followed by modulo is rather pointless. It should use square-and-multiply with reductions at each step to keep things small.

The other issue is that the patch only affects modular exponentation for one particular univariate polynomial class in Sage, but there are tons of other polynomial implementations that still don't implement this correctly. Presumably, the best way to go would be to implement a generic modular-exponentiation function (do we have that already?) and simply call it in each of the polynomial classes when the modulus argument is given.

I am going to begin work on a modular exponentiation function that is more efficient.

But before I begin, I had a few doubts:

• the polynomials are dealt by one of two C libraries: FLINT and NTL, so a general modular exponentiation function on polynomials must be written in python?
• Also, I am having trouble figuring where would I place the code for that, and understanding the organization of the polynomials code in general. If I could be pointed to a developer documentation on the same it will be very helpful

I started a mailing list thread here if that's more appropriate to respond to https://groups.google.com/g/sage-devel/c/5A8CDkj9jKo

[nbruin]

comment:19

If you're going to implement this on the level of __pow__ (which, according to python's interface specification is actually a reasonable thing to do), you will need to implement it for FLINT and NTL wrappers separately, since each implements __pow__, so any inheritance would be overridden.

Also note that "efficient" powering of polynomials modulo other polynomials over QQ is less useful than you might initially think: it helps a LOT to not carry the higher order terms around, but the coefficient swell will still be significant, bounding the feasible powers to a relatively small range. For polynomials over finite fields, where coefficient swell is not an issue, this is an entirely different story.

You probably also want to document what exactly the modulo operation that you're doing, since ZZ[x] by itself is not a Euclidean ring, so the concept of "remainder" may need further explanation. Are you using some kind of pseudo division?

As an example:

sage: R.<x>=ZZ['x']
sage: f=x^5-1
sage: m=3*x^2+1
sage: f % m
x^5 - 1
sage: (3*f) % m
-x^3 - 3
sage: (9*f) % m
x - 9

[yyyyx4]

comment:20

Replying to @nbruin:

You probably also want to document what exactly the modulo operation that you're doing, since ZZ[x] by itself is not a Euclidean ring, so the concept of "remainder" may need further explanation. Are you using some kind of pseudo division?

I think this is outside the scope here: The notion of reduction is simply whatever % does for the parent ring, which should not surprise anyone. Improving the documentation for %, if deemed insufficient, should be another ticket.

[yyyyx4]

comment:21

Replying to @shashwat1002:

• the polynomials are dealt by one of two C libraries: FLINT and NTL, so a general modular exponentiation function on polynomials must be written in python?

Yes. In fact, I discovered that we already have such a function in Sage: sage.arith.misc.power_mod(). The documentation (misleadingly) only talks about integers, but the implementation actually works for any ring where % is defined. I've created #33244 to change this.

This means all we really need to do for this ticket is adding code along the lines of

if modulus is not None:
    return power_mod(self, exponent, modulus)

to all the .__pow__() methods in polynomial-ring classes that do not yet handle the modulus argument.

[shashwat1002]

comment:22

Replying to @yyyyx4:

Replying to @shashwat1002:

• the polynomials are dealt by one of two C libraries: FLINT and NTL, so a general modular exponentiation function on polynomials must be written in python?

Yes. In fact, I discovered that we already have such a function in Sage: sage.arith.misc.power_mod(). The documentation (misleadingly) only talks about integers, but the implementation actually works for any ring where % is defined. I've created #33244 to change this.

This means all we really need to do for this ticket is adding code along the lines of

if modulus is not None:
    return power_mod(self, exponent, modulus)

to all the .__pow__() methods in polynomial-ring classes that do not yet handle the modulus argument.

Hi, I had already written binary exponentiation (modular) using FLINT C functions directly. I'll replace my work with this.

However, I am wondering, would there be a significant efficiency advantage writing it in cython (and making calls to FLINT functions as opposed to python as in power_mod)

[yyyyx4]

comment:23

There certainly are circumstances where saving on Python overhead is worth it, but my guess is that it won't give a noticeable speedup here. In this algorithm (except for toy-sized inputs), virtually all of the time will be spent on polynomial arithmetic, which is already implemented in low-level code even if you call it from Python.
But if your code is already functional, I figure there's no reason not to use it.

Out of curiousity: Do you have an application in mind where you need this operation to be really fast?

[sagetrac-git]

Changed commit from 8293154 to c37ccc8

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

c5054ad	Implement modular exponentiation in polynomial rings Flint.
c37ccc8	modular exponentiation using power_mod() in FLINT

[shashwat1002]

comment:25

Replying to @yyyyx4:

There certainly are circumstances where saving on Python overhead is worth it, but my guess is that it won't give a noticeable speedup here. In this algorithm (except for toy-sized inputs), virtually all of the time will be spent on polynomial arithmetic, which is already implemented in low-level code even if you call it from Python.
But if your code is already functional, I figure there's no reason not to use it.

Out of curiousity: Do you have an application in mind where you need this operation to be really fast?

Not at the moment, no.
I just did it with power_mod therefore now. Probably best to reuse functions.
If someone deems it necessary to make power_mod itself in cython then maybe I could use the code.

I couldn't find a __pow__ in the NTL implementation.

I am honestly struggling trying to find where to put appropriate code. Is there a detailed dev documentation somewhere?