Skip to content

Commit

Permalink
gh-36097: Add additional bindings from NTL to Polynomial_ZZ_pEX
Browse files Browse the repository at this point in the history
    
There are current cython bindings to NTL methods which are not made
available to polynomial rings over extension fields.

For some of these bindings, their introduction does not offer much speed
up, but there is a 50x improvement for using NTL reverse over that of
the current implementation, and about a 20% speed up using the NTL
`ZZ_pEX_InvTrunc` over what is currently implemented.

Method names and functionalities have been preserved, so this simply
improves performance for the `Polynomial_ZZ_pEX` type without effecting
other polynomial ring classes.

#### Before Patch

```py
sage: p = next_prime(2**1024)
sage: F.<i> = GF(p^2)
sage: R.<x> = PolynomialRing(F, implementation="NTL")
sage:
sage: f = R([p for p in primes(10)])
sage: %timeit f.reverse()
125 µs ± 582 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
sage: %timeit f.inverse_series_trunc(500)
4.96 ms ± 81.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
sage:
sage: f = R([p for p in primes(100)])
sage: %timeit f.reverse()
681 µs ± 4.21 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit f.inverse_series_trunc(500)
7.1 ms ± 62.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
sage:
sage: f = R([p for p in primes(1000)])
sage: %timeit f.reverse()
4.57 ms ± 77 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
sage: %timeit f.inverse_series_trunc(500)
10.9 ms ± 82.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

#### After Patch

```py
sage: p = next_prime(2**1024)
sage: F.<i> = GF(p^2)
sage: R.<x> = PolynomialRing(F, implementation="NTL")
sage:
sage: f = R([p for p in primes(10)])
sage: %timeit f.reverse()
836 ns ± 9.91 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops
each)
sage: %timeit f.inverse_series_trunc(500)
1.84 ms ± 77.1 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops
each)
sage:
sage: f = R([p for p in primes(100)])
sage: %timeit f.reverse()
3.63 µs ± 103 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops
each)
sage: %timeit f.inverse_series_trunc(500)
6.66 ms ± 441 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

sage: f = R([p for p in primes(1000)])
sage: %timeit f.reverse()
26.8 µs ± 1.48 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops
each)
sage: %timeit f.inverse_series_trunc(500)
8.75 ms ± 432 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```


### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [ ] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

The only additional import is `sig_on()` and `sig_off()` into the
`Polynomial_ZZ_pEX.pyx`

### Warning!

This is my first attempt to contribute new functions to sagemath, and
it's also my first few weeks of playing with cython. I have done my best
to do what is needed, but maybe I have done something incorrect. Please
let me know if there's more I can do.

In particular, I know an alternative edit to the one I proposed is to
modify the child `Polynomial_template` class instead of the
`Polynomial_ZZ_pEX` class, but I know less about all the other parents
of this class and I was concerned I may make breaking changes somewhere
downstream.

Currently, this just adds the NTL methods when I know they are avaialble
and they simply do what is done currently, but faster.
    
URL: #36097
Reported by: Giacomo Pope
Reviewer(s): Frédéric Chapoton
  • Loading branch information
Release Manager committed Aug 26, 2023
2 parents 53e8cbb + 1308baa commit 48ef6df
Showing 1 changed file with 129 additions and 0 deletions.
129 changes: 129 additions & 0 deletions src/sage/rings/polynomial/polynomial_zz_pex.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,10 @@ AUTHOR:
- Yann Laigle-Chapuy (2010-01) initial implementation
- Lorenz Panny (2023-01): :meth:`minpoly_mod`
- Giacomo Pope (2023-08): :meth:`reverse`, :meth:`inverse_series_trunc`
"""
from cysignals.signals cimport sig_on, sig_off

from sage.libs.ntl.ntl_ZZ_pEContext cimport ntl_ZZ_pEContext_class
from sage.libs.ntl.ZZ_pE cimport ZZ_pE_to_ZZ_pX
from sage.libs.ntl.ZZ_pX cimport ZZ_pX_deg, ZZ_pX_coeff
Expand Down Expand Up @@ -483,3 +486,129 @@ cdef class Polynomial_ZZ_pEX(Polynomial_template):
x^4 + x^3 + x
"""
return self.shift(-n)

def reverse(self, degree=None):
r"""
Return the polynomial obtained by reversing the coefficients
of this polynomial. If degree is set then this function behaves
as if this polynomial has degree ``degree``.
EXAMPLES::
sage: R.<x> = GF(101^2)[]
sage: f = x^13 + 11*x^10 + 32*x^6 + 4
sage: f.reverse()
4*x^13 + 32*x^7 + 11*x^3 + 1
sage: f.reverse(degree=15)
4*x^15 + 32*x^9 + 11*x^5 + x^2
sage: f.reverse(degree=2)
4*x^2
TESTS::
sage: R.<x> = GF(163^2)[]
sage: f = R([p for p in primes(20)])
sage: f.reverse()
2*x^7 + 3*x^6 + 5*x^5 + 7*x^4 + 11*x^3 + 13*x^2 + 17*x + 19
sage: f.reverse(degree=200)
2*x^200 + 3*x^199 + 5*x^198 + 7*x^197 + 11*x^196 + 13*x^195 + 17*x^194 + 19*x^193
sage: f.reverse(degree=0)
Traceback (most recent call last):
...
ValueError: degree argument must be a non-negative integer, got 0
sage: f.reverse(degree=-5)
Traceback (most recent call last):
...
ValueError: degree argument must be a non-negative integer, got -5
"""
self._parent._modulus.restore()

# Construct output polynomial
cdef Polynomial_ZZ_pEX r
r = Polynomial_ZZ_pEX.__new__(Polynomial_ZZ_pEX)
celement_construct(&r.x, (<Polynomial_template>self)._cparent)
r._parent = (<Polynomial_template>self)._parent
r._cparent = (<Polynomial_template>self)._cparent

# When a degree has been supplied, ensure it is a valid input
cdef unsigned long d
if degree is not None:
if degree <= 0:
raise ValueError("degree argument must be a non-negative integer, got %s" % (degree))
try:
d = degree
except ValueError:
raise ValueError("degree argument must be a non-negative integer, got %s" % (degree))
ZZ_pEX_reverse_hi(r.x, (<Polynomial_ZZ_pEX> self).x, d)
else:
ZZ_pEX_reverse(r.x, (<Polynomial_ZZ_pEX> self).x)
return r

def inverse_series_trunc(self, prec):
r"""
Compute and return the inverse of self modulo `x^prec`.
The constant term of self must be invertible.
EXAMPLES::
sage: R.<x> = GF(101^2)[]
sage: z2 = R.base_ring().gen()
sage: f = (3*z2 + 57)*x^3 + (13*z2 + 94)*x^2 + (7*z2 + 2)*x + 66*z2 + 15
sage: f.inverse_series_trunc(1)
51*z2 + 92
sage: f.inverse_series_trunc(2)
(30*z2 + 30)*x + 51*z2 + 92
sage: f.inverse_series_trunc(3)
(42*z2 + 94)*x^2 + (30*z2 + 30)*x + 51*z2 + 92
sage: f.inverse_series_trunc(4)
(99*z2 + 96)*x^3 + (42*z2 + 94)*x^2 + (30*z2 + 30)*x + 51*z2 + 92
TESTS::
sage: R.<x> = GF(163^2)[]
sage: f = R([p for p in primes(20)])
sage: f.inverse_series_trunc(1)
82
sage: f.inverse_series_trunc(2)
40*x + 82
sage: f.inverse_series_trunc(3)
61*x^2 + 40*x + 82
sage: f.inverse_series_trunc(0)
Traceback (most recent call last):
...
ValueError: the precision must be positive, got 0
sage: f.inverse_series_trunc(-1)
Traceback (most recent call last):
...
ValueError: the precision must be positive, got -1
sage: f = x + x^2 + x^3
sage: f.inverse_series_trunc(5)
Traceback (most recent call last):
...
ValueError: constant term 0 is not a unit
"""
self._parent._modulus.restore()

# Ensure precision is non-negative
if prec <= 0:
raise ValueError("the precision must be positive, got {}".format(prec))

# Ensure we can invert the constant term
const_term = self.get_coeff_c(0)
if not const_term.is_unit():
raise ValueError("constant term {} is not a unit".format(const_term))

# Construct output polynomial
cdef Polynomial_ZZ_pEX r
r = Polynomial_ZZ_pEX.__new__(Polynomial_ZZ_pEX)
celement_construct(&r.x, (<Polynomial_template>self)._cparent)
r._parent = (<Polynomial_template>self)._parent
r._cparent = (<Polynomial_template>self)._cparent

# Call to NTL for the inverse truncation
if prec > 0:
sig_on()
ZZ_pEX_InvTrunc(r.x, self.x, prec)
sig_off()
return r

0 comments on commit 48ef6df

Please sign in to comment.