diff --git a/src/flint/test/test.py b/src/flint/test/test.py index 7f6d5cd0..f41c9093 100644 --- a/src/flint/test/test.py +++ b/src/flint/test/test.py @@ -2090,6 +2090,18 @@ def test_fmpz_mod_poly(): assert f*f == f**2 assert f*f == f**fmpz(2) + # powmod + # 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 + 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.powmod(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) @@ -2121,6 +2133,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]) @@ -2606,8 +2625,9 @@ def setbad(obj, i, val): 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: + # TODO: this now fails as fmpz_mod_poly allows modulus + # assert raises(lambda: pow(P([1, 1]), 2, 3), NotImplementedError) assert P([1, 2, 1]).gcd(P([1, 1])) == P([1, 1]) assert raises(lambda: P([1, 2, 1]).gcd(None), TypeError) diff --git a/src/flint/types/fmpz_mod_poly.pyx b/src/flint/types/fmpz_mod_poly.pyx index f5b5d60b..5b1d8a3d 100644 --- a/src/flint/types/fmpz_mod_poly.pyx +++ b/src/flint/types/fmpz_mod_poly.pyx @@ -538,7 +538,7 @@ cdef class fmpz_mod_poly(flint_poly): def __pow__(self, e, mod=None): if mod is not None: - raise NotImplementedError + return self.powmod(e, mod) cdef fmpz_mod_poly res if e < 0: @@ -784,11 +784,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) @@ -800,12 +800,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): @@ -1110,10 +1143,14 @@ cdef class fmpz_mod_poly(flint_poly): ) return res - def powmod(self, e, modulus): + def powmod(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 ``powmod()`` can optimise + powering of polynomials with large exponents. >>> R = fmpz_mod_poly_ctx(163) >>> x = R.gen() @@ -1123,17 +1160,55 @@ cdef class fmpz_mod_poly(flint_poly): >>> >>> f.powmod(123, mod) 3*x^3 + 25*x^2 + 115*x + 161 + >>> f.powmod(2**64, mod) + 52*x^3 + 96*x^2 + 136*x + 9 + >>> mod_rev_inv = mod.reverse().inverse_series_trunc(4) + >>> f.powmod(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 powering 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 can use faster algorithms, first convert exp to fmpz type + e_fmpz = any_as_fmpz(e) + if e_fmpz is NotImplemented: + raise ValueError(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 windows 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):