From 99b91fbdd6dbb06e1bed0d83a7e796fc4cdcb740 Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Mon, 5 Aug 2024 10:40:54 +0100 Subject: [PATCH 1/7] add compose_mod and powmod with large exp --- src/flint/test/{test.py => test_all.py} | 538 ++++++++++++++++++++++-- src/flint/types/fmpz_mod_poly.pyx | 95 ++++- 2 files changed, 588 insertions(+), 45 deletions(-) rename src/flint/test/{test.py => test_all.py} (85%) diff --git a/src/flint/test/test.py b/src/flint/test/test_all.py similarity index 85% rename from src/flint/test/test.py rename to src/flint/test/test_all.py index 7f6d5cd0..3a78e2ac 100644 --- a/src/flint/test/test.py +++ b/src/flint/test/test_all.py @@ -5,13 +5,10 @@ import doctest import platform -from flint.utils.flint_exceptions import DomainError +from flint.utils.flint_exceptions import DomainError, IncompatibleContextError import flint -if sys.version_info[0] >= 3: - long = int - PYPY = platform.python_implementation() == "PyPy" ctx = flint.ctx @@ -35,7 +32,7 @@ def raises(f, exception): def test_pyflint(): - assert flint.__version__ == "0.6.0" + assert flint.__version__ == "0.7.0a4" ctx = flint.ctx assert str(ctx) == repr(ctx) == _default_ctx_string @@ -95,8 +92,8 @@ def test_fmpz(): assert raises(lambda: flint.fmpz([]), TypeError) for s in L: for t in L: - for ltype in (flint.fmpz, int, long): - for rtype in (flint.fmpz, int, long): + for ltype in (flint.fmpz, int): + for rtype in (flint.fmpz, int): assert (ltype(s) == rtype(t)) == (s == t) assert (ltype(s) != rtype(t)) == (s != t) @@ -297,6 +294,8 @@ def test_fmpz_functions(): [0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0]), (lambda n: flint.fmpz(n).is_perfect_power(), [T, T, T, F, F, T, F, F, F, T, T, F]), + (lambda n: flint.fmpz(n).is_square(), + [F, T, T, F, F, T, F, F, F, F, T, F]), (lambda n: flint.fmpz(n).partitions_p(), [0, 1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42]), (lambda n: flint.fmpz(n).moebius_mu(), @@ -370,7 +369,7 @@ def test_fmpz_poly(): assert raises(lambda: Z({}), TypeError) # XXX: This should probably be made to work: assert raises(lambda: Z((1,2,3)), TypeError) - for ztype in [int, long, flint.fmpz]: + for ztype in [int, flint.fmpz]: assert Z([1,2,3]) + ztype(5) == Z([6,2,3]) assert ztype(5) + Z([1,2,3]) == Z([6,2,3]) assert Z([1,2,3]) - ztype(5) == Z([-4,2,3]) @@ -526,8 +525,6 @@ def test_fmpz_mat(): assert raises(lambda: a + c, ValueError) assert (a * 3).entries() == [3,6,9,12,15,18] assert (3 * a).entries() == [3,6,9,12,15,18] - assert (a * long(3)).entries() == [3,6,9,12,15,18] - assert (long(3) * a).entries() == [3,6,9,12,15,18] assert (a * flint.fmpz(3)).entries() == [3,6,9,12,15,18] assert (flint.fmpz(3) * a).entries() == [3,6,9,12,15,18] assert M.randrank(5,7,3,10).rank() == 3 @@ -654,6 +651,7 @@ def set_bad(i,j): assert M6.minpoly() == flint.fmpz_poly([4,-4,1]) assert list(M6) == [2,0,0,0,2,1,0,0,2] + def test_fmpz_series(): Zp = flint.fmpz_poly Z = flint.fmpz_series @@ -756,7 +754,7 @@ def test_fmpq(): assert 1 != Q(2) assert Q(1) != () assert Q(1,2) != 1 - assert Q(2,3) == Q(flint.fmpz(2),long(3)) + assert Q(2,3) == Q(flint.fmpz(2),3) assert Q(-2,-4) == Q(1,2) assert Q("1") == Q(1) assert Q("1/2") == Q(1,2) @@ -1133,6 +1131,7 @@ def set_bad(i): assert M3.charpoly() == flint.fmpq_poly([-1,6,-12,8]) / 8 assert M3.minpoly() == flint.fmpq_poly([1,-4,4]) / 4 + def test_fmpq_series(): Qp = flint.fmpq_poly Q = flint.fmpq_series @@ -1481,8 +1480,6 @@ def test_nmod_mat(): assert raises(lambda: a + c, ValueError) assert (a * 3).entries() == [G(x,17) for x in [3,6,9,12,15,18]] assert (3 * a).entries() == [G(x,17) for x in [3,6,9,12,15,18]] - assert (a * long(3)).entries() == [G(x,17) for x in [3,6,9,12,15,18]] - assert (long(3) * a).entries() == [G(x,17) for x in [3,6,9,12,15,18]] assert (a * flint.fmpz(3)).entries() == [G(x,17) for x in [3,6,9,12,15,18]] assert (flint.fmpz(3) * a).entries() == [G(x,17) for x in [3,6,9,12,15,18]] assert M(2,2,[1,1,2,2],17).rank() == 1 @@ -1609,7 +1606,7 @@ def test_pickling(): def test_fmpz_mod(): from flint import fmpz_mod_ctx, fmpz, fmpz_mod - + p_sml = 163 p_med = 2**127 - 1 p_big = 2**255 - 19 @@ -1759,7 +1756,7 @@ def test_fmpz_mod(): assert raises(lambda: F_test(test_x) * "AAA", TypeError) assert raises(lambda: F_test(test_x) * F_other(test_x), ValueError) - # Exponentiation + # Exponentiation assert F_test(0)**0 == pow(0, 0, test_mod) assert F_test(0)**1 == pow(0, 1, test_mod) @@ -1809,7 +1806,7 @@ def test_fmpz_mod(): assert fmpz(test_y) / F_test(test_x) == (test_y * pow(test_x, -1, test_mod)) % test_mod assert test_y / F_test(test_x) == (test_y * pow(test_x, -1, test_mod)) % test_mod - + def test_fmpz_mod_dlog(): from flint import fmpz, fmpz_mod_ctx @@ -1831,7 +1828,7 @@ def test_fmpz_mod_dlog(): F = fmpz_mod_ctx(163) g = F(2) a = g**123 - + assert 123 == g.discrete_log(a) a_int = pow(2, 123, 163) @@ -1882,7 +1879,7 @@ def test_fmpz_mod_poly(): assert repr(R3) == "fmpz_mod_poly_ctx(13)" assert R1.modulus() == 11 - + assert R1.is_prime() assert R1.zero() == 0 assert R1.one() == 1 @@ -1951,7 +1948,7 @@ def test_fmpz_mod_poly(): assert str(f) == "8*x^3 + 7*x^2 + 6*x + 7" # TODO: currently repr does pretty printing - # just like str, we should address this. Mainly, + # just like str, we should address this. Mainly, # the issue is we want nice `repr` behaviour in # interactive shells, which currently is why this # choice has been made @@ -2008,14 +2005,14 @@ def test_fmpz_mod_poly(): f_bad = R_cmp([2,2,2,2,2]) for (F_test, R_test) in [(F_sml, R_sml), (F_med, R_med), (F_big, R_big)]: - + f = R_test([-1,-2]) g = R_test([-3,-4]) # pos, neg assert f is +f assert -f == R_test([1,2]) - + # add assert raises(lambda: f + f_cmp, ValueError) assert raises(lambda: f + "AAA", TypeError) @@ -2068,7 +2065,7 @@ def test_fmpz_mod_poly(): assert raises(lambda: f / "AAA", TypeError) assert raises(lambda: f / 0, ZeroDivisionError) assert raises(lambda: f_cmp / 2, ZeroDivisionError) - + assert (f + f) / 2 == f assert (f + f) / fmpz(2) == f assert (f + f) / F_test(2) == f @@ -2082,7 +2079,7 @@ def test_fmpz_mod_poly(): assert (f + f) // 2 == f assert (f + f) // fmpz(2) == f assert (f + f) // F_test(2) == f - assert 2 // R_test(2) == 1 + assert 2 // R_test(2) == 1 assert (f + 1) // f == 1 # pow @@ -2090,6 +2087,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 +2130,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]) @@ -2176,7 +2192,7 @@ def test_fmpz_mod_poly(): f1 = R_test([-3, 1]) f2 = R_test([-5, 1]) assert f1.resultant(f2) == (3 - 5) - assert raises(lambda: f.resultant("AAA"), TypeError) + assert raises(lambda: f.resultant("AAA"), TypeError) # sqrt f1 = R_test.random_element(irreducible=True) @@ -2433,14 +2449,14 @@ def _all_polys(): (flint.fmpz_poly, flint.fmpz, False), (flint.fmpq_poly, flint.fmpq, True), (lambda *a: flint.nmod_poly(*a, 17), lambda x: flint.nmod(x, 17), True), - (lambda *a: flint.fmpz_mod_poly(*a, flint.fmpz_mod_poly_ctx(163)), - lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(163)), + (lambda *a: flint.fmpz_mod_poly(*a, flint.fmpz_mod_poly_ctx(163)), + lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(163)), True), - (lambda *a: flint.fmpz_mod_poly(*a, flint.fmpz_mod_poly_ctx(2**127 - 1)), - lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(2**127 - 1)), + (lambda *a: flint.fmpz_mod_poly(*a, flint.fmpz_mod_poly_ctx(2**127 - 1)), + lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(2**127 - 1)), True), - (lambda *a: flint.fmpz_mod_poly(*a, flint.fmpz_mod_poly_ctx(2**255 - 19)), - lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(2**255 - 19)), + (lambda *a: flint.fmpz_mod_poly(*a, flint.fmpz_mod_poly_ctx(2**255 - 19)), + lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(2**255 - 19)), True), ] @@ -2472,6 +2488,28 @@ def test_polys(): assert (P([1]) == P([2])) is False assert (P([1]) != P([2])) is True + assert (P([1]) == 1) is True + assert (P([1]) != 1) is False + assert (P([1]) == 2) is False + assert (P([1]) != 2) is True + + assert (1 == P([1])) is True + assert (1 != P([1])) is False + assert (2 == P([1])) is False + assert (2 != P([1])) is True + + s1, s2 = S(1), S(2) + + assert (P([s1]) == s1) is True + assert (P([s1]) != s1) is False + assert (P([s1]) == s2) is False + assert (P([s1]) != s2) is True + + assert (s1 == P([s1])) is True + assert (s1 != P([s1])) is False + assert (s1 == P([s2])) is False + assert (s1 != P([s2])) is True + assert (P([1]) == None) is False assert (P([1]) != None) is True assert (None == P([1])) is False @@ -2505,12 +2543,17 @@ def setbad(obj, i, val): assert raises(lambda: setbad(p, -1, 1), ValueError) for v in [], [1], [1, 2]: - if P == flint.fmpz_poly: + p = P(v) + if type(p) == flint.fmpz_poly: assert P(v).repr() == f'fmpz_poly({v!r})' - elif P == flint.fmpq_poly: + elif type(p) == flint.fmpq_poly: assert P(v).repr() == f'fmpq_poly({v!r})' - elif P == flint.nmod_poly: + elif type(p) == flint.nmod_poly: assert P(v).repr() == f'nmod_poly({v!r}, 17)' + elif type(p) == flint.fmpz_mod_poly: + pass # fmpz_mod_poly does not have .repr() ... + else: + assert False assert repr(P([])) == '0' assert repr(P([1])) == '1' @@ -2526,6 +2569,12 @@ def setbad(obj, i, val): assert bool(P([])) is False assert bool(P([1])) is True + assert P([]).is_zero() is True + assert P([1]).is_zero() is False + + assert P([]).is_one() is False + assert P([1]).is_one() is True + assert +P([1, 2, 3]) == P([1, 2, 3]) assert -P([1, 2, 3]) == P([-1, -2, -3]) @@ -2606,8 +2655,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) @@ -2638,6 +2688,421 @@ def setbad(obj, i, val): 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), + (flint.fmpq_mpoly, flint.fmpq_mpoly_ctx, flint.fmpq, True), + ] + + +def test_mpolys(): + for P, C, S, is_field in _all_mpolys(): + + ctx = C.get_context(nvars=2) + + assert raises(lambda: C.get_context(nvars=2, ordering="bad"), TypeError) + assert raises(lambda: C.get_context(nvars=-1), ValueError) + assert raises(lambda: C(-1, flint.Ordering.lex, []), ValueError) + assert raises(lambda: ctx.constant("bad"), TypeError) + assert raises(lambda: ctx.from_dict("bad"), ValueError) + assert raises(lambda: ctx.from_dict({(0, 0): "bad"}), TypeError) + assert raises(lambda: ctx.from_dict({(0, "bad"): 1}), TypeError) + assert raises(lambda: ctx.from_dict({(0,): 1}), ValueError) + assert raises(lambda: ctx.gen(-1), IndexError) + assert raises(lambda: ctx.gen(10), IndexError) + + assert raises(lambda: P(val=C.get_context(nvars=1).constant(0), ctx=ctx), IncompatibleContextError) + assert raises(lambda: P(val={}, ctx=None), ValueError) + assert raises(lambda: P(val={"bad": 1}, ctx=None), ValueError) + assert raises(lambda: P(val="1", ctx=None), ValueError) + + assert P(val={(0, 0): 1}, ctx=ctx) == ctx.from_dict({(0, 0): 1}) + assert P(ctx=ctx).context() == ctx + assert P(1, ctx=ctx).is_one() + assert ctx.gen(1) == ctx.from_dict({(0, 1): 1}) + + def mpoly(x): + return ctx.from_dict(x) + + def quick_poly(): + return mpoly({(0, 0): 1, (0, 1): 2, (1, 0): 3, (2, 2): 4}) + + assert ctx.nvars() == 2 + assert ctx.ordering() == flint.Ordering.lex + + ctx1 = C.get_context(4) + assert [ctx1.name(i) for i in range(4)] == ['x0', 'x1', 'x2', 'x3'] + for order in list(flint.Ordering): + ctx1 = C.get_context(4, order) + assert ctx1.ordering() == order + + assert ctx.constant(1) == mpoly({(0, 0): 1}) == P(1, ctx=ctx) + + assert raises(lambda: P([None]), TypeError) + assert raises(lambda: P(object()), TypeError), f"{P(object()) = }" + assert raises(lambda: P(None), TypeError) + assert raises(lambda: P(None, None), TypeError) + assert raises(lambda: P([1,2], None), TypeError) + assert raises(lambda: P(1, None), ValueError) + + assert len(P(ctx=ctx)) == len(mpoly({(0, 0): 0})) == 0 + assert len(P(1, ctx=ctx)) == len(mpoly({(0, 0): 1})) == 1 + assert len(mpoly({(0, 0): 1, (0, 1): 1})) == 2 + assert len(mpoly({(0, 0): 1, (0, 1): 1, (1, 0): 1})) == 3 + + # degree is -1 when empty poly + assert P(ctx=ctx).degrees() == mpoly({(0, 0): 0}).degrees() == (-1, -1) + assert P(1, ctx=ctx).degrees() == mpoly({(0, 0): 1}).degrees() == (0, 0) + assert mpoly({(0, 0): 1, (0, 1): 1}).degrees() == (0, 1) + assert mpoly({(0, 0): 1, (0, 1): 1, (1, 0): 1}).degrees() == (1, 1) + assert mpoly({(0, 0): 1, (0, 1): 1, (1, 0): 1, (2, 2): 2}).degrees() == (2, 2) + + assert (P(1, ctx=ctx) == P(1, ctx=ctx)) is True + assert (P(1, ctx=ctx) != P(1, ctx=ctx)) is False + assert (P(1, ctx=ctx) == P(2, ctx=ctx)) is False + assert (P(1, ctx=ctx) != P(2, ctx=ctx)) is True + + assert (P(1, ctx=ctx) == 1) is False + assert (P(1, ctx=ctx) != 1) is True + assert (1 == P(1, ctx=ctx)) is False + assert (1 != P(1, ctx=ctx)) is True + + assert (P(1, ctx=ctx) == P(1, ctx=ctx1)) is False + assert (P(1, ctx=ctx) != P(1, ctx=ctx1)) is True + + assert (P(1, ctx=ctx) == None) is False + assert (P(1, ctx=ctx) != None) is True + assert (None == P(1, ctx=ctx)) is False + assert (None != P(1, ctx=ctx)) is True + + assert P(ctx.from_dict({(0, 1): 3})) == ctx.from_dict({(0, 1): 3}) + assert P({(0, 1): 3}, ctx=ctx) == ctx.from_dict({(0, 1): 3}) + + if P is flint.fmpq_mpoly: + ctx_z = flint.fmpz_mpoly_ctx.get_context(2) + assert quick_poly() == P(ctx_z.from_dict({(0, 0): 1, (0, 1): 2, (1, 0): 3, (2, 2): 4})) + assert P(ctx_z.from_dict({(0, 0): 1}), ctx=ctx) == P({(0, 0): 1}, ctx=ctx) + + assert raises(lambda: P(ctx=ctx) < P(ctx=ctx), TypeError) + assert raises(lambda: P(ctx=ctx) <= P(ctx=ctx), TypeError) + assert raises(lambda: P(ctx=ctx) > P(ctx=ctx), TypeError) + assert raises(lambda: P(ctx=ctx) >= P(ctx=ctx), TypeError) + assert raises(lambda: P(ctx=ctx) < None, TypeError) + assert raises(lambda: P(ctx=ctx) <= None, TypeError) + assert raises(lambda: P(ctx=ctx) > None, TypeError) + assert raises(lambda: P(ctx=ctx) >= None, TypeError) + assert raises(lambda: None < P(ctx=ctx), TypeError) + assert raises(lambda: None <= P(ctx=ctx), TypeError) + assert raises(lambda: None > P(ctx=ctx), TypeError) + assert raises(lambda: None >= P(ctx=ctx), TypeError) + + p = quick_poly() + assert p.coefficient(2) == S(2) + assert raises(lambda: p.coefficient(-1), IndexError) + assert raises(lambda: p.coefficient(10), IndexError) + + assert raises(lambda: p[-1], TypeError) + assert raises(lambda: p[4], TypeError) + + assert p[(2, 2)] == 4 + assert p[(0, 0)] == 1 + assert raises(lambda: p[(1,)], ValueError) + assert raises(lambda: p[(1, "bad")], TypeError) + assert raises(lambda: p["bad"], TypeError) + + p = quick_poly() + p[(1, 0)] = S(10) + assert p == mpoly({(0, 0): 1, (0, 1): 2, (1, 0): 10, (2, 2): 4}) + + p = quick_poly() + p[(1, 0)] = p[(1, 0)] + assert p == quick_poly() + assert (1, 0) in p + assert (100, 100) not in p + + assert raises(lambda: p.__setitem__((4,), 1), ValueError) + + assert raises(lambda: p.__setitem__((1,), 1), ValueError) + assert raises(lambda: p.__setitem__((1, "bad"), 1), TypeError) + assert raises(lambda: p.__setitem__(("bad", 1), 1), TypeError) + + assert raises(lambda: p.__setitem__((2, 1), None), TypeError) + + assert P(ctx=ctx).repr() == f"{ctx.__class__.__name__}(2, '', ('x0', 'x1')).from_dict({{}})" + assert P(1, ctx=ctx).repr() == f"{ctx.__class__.__name__}(2, '', ('x0', 'x1')).from_dict({{(0, 0): 1}})" + assert str(quick_poly()) == repr(quick_poly()) == '4*x0^2*x1^2 + 3*x0 + 2*x1 + 1' + + assert p.monomial(0) == (2, 2) + assert p.monomial(3) == (0, 0) + assert raises(lambda: p.monomial(-1), IndexError) + assert raises(lambda: p.monomial(4), IndexError) + + assert p.total_degree() == 4 + assert P(ctx=ctx).total_degree() == -1 + assert P(1, ctx=ctx).total_degree() == 0 + + p = quick_poly() + assert p(0, 0) == p(0, S(0)) == p(S(0), S(0)) == S(1) == 1 + assert p(1, 1) == S(10) == 10 + + p = quick_poly() + assert p.monoms() == [(2, 2), (1, 0), (0, 1), (0, 0)] + assert p.coeffs() == [4, 3, 2, 1] + assert list(p.terms()) == list(zip([(2, 2), (1, 0), (0, 1), (0, 0)], [4, 3, 2, 1])) + + assert p.subs({"x1": S(0), 0: S(0)}) == ctx.from_dict({(0, 0): 1}) + assert p.compose(p.subs({"x1": 0}), ctx.from_dict({(0, 1): 1})) == mpoly({ + (2, 2): 36, + (1, 2): 24, + (1, 0): 9, + (0, 2): 4, + (0, 1): 2, + (0, 0): 4 + }) + assert p.compose(ctx.from_dict({(1, 0): 1}), ctx.from_dict({(0, 1): 1})) == p + + assert raises(lambda: p(None, None), TypeError) + assert raises(lambda: p(1), ValueError) + assert raises(lambda: p(0, 1, 2), ValueError) + + assert raises(lambda: p.subs({"x0": None}), TypeError) + assert raises(lambda: p.subs({"x0": None, "x1": None}), TypeError) + assert raises(lambda: p.subs({"a": 1}), ValueError) + assert raises(lambda: p.subs({"x0": 0, "x1": 1, "x2": 2}), ValueError) + + no_gens_ctx = C.get_context(0) + no_gens_p = P("2", no_gens_ctx) + assert no_gens_p.compose(ctx=ctx1).context() is ctx1 + assert raises(lambda: no_gens_p.compose(), ValueError) + + assert raises(lambda: p.compose(p, P(ctx=ctx1)), IncompatibleContextError) + + assert bool(P(ctx=ctx)) is False + assert bool(P(1, ctx=ctx)) is True + + assert P(ctx=ctx).is_zero() is True + assert P(1, ctx=ctx).is_zero() is False + + assert P(ctx=ctx).is_one() is False + assert P(1, ctx=ctx).is_one() is True + + assert +quick_poly() \ + == quick_poly() + assert -quick_poly() \ + == mpoly({(0, 0): -1, (0, 1): -2, (1, 0): -3, (2, 2): -4}) + + assert quick_poly() \ + + mpoly({(0, 0): 5, (0, 1): 6, (1, 0): 7, (2, 2): 8}) \ + == mpoly({(0, 0): 6, (0, 1): 8, (1, 0): 10, (2, 2): 12}) + + for T in [int, S, flint.fmpz, lambda x: P(x, ctx=ctx)]: + p = quick_poly() + p += T(1) + q = quick_poly() + assert q.iadd(T(1)) is None + assert quick_poly() + T(1) \ + == p == q == mpoly({(0, 0): 2, (0, 1): 2, (1, 0): 3, (2, 2): 4}) + assert T(1) + quick_poly() \ + == mpoly({(0, 0): 2, (0, 1): 2, (1, 0): 3, (2, 2): 4}) + + assert raises(lambda: mpoly({(0, 0): 2, (0, 1): 2, (1, 0): 3, (2, 2): 4}) + None, TypeError) + assert raises(lambda: None + mpoly({(0, 0): 2, (0, 1): 2, (1, 0): 3, (2, 2): 4}), TypeError) + assert raises(lambda: quick_poly() + P(ctx=ctx1), IncompatibleContextError) + assert raises(lambda: quick_poly().iadd(P(ctx=ctx1)), IncompatibleContextError) + assert raises(lambda: quick_poly().iadd(None), NotImplementedError) + + assert quick_poly() - mpoly({(0, 0): 5, (0, 1): 6, (1, 0): 7, (2, 2): 8}) \ + == mpoly({(0, 0): -4, (0, 1): -4, (1, 0): -4, (2, 2): -4}) + + for T in [int, S, flint.fmpz, lambda x: P(x, ctx=ctx)]: + p = quick_poly() + p -= T(1) + q = quick_poly() + assert q.isub(T(1)) is None + assert quick_poly() - T(1) == p == q == mpoly({(0, 1): 2, (1, 0): 3, (2, 2): 4}) + assert T(1) - quick_poly() == mpoly({(0, 1): -2, (1, 0): -3, (2, 2): -4}) + + assert raises(lambda: quick_poly() - None, TypeError) + assert raises(lambda: None - quick_poly(), TypeError) + assert raises(lambda: quick_poly() - P(ctx=ctx1), IncompatibleContextError) + assert raises(lambda: quick_poly().isub(P(ctx=ctx1)), IncompatibleContextError) + assert raises(lambda: quick_poly().isub(None), NotImplementedError) + + assert quick_poly() * mpoly({(1, 0): 5, (0, 1): 6}) \ + == mpoly({ + (3, 2): 20, + (2, 3): 24, + (2, 0): 15, + (1, 1): 28, + (1, 0): 5, + (0, 2): 12, + (0, 1): 6 + }) + + for T in [int, S, flint.fmpz, lambda x: P(x, ctx=ctx)]: + p = quick_poly() + p *= T(2) + q = quick_poly() + assert q.imul(T(2)) is None + assert quick_poly() * T(2) == p == q == mpoly({(0, 0): 2, (0, 1): 4, (1, 0): 6, (2, 2): 8}) + assert T(2) * quick_poly() == mpoly({(0, 0): 2, (0, 1): 4, (1, 0): 6, (2, 2): 8}) + + assert raises(lambda: quick_poly() * None, TypeError) + assert raises(lambda: None * quick_poly(), TypeError) + assert raises(lambda: quick_poly() * P(ctx=ctx1), IncompatibleContextError) + assert raises(lambda: quick_poly().imul(P(ctx=ctx1)), IncompatibleContextError) + assert raises(lambda: quick_poly().imul(None), NotImplementedError) + + assert quick_poly() // mpoly({(1, 1): 1}) == mpoly({(1, 1): 4}) + assert quick_poly() % mpoly({(1, 1): 1}) \ + == mpoly({(1, 0): 3, (0, 1): 2, (0, 0): 1}) + assert divmod(quick_poly(), mpoly({(1, 1): 1})) \ + == (mpoly({(1, 1): 4}), mpoly({(1, 0): 3, (0, 1): 2, (0, 0): 1})) + + assert 1 / P(1, ctx=ctx) == P(1, ctx=ctx) + assert quick_poly() / 1 == quick_poly() + assert quick_poly() // 1 == quick_poly() + assert quick_poly() % 1 == P(ctx=ctx) + assert divmod(quick_poly(), 1) == (quick_poly(), P(ctx=ctx)) + + if is_field: + assert quick_poly() / 3 == mpoly({(0, 0): S(1, 3), (0, 1): S(2, 3), (1, 0): S(1), (2, 2): S(4, 3)}) + else: + assert raises(lambda: quick_poly() / 3, DomainError) + + assert 1 // quick_poly() == P(ctx=ctx) + assert 1 % quick_poly() == P(1, ctx=ctx) + assert divmod(1, quick_poly()) == (P(ctx=ctx), P(1, ctx=ctx)) + + assert raises(lambda: quick_poly() / None, TypeError) + assert raises(lambda: quick_poly() // None, TypeError) + assert raises(lambda: quick_poly() % None, TypeError) + assert raises(lambda: divmod(quick_poly(), None), TypeError) + + assert raises(lambda: None / quick_poly(), TypeError) + assert raises(lambda: None // quick_poly(), TypeError) + assert raises(lambda: None % quick_poly(), TypeError) + assert raises(lambda: divmod(None, quick_poly()), TypeError) + + assert raises(lambda: quick_poly() / 0, ZeroDivisionError) + assert raises(lambda: quick_poly() // 0, ZeroDivisionError) + assert raises(lambda: quick_poly() % 0, ZeroDivisionError) + assert raises(lambda: divmod(quick_poly(), 0), ZeroDivisionError) + + assert raises(lambda: 1 / P(ctx=ctx), ZeroDivisionError) + assert raises(lambda: 1 // P(ctx=ctx), ZeroDivisionError) + assert raises(lambda: 1 % P(ctx=ctx), ZeroDivisionError) + assert raises(lambda: divmod(1, P(ctx=ctx)), ZeroDivisionError) + + assert raises(lambda: quick_poly() / P(ctx=ctx), ZeroDivisionError) + assert raises(lambda: quick_poly() // P(ctx=ctx), ZeroDivisionError) + assert raises(lambda: quick_poly() % P(ctx=ctx), ZeroDivisionError) + assert raises(lambda: divmod(quick_poly(), P(ctx=ctx)), ZeroDivisionError) + + assert raises(lambda: quick_poly() / P(1, ctx=ctx1), IncompatibleContextError) + assert raises(lambda: quick_poly() // P(1, ctx=ctx1), IncompatibleContextError) + assert raises(lambda: quick_poly() % P(1, ctx=ctx1), IncompatibleContextError) + assert raises(lambda: divmod(quick_poly(), P(1, ctx=ctx1)), IncompatibleContextError) + + f = mpoly({(1, 1): 4, (0, 0): 1}) + g = mpoly({(0, 1): 2, (1, 0): 2}) + assert f * g / mpoly({(0, 1): 1, (1, 0): 1}) \ + == mpoly({(1, 1): 8, (0, 0): 2}) + + if not is_field: + assert raises(lambda: 1 / quick_poly(), DomainError) + assert raises(lambda: quick_poly() / P(2, ctx=ctx), DomainError) + + assert quick_poly() ** 0 == P(1, ctx=ctx) + assert quick_poly() ** 1 == quick_poly() + assert quick_poly() ** 2 == mpoly({ + (4, 4): 16, + (3, 2): 24, + (2, 3): 16, + (2, 2): 8, + (2, 0): 9, + (1, 1): 12, + (1, 0): 6, + (0, 2): 4, + (0, 1): 4, + (0, 0): 1, + }) + assert raises(lambda: P(ctx=ctx) ** -1, ValueError) + assert raises(lambda: P(ctx=ctx) ** None, TypeError) + + # # XXX: Not sure what this should do in general: + assert raises(lambda: pow(P(1, ctx=ctx), 2, 3), NotImplementedError) + + if is_field: + assert (f * g).gcd(f) == f / 4 + else: + assert (f * g).gcd(f) == f + assert raises(lambda: quick_poly().gcd(None), TypeError) + assert raises(lambda: quick_poly().gcd(P(ctx=ctx1)), IncompatibleContextError) + + assert (f * g).factor() == (S(2), [(mpoly({(0, 1): 1, (1, 0): 1}), 1), (f, 1)]) + + assert (f * f).sqrt() == f + if P is flint.fmpz_mpoly: + assert (f * f).sqrt(assume_perfect_square=True) == f + assert raises(lambda: quick_poly().sqrt(), ValueError) + + p = quick_poly() + assert p.derivative(0) == p.derivative("x0") == mpoly({(0, 0): 3, (1, 2): 8}) + assert p.derivative(1) == p.derivative("x1") == mpoly({(0, 0): 2, (2, 1): 8}) + + assert raises(lambda: p.derivative("x3"), ValueError) + assert raises(lambda: p.derivative(3), IndexError) + assert raises(lambda: p.derivative(None), TypeError) + + if is_field: + assert quick_poly().integral(0) == quick_poly().integral("x0") == \ + mpoly({(3, 2): S(4, 3), (2, 0): S(3, 2), (1, 1): S(2), (1, 0): S(1)}) + assert quick_poly().integral(1) == quick_poly().integral("x1") == \ + mpoly({(2, 3): S(4, 3), (1, 1): S(3), (0, 2): S(1), (0, 1): S(1)}) + else: + assert quick_poly().integral(0) == quick_poly().integral("x0") == \ + (6, mpoly({(3, 2): 8, (2, 0): 9, (1, 1): 12, (1, 0): 6})) + assert quick_poly().integral(1) == quick_poly().integral("x1") == \ + (3, mpoly({(2, 3): 4, (1, 1): 9, (0, 2): 3, (0, 1): 3})) + + assert raises(lambda: p.integral("x3"), ValueError) + assert raises(lambda: p.integral(3), IndexError) + assert raises(lambda: p.integral(None), TypeError) + + +def _all_mpoly_vecs(): + return [(flint.fmpz_mpoly_ctx, flint.fmpz_mpoly_vec), (flint.fmpq_mpoly_ctx, flint.fmpq_mpoly_vec)] + + +def test_fmpz_mpoly_vec(): + for context, mpoly_vec in _all_mpoly_vecs(): + ctx = context.get_context(nvars=2) + ctx1 = context.get_context(nvars=4) + + vec = mpoly_vec(3, ctx) + + assert vec[0] == ctx.from_dict({}) + assert vec[1] == ctx.from_dict({}) + + assert vec[2] == ctx.from_dict({}) + + assert len(vec) == 3 + assert str(vec) == f"[{', '.join(str(ctx.from_dict({})) for _ in range(3))}]" + assert repr(vec) == f"{mpoly_vec.__name__}([{', '.join(str(ctx.from_dict({})) for _ in range(3))}], ctx={str(ctx)})" + + assert raises(lambda: vec[None], TypeError) + assert raises(lambda: vec[-1], IndexError) + + vec[1] = ctx.from_dict({(1, 1): 1}) + assert vec.to_tuple() == mpoly_vec([ctx.from_dict({}), ctx.from_dict({(1, 1): 1}), ctx.from_dict({})], ctx).to_tuple() + assert raises(lambda: vec.__setitem__(None, 0), TypeError) + assert raises(lambda: vec.__setitem__(-1, 0), IndexError) + assert raises(lambda: vec.__setitem__(0, 0), TypeError) + assert raises(lambda: vec.__setitem__(0, ctx1.from_dict({})), IncompatibleContextError) + + def _all_matrices(): """Return a list of matrix types and scalar types.""" R163 = flint.fmpz_mod_ctx(163) @@ -3181,6 +3646,9 @@ def test_all_tests(): test_division_matrix, test_polys, + test_mpolys, + + test_fmpz_mpoly_vec, test_matrices_eq, test_matrices_constructor, diff --git a/src/flint/types/fmpz_mod_poly.pyx b/src/flint/types/fmpz_mod_poly.pyx index f5b5d60b..d914c807 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 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 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 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): From bc1e1103ea4ce02877f0e5ccfe19c74c5db74b92 Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Mon, 5 Aug 2024 19:13:00 +0100 Subject: [PATCH 2/7] add the same methods for nmod_poly --- src/flint/flintlib/nmod_poly.pxd | 1 + src/flint/test/test_all.py | 10 +- src/flint/types/fmpz_mod_poly.pyx | 2 +- src/flint/types/nmod_poly.pyx | 186 +++++++++++++++++++++++++++++- 4 files changed, 196 insertions(+), 3 deletions(-) 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 71644721..d8ff37b3 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) diff --git a/src/flint/types/fmpz_mod_poly.pyx b/src/flint/types/fmpz_mod_poly.pyx index c4a5d0d9..32fb9943 100644 --- a/src/flint/types/fmpz_mod_poly.pyx +++ b/src/flint/types/fmpz_mod_poly.pyx @@ -1182,7 +1182,7 @@ cdef class fmpz_mod_poly(flint_poly): # For larger exponents we need to cast e to an fmpz first e_fmpz = any_as_fmpz(e) if e_fmpz is NotImplemented: - raise ValueError(f"exponent cannot be cast to an fmpz type: {e = }") + 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: diff --git a/src/flint/types/nmod_poly.pyx b/src/flint/types/nmod_poly.pyx index ca5421dd..e1c26188 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,115 @@ 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 inverse_series_trunc(self, slong n): + """ + Returns the inverse of ``self`` modulo `x^n`. + + >>> 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 + """ + 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): @@ -358,7 +474,7 @@ cdef class nmod_poly(flint_poly): def __pow__(nmod_poly self, exp, mod): cdef nmod_poly res if mod is not None: - raise NotImplementedError("nmod_poly modular exponentiation") + return self.powmod(exp, mod) if exp < 0: raise ValueError("negative exponent") res = nmod_poly.__new__(nmod_poly) @@ -366,6 +482,74 @@ cdef class nmod_poly(flint_poly): nmod_poly_pow(res.val, self.val, exp) return res + def powmod(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 ``powmod()`` 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.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 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. From 9589530b3f6118fe53b13c21b7ccc30cc15cf17f Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Mon, 5 Aug 2024 19:18:43 +0100 Subject: [PATCH 3/7] make mod an optional arg --- src/flint/types/nmod_poly.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/flint/types/nmod_poly.pyx b/src/flint/types/nmod_poly.pyx index e1c26188..9cd82398 100644 --- a/src/flint/types/nmod_poly.pyx +++ b/src/flint/types/nmod_poly.pyx @@ -471,7 +471,7 @@ 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: return self.powmod(exp, mod) From ae31437cdec4667a12367977c173c60bbb7e09ec Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Tue, 6 Aug 2024 17:59:24 +0100 Subject: [PATCH 4/7] powmod to pow_mod --- src/flint/test/test_all.py | 12 ++++++------ src/flint/types/fmpz_mod_poly.pyx | 12 ++++++------ src/flint/types/nmod_poly.pyx | 12 ++++++------ 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py index d8ff37b3..e7c45fd8 100644 --- a/src/flint/test/test_all.py +++ b/src/flint/test/test_all.py @@ -2095,17 +2095,17 @@ def test_fmpz_mod_poly(): assert f*f == f**2 assert f*f == f**fmpz(2) - # powmod + # 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 powmod + # 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.powmod(2**32, g, mod_rev_inv="A"), 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) @@ -2162,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) diff --git a/src/flint/types/fmpz_mod_poly.pyx b/src/flint/types/fmpz_mod_poly.pyx index 32fb9943..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: - return self.powmod(e, mod) + return self.pow_mod(e, mod) cdef fmpz_mod_poly res if e < 0: @@ -1137,13 +1137,13 @@ cdef class fmpz_mod_poly(flint_poly): ) return res - def powmod(self, e, modulus, mod_rev_inv=None): + 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 ``powmod()`` can optimise + precomputing it and passing it to ``pow_mod()`` can optimise powering of polynomials with large exponents. >>> R = fmpz_mod_poly_ctx(163) @@ -1152,12 +1152,12 @@ 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.powmod(2**64, mod) + >>> 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.powmod(2**64, mod, mod_rev_inv) + >>> f.pow_mod(2**64, mod, mod_rev_inv) 52*x^3 + 96*x^2 + 136*x + 9 """ cdef fmpz_mod_poly res diff --git a/src/flint/types/nmod_poly.pyx b/src/flint/types/nmod_poly.pyx index 9cd82398..b7c5b53b 100644 --- a/src/flint/types/nmod_poly.pyx +++ b/src/flint/types/nmod_poly.pyx @@ -474,7 +474,7 @@ cdef class nmod_poly(flint_poly): def __pow__(nmod_poly self, exp, mod=None): cdef nmod_poly res if mod is not None: - return self.powmod(exp, mod) + return self.pow_mod(exp, mod) if exp < 0: raise ValueError("negative exponent") res = nmod_poly.__new__(nmod_poly) @@ -482,13 +482,13 @@ cdef class nmod_poly(flint_poly): nmod_poly_pow(res.val, self.val, exp) return res - def powmod(self, e, modulus, mod_rev_inv=None): + 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 ``powmod()`` can optimise + precomputing it and passing it to ``pow_mod()`` can optimise powering of polynomials with large exponents. >>> x = nmod_poly([0,1], 163) @@ -496,12 +496,12 @@ cdef class nmod_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.powmod(2**64, mod) + >>> 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.powmod(2**64, mod, mod_rev_inv) + >>> f.pow_mod(2**64, mod, mod_rev_inv) 52*x^3 + 96*x^2 + 136*x + 9 """ cdef nmod_poly res From e0312f883f5dba72c5140234ccb4e12765606311 Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Tue, 6 Aug 2024 18:08:59 +0100 Subject: [PATCH 5/7] add leading coefficent --- src/flint/types/nmod_poly.pyx | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/flint/types/nmod_poly.pyx b/src/flint/types/nmod_poly.pyx index b7c5b53b..86f039f9 100644 --- a/src/flint/types/nmod_poly.pyx +++ b/src/flint/types/nmod_poly.pyx @@ -233,9 +233,23 @@ cdef class nmod_poly(flint_poly): 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`. + 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) @@ -245,6 +259,9 @@ cdef class nmod_poly(flint_poly): >>> f.inverse_series_trunc(5) 45*x^4 + 23*x^3 + 159*x^2 + 151*x + 110 """ + if n <= 0: + raise ValueError("n must be positive") + cdef nmod_poly res res = nmod_poly.__new__(nmod_poly) nmod_poly_init_preinv(res.val, self.val.mod.n, self.val.mod.ninv) From 462318ea1920bb0d24263ccc0319874afef851d5 Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Tue, 6 Aug 2024 20:59:56 +0100 Subject: [PATCH 6/7] ensure self is not zero --- src/flint/types/nmod_poly.pyx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/flint/types/nmod_poly.pyx b/src/flint/types/nmod_poly.pyx index 86f039f9..a5b7fadd 100644 --- a/src/flint/types/nmod_poly.pyx +++ b/src/flint/types/nmod_poly.pyx @@ -260,7 +260,10 @@ cdef class nmod_poly(flint_poly): 45*x^4 + 23*x^3 + 159*x^2 + 151*x + 110 """ if n <= 0: - raise ValueError("n must be positive") + 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) From e3c9d030d5faab3778cc705584eaefb161d295c5 Mon Sep 17 00:00:00 2001 From: Giacomo Pope Date: Tue, 6 Aug 2024 21:01:50 +0100 Subject: [PATCH 7/7] modify test for generic poly --- src/flint/test/test_all.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py index e7c45fd8..e1e957bf 100644 --- a/src/flint/test/test_all.py +++ b/src/flint/test/test_all.py @@ -2664,8 +2664,12 @@ def setbad(obj, i, val): assert raises(lambda: P([1, 1]) ** None, TypeError) # 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) + 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)