diff --git a/src/flint/flintlib/nmod_poly.pxd b/src/flint/flintlib/nmod_poly.pxd index 53d6580c..3448cbd6 100644 --- a/src/flint/flintlib/nmod_poly.pxd +++ b/src/flint/flintlib/nmod_poly.pxd @@ -55,6 +55,7 @@ cdef extern from "flint/nmod_poly.h": int nmod_poly_equal_trunc(const nmod_poly_t poly1, const nmod_poly_t poly2, slong n) int nmod_poly_is_zero(const nmod_poly_t poly) int nmod_poly_is_one(const nmod_poly_t poly) + int nmod_poly_is_gen(const nmod_poly_t poly) void _nmod_poly_shift_left(mp_ptr res, mp_srcptr poly, slong len, slong k) void nmod_poly_shift_left(nmod_poly_t res, const nmod_poly_t poly, slong k) void _nmod_poly_shift_right(mp_ptr res, mp_srcptr poly, slong len, slong k) diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py index aca7668e..e1e957bf 100644 --- a/src/flint/test/test_all.py +++ b/src/flint/test/test_all.py @@ -1422,7 +1422,15 @@ def test_nmod_poly(): assert raises(lambda: [] * s, TypeError) assert raises(lambda: [] // s, TypeError) assert raises(lambda: [] % s, TypeError) - assert raises(lambda: pow(P([1,2],3), 3, 4), NotImplementedError) + assert raises(lambda: [] % s, TypeError) + assert raises(lambda: s.reverse(-1), ValueError) + assert raises(lambda: s.compose("A"), TypeError) + assert raises(lambda: s.compose_mod(s, "A"), TypeError) + assert raises(lambda: s.compose_mod("A", P([3,6,9],17)), TypeError) + assert raises(lambda: s.compose_mod(s, P([0], 17)), ZeroDivisionError) + assert raises(lambda: pow(s, -1, P([3,6,9],17)), ValueError) + assert raises(lambda: pow(s, 1, "A"), TypeError) + assert raises(lambda: pow(s, "A", P([3,6,9],17)), TypeError) assert str(P([1,2,3],17)) == "3*x^2 + 2*x + 1" assert P([1,2,3],17).repr() == "nmod_poly([1, 2, 3], 17)" p = P([3,4,5],17) @@ -2087,6 +2095,18 @@ def test_fmpz_mod_poly(): assert f*f == f**2 assert f*f == f**fmpz(2) + # pow_mod + # assert ui and fmpz exp agree for polynomials and generators + R_gen = R_test.gen() + assert pow(f, 2**60, g) == pow(pow(f, 2**30, g), 2**30, g) + assert pow(R_gen, 2**60, g) == pow(pow(R_gen, 2**30, g), 2**30, g) + + # Check other typechecks for pow_mod + assert raises(lambda: pow(f, -2, g), ValueError) + assert raises(lambda: pow(f, 1, "A"), TypeError) + assert raises(lambda: pow(f, "A", g), TypeError) + assert raises(lambda: f.pow_mod(2**32, g, mod_rev_inv="A"), TypeError) + # Shifts assert raises(lambda: R_test([1,2,3]).left_shift(-1), ValueError) assert raises(lambda: R_test([1,2,3]).right_shift(-1), ValueError) @@ -2118,6 +2138,13 @@ def test_fmpz_mod_poly(): # compose assert raises(lambda: h.compose("AAA"), TypeError) + # compose mod + mod = R_test([1,2,3,4]) + assert f.compose(h) % mod == f.compose_mod(h, mod) + assert raises(lambda: h.compose_mod("AAA", mod), TypeError) + assert raises(lambda: h.compose_mod(f, "AAA"), TypeError) + assert raises(lambda: h.compose_mod(f, R_test(0)), ZeroDivisionError) + # Reverse assert raises(lambda: h.reverse(degree=-100), ValueError) assert R_test([-1,-2,-3]).reverse() == R_test([-3,-2,-1]) @@ -2135,9 +2162,9 @@ def test_fmpz_mod_poly(): assert raises(lambda: f.mulmod(f, "AAA"), TypeError) assert raises(lambda: f.mulmod("AAA", g), TypeError) - # powmod - assert f.powmod(2, g) == (f*f) % g - assert raises(lambda: f.powmod(2, "AAA"), TypeError) + # pow_mod + assert f.pow_mod(2, g) == (f*f) % g + assert raises(lambda: f.pow_mod(2, "AAA"), TypeError) # divmod S, T = f.divmod(g) @@ -2635,9 +2662,14 @@ def setbad(obj, i, val): assert P([1, 1]) ** 2 == P([1, 2, 1]) assert raises(lambda: P([1, 1]) ** -1, ValueError) assert raises(lambda: P([1, 1]) ** None, TypeError) - - # # XXX: Not sure what this should do in general: - assert raises(lambda: pow(P([1, 1]), 2, 3), NotImplementedError) + + # XXX: Not sure what this should do in general: + p = P([1, 1]) + mod = P([1, 1]) + if type(p) not in [flint.fmpz_mod_poly, flint.nmod_poly]: + assert raises(lambda: pow(p, 2, mod), NotImplementedError) + else: + assert p * p % mod == pow(p, 2, mod) assert P([1, 2, 1]).gcd(P([1, 1])) == P([1, 1]) assert raises(lambda: P([1, 2, 1]).gcd(None), TypeError) @@ -2667,7 +2699,6 @@ def setbad(obj, i, val): if is_field: assert P([1, 2, 1]).integral() == P([0, 1, 1, S(1)/3]) - def _all_mpolys(): return [ (flint.fmpz_mpoly, flint.fmpz_mpoly_ctx, flint.fmpz, False), diff --git a/src/flint/types/fmpz_mod_poly.pyx b/src/flint/types/fmpz_mod_poly.pyx index fa334c63..8bac8772 100644 --- a/src/flint/types/fmpz_mod_poly.pyx +++ b/src/flint/types/fmpz_mod_poly.pyx @@ -536,7 +536,7 @@ cdef class fmpz_mod_poly(flint_poly): def __pow__(self, e, mod=None): if mod is not None: - raise NotImplementedError + return self.pow_mod(e, mod) cdef fmpz_mod_poly res if e < 0: @@ -778,11 +778,11 @@ cdef class fmpz_mod_poly(flint_poly): return evaluations - def compose(self, input): + def compose(self, other): """ Returns the composition of two polynomials - To be precise about the order of composition, given ``self``, and ``input`` + To be precise about the order of composition, given ``self``, and ``other`` by `f(x)`, `g(x)`, returns `f(g(x))`. >>> R = fmpz_mod_poly_ctx(163) @@ -794,12 +794,45 @@ cdef class fmpz_mod_poly(flint_poly): 9*x^4 + 12*x^3 + 10*x^2 + 4*x + 1 """ cdef fmpz_mod_poly res - val = self.ctx.any_as_fmpz_mod_poly(input) + val = self.ctx.any_as_fmpz_mod_poly(other) if val is NotImplemented: - raise TypeError(f"Cannot compose the polynomial with input: {input}") + raise TypeError(f"Cannot compose the polynomial with input: {other}") res = self.ctx.new_ctype_poly() fmpz_mod_poly_compose(res.val, self.val, (val).val, self.ctx.mod.val) + return res + + def compose_mod(self, other, modulus): + """ + Returns the composition of two polynomials modulo a third. + + To be precise about the order of composition, given ``self``, and ``other`` + and ``modulus`` by `f(x)`, `g(x)` and `h(x)`, returns `f(g(x)) \mod h(x)`. + We require that `h(x)` is non-zero. + + >>> R = fmpz_mod_poly_ctx(163) + >>> f = R([1,2,3,4,5]) + >>> g = R([3,2,1]) + >>> h = R([1,0,1,0,1]) + >>> f.compose_mod(g, h) + 63*x^3 + 100*x^2 + 17*x + 63 + >>> g.compose_mod(f, h) + 147*x^3 + 159*x^2 + 4*x + 7 + """ + cdef fmpz_mod_poly res + val = self.ctx.any_as_fmpz_mod_poly(other) + if val is NotImplemented: + raise TypeError(f"cannot compose the polynomial with input: {other}") + + h = self.ctx.any_as_fmpz_mod_poly(modulus) + if h is NotImplemented: + raise TypeError(f"cannot reduce the polynomial with input: {modulus}") + + if h.is_zero(): + raise ZeroDivisionError("cannot reduce modulo zero") + + res = self.ctx.new_ctype_poly() + fmpz_mod_poly_compose_mod(res.val, self.val, (val).val, (h).val, self.ctx.mod.val) return res cpdef long length(self): @@ -1104,10 +1137,14 @@ cdef class fmpz_mod_poly(flint_poly): ) return res - def powmod(self, e, modulus): + def pow_mod(self, e, modulus, mod_rev_inv=None): """ Returns ``self`` raised to the power ``e`` modulo ``modulus``: - :math:`f^e \mod g` + :math:`f^e \mod g`/ + + ``mod_rev_inv`` is the inverse of the reverse of the modulus, + precomputing it and passing it to ``pow_mod()`` can optimise + powering of polynomials with large exponents. >>> R = fmpz_mod_poly_ctx(163) >>> x = R.gen() @@ -1115,19 +1152,57 @@ cdef class fmpz_mod_poly(flint_poly): >>> g = 43*x**6 + 91*x**5 + 77*x**4 + 113*x**3 + 71*x**2 + 132*x + 60 >>> mod = x**4 + 93*x**3 + 78*x**2 + 72*x + 149 >>> - >>> f.powmod(123, mod) + >>> f.pow_mod(123, mod) 3*x^3 + 25*x^2 + 115*x + 161 + >>> f.pow_mod(2**64, mod) + 52*x^3 + 96*x^2 + 136*x + 9 + >>> mod_rev_inv = mod.reverse().inverse_series_trunc(4) + >>> f.pow_mod(2**64, mod, mod_rev_inv) + 52*x^3 + 96*x^2 + 136*x + 9 """ cdef fmpz_mod_poly res + if e < 0: + raise ValueError("Exponent must be non-negative") + modulus = self.ctx.any_as_fmpz_mod_poly(modulus) if modulus is NotImplemented: raise TypeError(f"Cannot interpret {modulus} as a polynomial") + # Output polynomial res = self.ctx.new_ctype_poly() - fmpz_mod_poly_powmod_ui_binexp( - res.val, self.val, e, (modulus).val, res.ctx.mod.val - ) + + # For small exponents, use a simple binary exponentiation method + if e.bit_length() < 32: + fmpz_mod_poly_powmod_ui_binexp( + res.val, self.val, e, (modulus).val, res.ctx.mod.val + ) + return res + + # For larger exponents we need to cast e to an fmpz first + e_fmpz = any_as_fmpz(e) + if e_fmpz is NotImplemented: + raise TypeError(f"exponent cannot be cast to an fmpz type: {e = }") + + # To optimise powering, we precompute the inverse of the reverse of the modulus + if mod_rev_inv is not None: + mod_rev_inv = self.ctx.any_as_fmpz_mod_poly(mod_rev_inv) + if mod_rev_inv is NotImplemented: + raise TypeError(f"Cannot interpret {mod_rev_inv} as a polynomial") + else: + mod_rev_inv = modulus.reverse().inverse_series_trunc(modulus.length()) + + # Use windowed exponentiation optimisation when self = x + if self.is_gen(): + fmpz_mod_poly_powmod_x_fmpz_preinv( + res.val, (e_fmpz).val, (modulus).val, (mod_rev_inv).val, res.ctx.mod.val + ) + return res + + # Otherwise using binary exponentiation for all other inputs + fmpz_mod_poly_powmod_fmpz_binexp_preinv( + res.val, self.val, (e_fmpz).val, (modulus).val, (mod_rev_inv).val, res.ctx.mod.val + ) return res def divmod(self, other): diff --git a/src/flint/types/nmod_poly.pyx b/src/flint/types/nmod_poly.pyx index ca5421dd..a5b7fadd 100644 --- a/src/flint/types/nmod_poly.pyx +++ b/src/flint/types/nmod_poly.pyx @@ -1,6 +1,7 @@ from cpython.list cimport PyList_GET_SIZE from flint.flint_base.flint_base cimport flint_poly from flint.utils.typecheck cimport typecheck +from flint.types.fmpz cimport fmpz, any_as_fmpz from flint.types.fmpz_poly cimport any_as_fmpz_poly from flint.types.fmpz_poly cimport fmpz_poly from flint.types.nmod cimport any_as_nmod @@ -182,6 +183,12 @@ cdef class nmod_poly(flint_poly): else: raise TypeError("cannot set element of type %s" % type(x)) + def degree(self): + return nmod_poly_degree(self.val) + + def length(self): + return nmod_poly_length(self.val) + def __bool__(self): return not nmod_poly_is_zero(self.val) @@ -191,6 +198,135 @@ cdef class nmod_poly(flint_poly): def is_one(self): return nmod_poly_is_one(self.val) + def is_gen(self): + return nmod_poly_is_gen(self.val) + + def reverse(self, degree=None): + """ + Return a polynomial with the coefficients of this polynomial + reversed. + + If ``degree`` is not None, the output polynomial will be zero-padded + or truncated before being reversed. NOTE: degree must be non-negative. + + >>> f = nmod_poly([1,2,3,4,5], 163) + >>> f.reverse() + x^4 + 2*x^3 + 3*x^2 + 4*x + 5 + >>> f.reverse(degree=1) + x + 2 + >>> f.reverse(degree=100) + x^100 + 2*x^99 + 3*x^98 + 4*x^97 + 5*x^96 + """ + cdef nmod_poly res + cdef slong d + + if degree is not None: + d = degree + if d != degree or d < 0: + raise ValueError(f"degree argument must be a non-negative integer, got {degree}") + length = d + 1 + else: + length = nmod_poly_length(self.val) + + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init_preinv(res.val, self.val.mod.n, self.val.mod.ninv) + nmod_poly_reverse(res.val, self.val, length) + return res + + def leading_coefficient(self): + """ + Return the leading coefficient of this polynomial. + + >>> f = nmod_poly([123, 129, 63, 14, 51, 76, 133], 163) + >>> f.leading_coefficient() + 133 + """ + d = self.degree() + if d < 0: + return 0 + return nmod_poly_get_coeff_ui(self.val, d) + + def inverse_series_trunc(self, slong n): + """ + Returns the inverse of ``self`` modulo `x^n`. Assumes the leading + coefficient of the polynomial is invertible. + + >>> f = nmod_poly([123, 129, 63, 14, 51, 76, 133], 163) + >>> f.inverse_series_trunc(3) + 159*x^2 + 151*x + 110 + >>> f.inverse_series_trunc(4) + 23*x^3 + 159*x^2 + 151*x + 110 + >>> f.inverse_series_trunc(5) + 45*x^4 + 23*x^3 + 159*x^2 + 151*x + 110 + """ + if n <= 0: + raise ValueError(f"{n = } must be positive") + + if self.is_zero(): + raise ValueError(f"cannot invert the zero element") + + cdef nmod_poly res + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init_preinv(res.val, self.val.mod.n, self.val.mod.ninv) + nmod_poly_inv_series(res.val, self.val, n) + return res + + def compose(self, other): + """ + Returns the composition of two polynomials + + To be precise about the order of composition, given ``self``, and ``other`` + by `f(x)`, `g(x)`, returns `f(g(x))`. + + >>> f = nmod_poly([1,2,3], 163) + >>> g = nmod_poly([0,0,1], 163) + >>> f.compose(g) + 3*x^4 + 2*x^2 + 1 + >>> g.compose(f) + 9*x^4 + 12*x^3 + 10*x^2 + 4*x + 1 + """ + cdef nmod_poly res + other = any_as_nmod_poly(other, (self).val.mod) + if other is NotImplemented: + raise TypeError("cannot convert input to nmod_poly") + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init_preinv(res.val, self.val.mod.n, self.val.mod.ninv) + nmod_poly_compose(res.val, self.val, (other).val) + return res + + def compose_mod(self, other, modulus): + """ + Returns the composition of two polynomials modulo a third. + + To be precise about the order of composition, given ``self``, and ``other`` + and ``modulus`` by `f(x)`, `g(x)` and `h(x)`, returns `f(g(x)) \mod h(x)`. + We require that `h(x)` is non-zero. + + >>> f = nmod_poly([1,2,3,4,5], 163) + >>> g = nmod_poly([3,2,1], 163) + >>> h = nmod_poly([1,0,1,0,1], 163) + >>> f.compose_mod(g, h) + 63*x^3 + 100*x^2 + 17*x + 63 + >>> g.compose_mod(f, h) + 147*x^3 + 159*x^2 + 4*x + 7 + """ + cdef nmod_poly res + g = any_as_nmod_poly(other, self.val.mod) + if g is NotImplemented: + raise TypeError(f"cannot convert {other = } to nmod_poly") + + h = any_as_nmod_poly(modulus, self.val.mod) + if h is NotImplemented: + raise TypeError(f"cannot convert {modulus = } to nmod_poly") + + if modulus.is_zero(): + raise ZeroDivisionError("cannot reduce modulo zero") + + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init_preinv(res.val, self.val.mod.n, self.val.mod.ninv) + nmod_poly_compose_mod(res.val, self.val, (other).val, (modulus).val) + return res + def __call__(self, other): cdef mp_limb_t c if any_as_nmod(&c, other, self.val.mod): @@ -355,10 +491,10 @@ cdef class nmod_poly(flint_poly): def __rmod__(s, t): return divmod(t, s)[1] # XXX - def __pow__(nmod_poly self, exp, mod): + def __pow__(nmod_poly self, exp, mod=None): cdef nmod_poly res if mod is not None: - raise NotImplementedError("nmod_poly modular exponentiation") + return self.pow_mod(exp, mod) if exp < 0: raise ValueError("negative exponent") res = nmod_poly.__new__(nmod_poly) @@ -366,6 +502,74 @@ cdef class nmod_poly(flint_poly): nmod_poly_pow(res.val, self.val, exp) return res + def pow_mod(self, e, modulus, mod_rev_inv=None): + """ + Returns ``self`` raised to the power ``e`` modulo ``modulus``: + :math:`f^e \mod g`/ + + ``mod_rev_inv`` is the inverse of the reverse of the modulus, + precomputing it and passing it to ``pow_mod()`` can optimise + powering of polynomials with large exponents. + + >>> x = nmod_poly([0,1], 163) + >>> f = 30*x**6 + 104*x**5 + 76*x**4 + 33*x**3 + 70*x**2 + 44*x + 65 + >>> g = 43*x**6 + 91*x**5 + 77*x**4 + 113*x**3 + 71*x**2 + 132*x + 60 + >>> mod = x**4 + 93*x**3 + 78*x**2 + 72*x + 149 + >>> + >>> f.pow_mod(123, mod) + 3*x^3 + 25*x^2 + 115*x + 161 + >>> f.pow_mod(2**64, mod) + 52*x^3 + 96*x^2 + 136*x + 9 + >>> mod_rev_inv = mod.reverse().inverse_series_trunc(4) + >>> f.pow_mod(2**64, mod, mod_rev_inv) + 52*x^3 + 96*x^2 + 136*x + 9 + """ + cdef nmod_poly res + + if e < 0: + raise ValueError("Exponent must be non-negative") + + modulus = any_as_nmod_poly(modulus, (self).val.mod) + if modulus is NotImplemented: + raise TypeError("cannot convert input to nmod_poly") + + # Output polynomial + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init_preinv(res.val, self.val.mod.n, self.val.mod.ninv) + + # For small exponents, use a simple binary exponentiation method + if e.bit_length() < 32: + nmod_poly_powmod_ui_binexp( + res.val, self.val, e, (modulus).val + ) + return res + + # For larger exponents we need to cast e to an fmpz first + e_fmpz = any_as_fmpz(e) + if e_fmpz is NotImplemented: + raise TypeError(f"exponent cannot be cast to an fmpz type: {e = }") + + # To optimise powering, we precompute the inverse of the reverse of the modulus + if mod_rev_inv is not None: + mod_rev_inv = any_as_nmod_poly(mod_rev_inv, (self).val.mod) + if mod_rev_inv is NotImplemented: + raise TypeError(f"Cannot interpret {mod_rev_inv} as a polynomial") + else: + mod_rev_inv = modulus.reverse().inverse_series_trunc(modulus.length()) + + # Use windowed exponentiation optimisation when self = x + if self.is_gen(): + nmod_poly_powmod_x_fmpz_preinv( + res.val, (e_fmpz).val, (modulus).val, (mod_rev_inv).val + ) + return res + + # Otherwise using binary exponentiation for all other inputs + nmod_poly_powmod_fmpz_binexp_preinv( + res.val, self.val, (e_fmpz).val, (modulus).val, (mod_rev_inv).val + ) + return res + def gcd(self, other): """ Returns the monic greatest common divisor of self and other.