diff --git a/src/flint/flint_base/flint_base.pyi b/src/flint/flint_base/flint_base.pyi
index ec2fbdf4..784cddba 100644
--- a/src/flint/flint_base/flint_base.pyi
+++ b/src/flint/flint_base/flint_base.pyi
@@ -1,7 +1,8 @@
 #
 #
 #
-from typing import Generic, TypeVar, Iterator, Iterable, Any, Final, Self, Mapping
+from typing import Generic, TypeVar, Iterator, Iterable, Any, Final, Self, Mapping, Sequence
+from abc import abstractmethod
 import enum
 
 
@@ -10,12 +11,14 @@ FLINT_RELEASE: Final[int]
 
 
 Telem = TypeVar('Telem', bound=flint_scalar)
-Telem_coerce = TypeVar('Telem_coerce')
+Telem_coerce = TypeVar('Telem_coerce', contravariant=True)
 Tmpoly = TypeVar('Tmpoly', bound=flint_mpoly)
 Tctx = TypeVar('Tctx', bound=flint_mpoly_context)
 
 Sctx = TypeVar('Sctx', bound=flint_mpoly_context)
 
+_str = str
+
 
 class flint_elem:
     pass
@@ -27,7 +30,7 @@ class flint_scalar(flint_elem):
 
 class flint_poly(flint_elem, Generic[Telem]):
     def __iter__(self) -> Iterator[Telem]: ...
-    def __getitem__(self, index: int) -> Telem: ...
+    def __getitem__(self, index: int, /) -> Telem: ...
     def coeffs(self) -> list[Telem]: ...
     def str(self, ascending: bool = False, var: str = "x", *args: Any, **kwargs: Any): ...
     def roots(self) -> list[tuple[Telem, int]]: ...
@@ -49,8 +52,8 @@ class flint_mpoly(flint_elem, Generic[Tctx, Telem, Telem_coerce]):
         ctx: Tctx | None = None
     ) -> None: ...
 
-    def str(self) -> str: ...
-    def repr(self) -> str: ...
+    def str(self) -> _str: ...
+    def repr(self) -> _str: ...
 
     def context(self) -> Tctx: ...
 
@@ -59,25 +62,22 @@ class flint_mpoly(flint_elem, Generic[Tctx, Telem, Telem_coerce]):
 
     def leading_coefficient(self) -> Telem: ...
     def to_dict(self) -> dict[tuple[int, ...], Telem]: ...
-    def terms(self) -> Iterable[tuple[tuple[int, ...], Telem]]: ...
 
     def is_one(self) -> bool: ...
     def is_zero(self) -> bool: ...
     def is_constant(self) -> bool: ...
 
     def __len__(self) -> int: ...
+    def __getitem__(self, index: tuple[int, ...]) -> Telem: ...
+    def __setitem__(self, index: tuple[int, ...], coeff: Telem | Telem_coerce | int) -> None: ...
+    def __iter__(self) -> Iterable[tuple[int, ...]]: ...
+    def __contains__(self, index: tuple[int, ...]) -> bool: ...
+
     def coefficient(self, i: int) -> Telem: ...
     def monomial(self, i: int) -> tuple[int, ...]: ...
-
+    def terms(self) -> Iterable[tuple[tuple[int, ...], Telem]]: ...
     def monoms(self) -> list[tuple[int, ...]]: ...
     def coeffs(self) -> list[Telem]: ...
-    def __getitem__(self, index: tuple[int, ...]) -> Telem: ...
-    def __setitem__(self, index: tuple[int, ...], coeff: Telem | Telem_coerce | int) -> None: ...
-
-    def subs(self, mapping: dict[str | int, Telem | Telem_coerce | int]) -> Self: ...
-    def compose(self, *args: Self, ctx: Tctx | None = None) -> Self: ...
-
-    def __call__(self, *args: Telem | Telem_coerce) -> Telem: ...
 
     def __pos__(self) -> Self: ...
     def __neg__(self) -> Self: ...
@@ -97,15 +97,41 @@ class flint_mpoly(flint_elem, Generic[Tctx, Telem, Telem_coerce]):
     def __rdivmod__(self, other: Telem | Telem_coerce | int) -> tuple[Self, Self]: ...
     def __pow__(self, other: Telem | Telem_coerce | int) -> Self: ...
     def __rpow__(self, other: Telem | Telem_coerce | int) -> Self: ...
-    def __iter__(self) -> Iterable[tuple[int, ...]]: ...
-    def __contains__(self, index: tuple[int, ...]) -> bool: ...
 
-    def unused_gens(self) -> tuple[str, ...]: ...
+    def iadd(self, other: Telem | Telem_coerce | int) -> None: ...
+    def isub(self, other: Telem | Telem_coerce | int) -> None: ...
+    def imul(self, other: Telem | Telem_coerce | int) -> None: ...
+
+    def gcd(self, other: Self) -> Self: ...
+    def term_content(self) -> Self: ...
+
+    def factor(self) -> tuple[Telem, Sequence[tuple[Self, int]]]: ...
+    def factor_squarefree(self) -> tuple[Telem, Sequence[tuple[Self, int]]]: ...
+
+    def sqrt(self) -> Self: ...
+
+    def resultant(self, other: Self, var: _str | int) -> Self: ...
+    def discriminant(self, var: _str | int) -> Self: ...
+
+    def deflation_index(self) -> tuple[list[int], list[int]]: ...
+    def deflation(self) -> tuple[Self, list[int]]: ...
+    def deflation_monom(self) -> tuple[Self, list[int], Self]: ...
+
+    def inflate(self, N: list[int]) -> Self: ...
+    def deflate(self, N: list[int]) -> Self: ...
+
+    def subs(self, mapping: dict[_str | int, Telem | Telem_coerce | int]) -> Self: ...
+    def compose(self, *args: Self, ctx: Tctx | None = None) -> Self: ...
+    def __call__(self, *args: Telem | Telem_coerce) -> Telem: ...
+
+    def derivative(self, var: _str | int) -> Self: ...
+
+    def unused_gens(self) -> tuple[_str, ...]: ...
 
     def project_to_context(
         self,
         other_ctx: Tctx,
-        mapping: dict[str | int, str | int] | None = None
+        mapping: dict[_str | int, _str | int] | None = None
     ) -> Self: ...
 
 
@@ -114,18 +140,18 @@ class flint_mpoly_context(flint_elem, Generic[Tmpoly, Telem, Telem_coerce]):
     def nvars(self) -> int: ...
     def ordering(self) -> Ordering: ...
 
-    def gen(self, i: int) -> Tmpoly: ...
-    def from_dict(self, d: Mapping[tuple[int, ...], Telem_coerce]) -> Tmpoly: ...
-    def constant(self, z: Telem_coerce) -> Tmpoly: ...
+    def gen(self, i: int, /) -> Tmpoly: ...
+    def from_dict(self, d: Mapping[tuple[int, ...], Telem_coerce], /) -> Tmpoly: ...
+    def constant(self, z: Telem_coerce, /) -> Tmpoly: ...
 
-    def name(self, i: int) -> str: ...
+    def name(self, i: int, /) -> str: ...
     def names(self) -> tuple[str]: ...
     def gens(self) -> tuple[Tmpoly, ...]: ...
-    def variable_to_index(self, var: str) -> int: ...
+    def variable_to_index(self, var: str, /) -> int: ...
     def term(self, coeff: Telem_coerce | None = None, exp_vec: Iterable[int] | None = None) -> Tmpoly: ...
-    def drop_gens(self, gens: Iterable[str | int]) -> Self: ...
-    def append_gens(self, gens: Iterable[str | int]) -> Self: ...
-    def infer_generator_mapping(self, ctx: flint_mpoly_context) -> dict[int, int]: ...
+    def drop_gens(self, gens: Iterable[str | int], /) -> Self: ...
+    def append_gens(self, gens: Iterable[str | int], /) -> Self: ...
+    def infer_generator_mapping(self, ctx: flint_mpoly_context, /) -> dict[int, int]: ...
 
     @classmethod
     def from_context(cls,
@@ -134,3 +160,14 @@ class flint_mpoly_context(flint_elem, Generic[Tmpoly, Telem, Telem_coerce]):
         ordering: Ordering | str = Ordering.lex,
     ) -> Sctx:
         ...
+
+
+class flint_mod_mpoly_context(flint_mpoly_context[Tmpoly, Telem, Telem_coerce]):
+    @abstractmethod
+    def modulus(self) -> int: ...
+
+
+class flint_series(flint_elem, Generic[Telem]):
+    """Base class for power series."""
+    def __iter__(self) -> Iterator[Telem]: ...
+    def coeffs(self) -> list[Telem]: ...
diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py
index bee2d5d0..98e14e20 100644
--- a/src/flint/test/test_all.py
+++ b/src/flint/test/test_all.py
@@ -1,19 +1,22 @@
-from typing import Any, Callable, TypeVar, Iterable, Protocol
+from __future__ import annotations
+from typing import Any, Callable, TypeVar, Iterable, Protocol, TYPE_CHECKING
+
 import math
 import operator
 import pickle
 import platform
 import random
 
+import flint
+import flint.flint_base.flint_base as flint_base
 from flint.utils.flint_exceptions import DomainError, IncompatibleContextError
 
-import flint
 
 PYPY = platform.python_implementation() == "PyPy"
 
 ctx = flint.ctx
 
-def raises(f, exception):
+def raises(f, exception) -> bool:
     try:
         f()
     except exception:
@@ -21,6 +24,16 @@ def raises(f, exception):
     return False
 
 
+if TYPE_CHECKING:
+    from typing import TypeIs
+
+
+Tscalar = TypeVar('Tscalar', bound=flint_base.flint_scalar)
+Tscalar_co = TypeVar('Tscalar_co', bound=flint_base.flint_scalar, covariant=True)
+Tmpoly = TypeVar('Tmpoly', bound=flint_base.flint_mpoly)
+Tmpolyctx_co = TypeVar('Tmpolyctx_co', bound=flint_base.flint_mpoly_context, covariant=True)
+
+
 _default_ctx_string = """\
 pretty = True      # pretty-print repr() output
 unicode = False    # use unicode characters in output
@@ -368,10 +381,10 @@ def test_fmpz_poly():
     assert Z([1]) == Z([1])
     assert Z([1]) == Z([flint.fmpz(1)])
     assert Z(Z([1,2])) == Z([1,2])
-    assert raises(lambda: Z([1,2,[]]), TypeError)
-    assert raises(lambda: Z({}), TypeError)
+    assert raises(lambda: Z([1,2,[]]), TypeError) # type: ignore
+    assert raises(lambda: Z({}), TypeError) # type: ignore
     # XXX: This should probably be made to work:
-    assert raises(lambda: Z((1,2,3)), TypeError)
+    assert raises(lambda: Z((1,2,3)), TypeError) # type: ignore
     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])
@@ -389,32 +402,32 @@ def test_fmpz_poly():
         assert Z([1,2,3]) ** ztype(1) == Z([1,2,3])
         assert Z([1,2,3]) ** ztype(2) == Z([1,4,10,12,9])
     assert divmod(Z([11,6,2]), Z([1,2])) == (Z([2,1]), Z([9,1]))
-    assert raises(lambda: pow(Z([1,2]), 2, 3), NotImplementedError)
+    assert raises(lambda: pow(Z([1,2]), 2, 3), NotImplementedError) # type: ignore
     assert +Z([1,2]) == Z([1,2])
     assert -Z([1,2]) == Z([-1,-2])
-    assert raises(lambda: Z([1,2]) + [], TypeError)
-    assert raises(lambda: Z([1,2]) - [], TypeError)
-    assert raises(lambda: Z([1,2]) * [], TypeError)
-    assert raises(lambda: Z([1,2]) / [], TypeError)
-    assert raises(lambda: Z([1,2]) // [], TypeError)
-    assert raises(lambda: Z([1,2]) % [], TypeError)
-    assert raises(lambda: divmod(Z([1,2]), []), TypeError)
-    assert raises(lambda: [] + Z([1,2]), TypeError)
-    assert raises(lambda: [] - Z([1,2]), TypeError)
-    assert raises(lambda: [] * Z([1,2]), TypeError)
-    assert raises(lambda: [] / Z([1,2]), TypeError)
-    assert raises(lambda: [] // Z([1,2]), TypeError)
-    assert raises(lambda: [] % Z([1,2]), TypeError)
-    assert raises(lambda: divmod([], Z([1,2])), TypeError)
-    assert raises(lambda: Z([1,2,3]) ** -1, DomainError)
-    assert raises(lambda: Z([1,2,3]) ** Z([1,2]), TypeError)
+    assert raises(lambda: Z([1,2]) + [], TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) - [], TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) * [], TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) / [], TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) // [], TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) % [], TypeError) # type: ignore
+    assert raises(lambda: divmod(Z([1,2]), []), TypeError) # type: ignore
+    assert raises(lambda: [] + Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: [] - Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: [] * Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: [] / Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: [] // Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: [] % Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: divmod([], Z([1,2])), TypeError) # type: ignore
+    assert raises(lambda: Z([1,2,3]) ** -1, DomainError) # type: ignore
+    assert raises(lambda: Z([1,2,3]) ** Z([1,2]), TypeError) # type: ignore
     assert raises(lambda: Z([1,2]) // Z([]), ZeroDivisionError)
     assert raises(lambda: Z([]) // Z([]), ZeroDivisionError)
     assert raises(lambda: Z([1,2]) % Z([]), ZeroDivisionError)
     assert raises(lambda: divmod(Z([1,2]), Z([])), ZeroDivisionError)
-    assert raises(lambda: Z([1,2]) < Z([1,2]), TypeError)
-    assert raises(lambda: Z([1,2]) < [], TypeError)
-    assert raises(lambda: [] < Z([1,2]), TypeError)
+    assert raises(lambda: Z([1,2]) < Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) < [], TypeError) # type: ignore
+    assert raises(lambda: [] < Z([1,2]), TypeError) # type: ignore
     assert Z([]).degree() == -1
     assert Z([]).length() == 0
     p = Z([1,2])
@@ -452,7 +465,7 @@ def test_fmpz_poly():
     assert p(flint.fmpq_poly([2,3],5)) == flint.fmpq_poly([27,24,9],5)
     assert p(flint.arb("1.1")).overlaps(flint.arb("13.45"))
     assert p(flint.acb("1.1", "1.2")).overlaps(flint.acb("6.25", "18.00"))
-    assert raises(lambda: p(None), TypeError)
+    assert raises(lambda: p(None), TypeError) # type: ignore
     assert Z([1,2,3]).derivative() == Z([2,6])
     assert Z([1,2,-4]).height_bits() == 3
     assert Z([1,2,-4]).height_bits(signed=True) == -3
@@ -477,7 +490,7 @@ def test_fmpz_poly_factor():
     Z = flint.fmpz_poly
     assert Z([1,2]).gcd(Z([3,4])) == 1
     assert Z([1,2,1]).gcd(Z([1,1])) == Z([1,1])
-    assert raises(lambda: Z([1,2,1]).gcd([]), TypeError)
+    assert raises(lambda: Z([1,2,1]).gcd([]), TypeError) # type: ignore
     assert Z([1,2,1]).factor() == (1, [(Z([1,1]), 2)])
 
 def test_fmpz_poly_functions():
@@ -686,7 +699,7 @@ def test_fmpz_series():
     assert Z(Z([1]))._equal_repr(Z([1]))
     assert Z(Zp([1]))._equal_repr(Z([1]))
     assert Z(1)._equal_repr(Z([1]))
-    assert raises(lambda: Z([1]) < Z([1]), TypeError)
+    assert raises(lambda: Z([1]) < Z([1]), TypeError) # type: ignore
     assert len(Z([1,2])) == 2
     assert Z([1,2]).length() == 2
     s5 = Z([1,2])
@@ -714,18 +727,18 @@ def set_bad():
     assert (Z([1,2]) + Zp([3,4,5]))._equal_repr(Z([4,6,5]))
     assert (1 + Z([1,2]))._equal_repr(Z([2,2]))
     assert (Zp([1,2]) + Z([3,4,5]))._equal_repr(Z([4,6,5]))
-    assert raises(lambda: Z([1,2]) + [], TypeError)
-    assert raises(lambda: [] + Z([1,2]), TypeError)
+    assert raises(lambda: Z([1,2]) + [], TypeError) # type: ignore
+    assert raises(lambda: [] + Z([1,2]), TypeError) # type: ignore
     assert (Z([1,2]) - Z([3,5]))._equal_repr(Z([-2,-3]))
     assert (Z([1,2]) - 1)._equal_repr(Z([0,2]))
     assert (1 - Z([1,2]))._equal_repr(Z([0,-2]))
-    assert raises(lambda: [] - Z([1,2]), TypeError)
-    assert raises(lambda: Z([1,2]) - [], TypeError)
+    assert raises(lambda: [] - Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) - [], TypeError) # type: ignore
     assert (Z([1,2]) * Z([1,2]))._equal_repr(Z([1,4,4]))
     assert (2 * Z([1,2]))._equal_repr(Z([2,4]))
     assert (Z([1,2]) * 2)._equal_repr(Z([2,4]))
-    assert raises(lambda: [] * Z([1,2]), TypeError)
-    assert raises(lambda: Z([1,2]) * [], TypeError)
+    assert raises(lambda: [] * Z([1,2]), TypeError) # type: ignore
+    assert raises(lambda: Z([1,2]) * [], TypeError) # type: ignore
     assert Z([1,2]).valuation() == 0
     assert Z([0,2]).valuation() == 1
     assert Z([0,0]).valuation() == -1
@@ -740,13 +753,13 @@ def set_bad():
     # XXX: This gives leading term not a unit:
     # assert Z([2,4]) / 2 == Z([1,2])
     assert (Z([1,1]) / -1)._equal_repr(Z([-1,-1]))
-    assert raises(lambda: Z([1,1]) / [], TypeError)
-    assert raises(lambda: [] / Z([1,1]), TypeError)
+    assert raises(lambda: Z([1,1]) / [], TypeError) # type: ignore
+    assert raises(lambda: [] / Z([1,1]), TypeError) # type: ignore
     assert (Z([1,2]) ** 2)._equal_repr(Z([1,4,4]))
-    assert raises(lambda: pow(Z([1,2]), 3, 5), NotImplementedError)
+    assert raises(lambda: pow(Z([1,2]), 3, 5), NotImplementedError) # type: ignore
     assert (Z([1,2])(Z([0,1,2])))._equal_repr(Z([1,2,4]))
     assert raises(lambda: Z([1,2])(Z([1,2])), ValueError)
-    assert raises(lambda: Z([1,2])([]), TypeError)
+    assert raises(lambda: Z([1,2])([]), TypeError) # type: ignore
     coeffs = [0, 1, -2, 8, -40, 224, -1344, 8448, -54912, 366080]
     assert Z([0,1,2]).reversion()._equal_repr(Z(coeffs))
     # power series reversion must have valuation 1
@@ -936,9 +949,9 @@ def test_fmpq_poly():
     assert Q() != Q([1])
     assert Q([1]) == Q([1])
     assert Q([1], 2) == Q([flint.fmpq(1,2)])
-    assert raises(lambda: Q([1,[]]), TypeError)
-    assert raises(lambda: Q({}), TypeError)
-    assert raises(lambda: Q([1], []), TypeError)
+    assert raises(lambda: Q([1,[]]), TypeError) # type: ignore
+    assert raises(lambda: Q({}), TypeError) # type: ignore
+    assert raises(lambda: Q([1], []), TypeError) # type: ignore
     assert raises(lambda: Q([1], 0), ZeroDivisionError)
     assert bool(Q()) is False
     assert bool(Q([1])) is True
@@ -948,21 +961,21 @@ def test_fmpq_poly():
     assert 3 + Q([1,2]) == Q([4,2])
     assert Q([1,2]) - 3 == Q([-2,2])
     assert 3 - Q([1,2]) == Q([2,-2])
-    assert raises(lambda: Q([1]) + [], TypeError)
-    assert raises(lambda: [] + Q([1]), TypeError)
-    assert raises(lambda: Q([1]) - [], TypeError)
-    assert raises(lambda: [] - Q([1]), TypeError)
-    assert raises(lambda: Q([1]) * [], TypeError)
-    assert raises(lambda: [] * Q([1]), TypeError)
+    assert raises(lambda: Q([1]) + [], TypeError) # type: ignore
+    assert raises(lambda: [] + Q([1]), TypeError) # type: ignore
+    assert raises(lambda: Q([1]) - [], TypeError) # type: ignore
+    assert raises(lambda: [] - Q([1]), TypeError) # type: ignore
+    assert raises(lambda: Q([1]) * [], TypeError) # type: ignore
+    assert raises(lambda: [] * Q([1]), TypeError) # type: ignore
     assert Q([1,2,1]) // Q([1,1]) == Q([1,1])
     assert Q([1,2,1]) % Q([1,1]) == 0
     assert divmod(Q([1,2,1]), Q([1,1])) == (Q([1,1]), 0)
-    assert raises(lambda: Q([1,2,1]) // [], TypeError)
-    assert raises(lambda: [] // Q([1,2,1]), TypeError)
-    assert raises(lambda: Q([1,2,1]) % [], TypeError)
-    assert raises(lambda: [] % Q([1,2,1]), TypeError)
-    assert raises(lambda: divmod(Q([1,2,1]), []), TypeError)
-    assert raises(lambda: divmod([], Q([1,2,1])), TypeError)
+    assert raises(lambda: Q([1,2,1]) // [], TypeError) # type: ignore
+    assert raises(lambda: [] // Q([1,2,1]), TypeError) # type: ignore
+    assert raises(lambda: Q([1,2,1]) % [], TypeError) # type: ignore
+    assert raises(lambda: [] % Q([1,2,1]), TypeError) # type: ignore
+    assert raises(lambda: divmod(Q([1,2,1]), []), TypeError) # type: ignore
+    assert raises(lambda: divmod([], Q([1,2,1])), TypeError) # type: ignore
     assert raises(lambda: Q([1,2,1]) / 0, ZeroDivisionError)
     assert raises(lambda: Q([1,2,1]) // 0, ZeroDivisionError)
     assert raises(lambda: Q([1,2,1]) % 0, ZeroDivisionError)
@@ -975,7 +988,7 @@ def test_fmpq_poly():
     assert Q([1,2]) != Z([3,2])
     assert Z([1,2]) != Q([3,2])
     assert Q([1,2]) != []
-    assert raises(lambda: Q([1,2]) < Q([1,2]), TypeError)
+    assert raises(lambda: Q([1,2]) < Q([1,2]), TypeError) # type: ignore
     assert Q([1,2,3])*Q([1,2]) == Q([1,4,7,6])
     assert Q([1,2,3])*Z([1,2]) == Q([1,4,7,6])
     assert Q([1,2,3]) * 3 == Q([3,6,9])
@@ -986,7 +999,7 @@ def test_fmpq_poly():
     assert raises(lambda: Q([1,2]) / Q([2,2]), DomainError)
     assert Q([1,2,3]) / flint.fmpq(2,3) == Q([1,2,3]) * flint.fmpq(3,2)
     assert Q([1,2,3]) ** 2 == Q([1,2,3]) * Q([1,2,3])
-    assert raises(lambda: pow(Q([1,2]), 3, 5), NotImplementedError)
+    assert raises(lambda: pow(Q([1,2]), 3, 5), NotImplementedError) # type: ignore
     assert Q([1,2,flint.fmpq(1,2)]).coeffs() == [1,2,flint.fmpq(1,2)]
     assert Q().coeffs() == []
     assert Q().degree() == -1
@@ -1017,13 +1030,13 @@ def set_bad():
     assert p(flint.fmpq(2,3)) == flint.fmpq(71,63)
     assert p(Z([1,-1])) == Q([12,-14,5],7)
     assert p(flint.fmpq_poly([2,3],5)) == flint.fmpq_poly([27,24,9],7*5)
-    assert raises(lambda: p([]), TypeError)
+    assert raises(lambda: p([]), TypeError) # type: ignore
     assert Q([1,2,3]).derivative() == Q([2,6])
     assert Q([1,2,3]).integral() == Q([0,1,1,1])
     assert Q([1,2,1]).gcd(Q([1,1])) == Q([1,1])
     assert Q([1,2,1]).xgcd(Q([1,1])) == (Q([1,1]), 0, 1)
-    assert raises(lambda: Q([1,2,1]).gcd([]), TypeError)
-    assert raises(lambda: Q([1,2,1]).xgcd([]), TypeError)
+    assert raises(lambda: Q([1,2,1]).gcd([]), TypeError) # type: ignore
+    assert raises(lambda: Q([1,2,1]).xgcd([]), TypeError) # type: ignore
     assert Q([1,2,1]).factor() == (1, [(Q([1,1]), 2)])
     assert Q.bernoulli_poly(3) == Q([0,1,-3,2],2)
     assert Q.euler_poly(3) == Q([1,0,-6,4],4)
@@ -1193,9 +1206,9 @@ def test_fmpq_series():
     assert Q([1],2)._equal_repr(Q([flint.fmpq(1,2)]))
     assert not Q([1],1,10)._equal_repr(Q([1],1,11))
     assert Q([1,2],3)._equal_repr(Q([flint.fmpq(1,3), flint.fmpq(2,3)]))
-    assert raises(lambda: Q([1],[]), TypeError)
+    assert raises(lambda: Q([1],[]), TypeError) # type: ignore
     assert raises(lambda: Q([1],0), ZeroDivisionError)
-    assert raises(lambda: Q([1]) < Q([1]), TypeError)
+    assert raises(lambda: Q([1]) < Q([1]), TypeError) # type: ignore
     assert len(Q([1,2])) == 2
     assert Q([1,2]).length() == 2
     s5 = Q([1,2])
@@ -1225,18 +1238,18 @@ def set_bad():
     assert (Q([1,2]) + Qp([3,4,5]))._equal_repr(Q([4,6,5]))
     assert (1 + Q([1,2]))._equal_repr(Q([2,2]))
     assert (Qp([1,2]) + Q([3,4,5]))._equal_repr(Q([4,6,5]))
-    assert raises(lambda: Q([1,2]) + [], TypeError)
-    assert raises(lambda: [] + Q([1,2]), TypeError)
+    assert raises(lambda: Q([1,2]) + [], TypeError) # type: ignore
+    assert raises(lambda: [] + Q([1,2]), TypeError) # type: ignore
     assert (Q([1,2]) - Q([3,5]))._equal_repr(Q([-2,-3]))
     assert (Q([1,2]) - 1)._equal_repr(Q([0,2]))
     assert (1 - Q([1,2]))._equal_repr(Q([0,-2]))
-    assert raises(lambda: [] - Q([1,2]), TypeError)
-    assert raises(lambda: Q([1,2]) - [], TypeError)
+    assert raises(lambda: [] - Q([1,2]), TypeError) # type: ignore
+    assert raises(lambda: Q([1,2]) - [], TypeError) # type: ignore
     assert (Q([1,2]) * Q([1,2]))._equal_repr(Q([1,4,4]))
     assert (2 * Q([1,2]))._equal_repr(Q([2,4]))
     assert (Q([1,2]) * 2)._equal_repr(Q([2,4]))
-    assert raises(lambda: [] * Q([1,2]), TypeError)
-    assert raises(lambda: Q([1,2]) * [], TypeError)
+    assert raises(lambda: [] * Q([1,2]), TypeError) # type: ignore
+    assert raises(lambda: Q([1,2]) * [], TypeError) # type: ignore
     assert Q([1,2]).valuation() == 0
     assert Q([0,2]).valuation() == 1
     assert Q([0,0]).valuation() == -1
@@ -1250,8 +1263,8 @@ def set_bad():
     assert (Q([2,4]) / 2)._equal_repr(Q([1,2]))
     assert (Q([1,4]) / 2)._equal_repr(Q([1,4],2))
     assert (Q([1,1]) / -1)._equal_repr(Q([-1,-1]))
-    assert raises(lambda: Q([1,1]) / [], TypeError)
-    assert raises(lambda: [] / Q([1,1]), TypeError)
+    assert raises(lambda: Q([1,1]) / [], TypeError) # type: ignore
+    assert raises(lambda: [] / Q([1,1]), TypeError) # type: ignore
     q = Q([1,2],3)
     assert (q ** 0)._equal_repr(Q([1]))
     assert (q ** 1)._equal_repr(q)
@@ -1260,10 +1273,10 @@ def set_bad():
     assert (q ** 4)._equal_repr(Q([1, 8, 24, 32, 16], 81))
     assert (Q([1,2]) ** 2)._equal_repr(Q([1,4,4]))
     assert (Q([1,2]) ** 2)._equal_repr(Q([1,4,4]))
-    assert raises(lambda: pow(Q([1,2]), 3, 5), NotImplementedError)
+    assert raises(lambda: pow(Q([1,2]), 3, 5), NotImplementedError) # type: ignore
     assert (Q([1,2])(Q([0,1,2])))._equal_repr(Q([1,2,4]))
     assert raises(lambda: Q([1,2])(Q([1,2])), ValueError)
-    assert raises(lambda: Q([1,2])([]), TypeError)
+    assert raises(lambda: Q([1,2])([]), TypeError) # type: ignore
     coeffs = [0, 1, -2, 8, -40, 224, -1344, 8448, -54912, 366080]
     assert Q([0,1,2]).reversion()._equal_repr(Q(coeffs))
     # power series reversion must have valuation 1
@@ -1669,7 +1682,7 @@ def test_fmpz_mod():
     assert F_cmp.is_prime() is False
 
     # Context tests
-    assert raises(lambda: fmpz_mod_ctx("AAA"), TypeError)
+    assert raises(lambda: fmpz_mod_ctx("AAA"), TypeError) # type: ignore
     assert raises(lambda: fmpz_mod_ctx(-1), ValueError)
     assert F_sml.modulus() == p_sml
     assert F_med.modulus() == p_med
@@ -1687,7 +1700,7 @@ def test_fmpz_mod():
     assert str(F_sml) == "Context for fmpz_mod with modulus: 163"
 
     # Type tests
-    assert raises(lambda: fmpz_mod(1, "AAA"), TypeError)
+    assert raises(lambda: fmpz_mod(1, "AAA"), TypeError) # type: ignore
 
     # Test for small, medium and large char.
     for F_test in [F_sml, F_med, F_big]:
@@ -1700,15 +1713,15 @@ def test_fmpz_mod():
         assert int(test_z) >= 0
 
         assert raises(lambda: F_test(F_cmp(1)), ValueError)
-        assert raises(lambda: F_test("abc"), NotImplementedError)
+        assert raises(lambda: F_test("abc"), NotImplementedError) # type: ignore
 
         F_test_copy = fmpz_mod_ctx(test_mod)
         F_other = fmpz_mod_ctx(11)
 
-        assert raises(lambda: F_test(test_x) > 0, TypeError)
-        assert raises(lambda: F_test(test_x) >= 0, TypeError)
-        assert raises(lambda: F_test(test_x) < 0, TypeError)
-        assert raises(lambda: F_test(test_x) <= 0, TypeError)
+        assert raises(lambda: F_test(test_x) > 0, TypeError) # type: ignore
+        assert raises(lambda: F_test(test_x) >= 0, TypeError) # type: ignore
+        assert raises(lambda: F_test(test_x) < 0, TypeError) # type: ignore
+        assert raises(lambda: F_test(test_x) <= 0, TypeError) # type: ignore
 
         assert (test_x == F_test(test_x)) is True, f"{test_x}, {F_test(test_x)}"
         assert (124 != F_test(test_x)) is True
@@ -1759,7 +1772,7 @@ def test_fmpz_mod():
         assert F_test(test_x) + F_test(test_y) == F_test(test_x + test_y)
         assert F_test(test_x) + F_test_copy(test_y) == F_test(test_x + test_y)
         assert F_test(test_x) + F_test(test_y) == F_test_copy(test_x + test_y)
-        assert raises(lambda: F_test(test_x) + "AAA", TypeError)
+        assert raises(lambda: F_test(test_x) + "AAA", TypeError) # type: ignore
 
         assert F_test(test_x) + F_test(test_y) == F_test(test_y) + F_test(test_x)
         assert F_test(test_x) + test_y == F_test(test_x + test_y)
@@ -1777,8 +1790,8 @@ def test_fmpz_mod():
         assert test_y - F_test(test_x) == F_test(test_y) - F_test(test_x)
         assert F_test(test_x) - fmpz(test_y) == F_test(test_x) - F_test(test_y)
         assert raises(lambda: F_test(test_x) - F_other(test_y), ValueError)
-        assert raises(lambda: F_test(test_x) - "AAA", TypeError)
-        assert raises(lambda: "AAA" - F_test(test_x), TypeError)
+        assert raises(lambda: F_test(test_x) - "AAA", TypeError) # type: ignore
+        assert raises(lambda: "AAA" - F_test(test_x), TypeError) # type: ignore
 
         # Multiplication
 
@@ -1798,7 +1811,7 @@ def test_fmpz_mod():
         assert fmpz(1) * F_test(test_x) == F_test(1 * test_x)
         assert fmpz(2) * F_test(test_x) == F_test(2 * test_x)
         assert fmpz(3) * F_test(test_x) == F_test(3 * test_x)
-        assert raises(lambda: F_test(test_x) * "AAA", TypeError)
+        assert raises(lambda: F_test(test_x) * "AAA", TypeError) # type: ignore
         assert raises(lambda: F_test(test_x) * F_other(test_x), ValueError)
 
         # Exponentiation
@@ -1807,7 +1820,7 @@ def test_fmpz_mod():
         assert F_test(0)**1 == pow(0, 1, test_mod)
         assert F_test(0)**2 == pow(0, 2, test_mod)
         assert raises(lambda: F_test(0)**(-1), ZeroDivisionError)
-        assert raises(lambda: F_test(0)**("AA"), NotImplementedError)
+        assert raises(lambda: F_test(0)**("AA"), TypeError) # type: ignore
 
         assert F_test(test_x)**fmpz(0) == pow(test_x, 0, test_mod)
         assert F_test(test_x)**fmpz(1) == pow(test_x, 1, test_mod)
@@ -1841,10 +1854,10 @@ def test_fmpz_mod():
         assert F_test(test_x) / F_test(test_y) == (test_x * pow(test_y, -1, test_mod)) % test_mod
         assert F_test(test_x) / fmpz(test_y) == (test_x * pow(test_y, -1, test_mod)) % test_mod
         assert F_test(test_x) / test_y == (test_x * pow(test_y, -1, test_mod)) % test_mod
-        assert raises(lambda: F_test(test_x) / "AAA", TypeError)
-        assert raises(lambda: "AAA" / F_test(test_x), TypeError)
+        assert raises(lambda: F_test(test_x) / "AAA", TypeError) # type: ignore
+        assert raises(lambda: "AAA" / F_test(test_x), TypeError) # type: ignore
         assert raises(lambda: F_other(test_x) / F_test(test_x), ValueError)
-        assert raises(lambda: F_test(test_x) // F_test(test_x), TypeError)
+        assert raises(lambda: F_test(test_x) // F_test(test_x), TypeError) # type: ignore
         assert 1 / F_test(2) ==  pow(2, -1, test_mod)
         assert 1 / F_test(test_x) ==  pow(test_x, -1, test_mod)
         assert 1 / F_test(test_y) ==  pow(test_y, -1, test_mod)
@@ -1869,7 +1882,7 @@ def test_fmpz_mod_dlog():
 
     # Need to use either fmpz_mod or something which can be case to
     # fmpz
-    assert raises(lambda: g.discrete_log("A"), TypeError)
+    assert raises(lambda: g.discrete_log("A"), TypeError) # type: ignore
 
     F = fmpz_mod_ctx(163)
     g = F(2)
@@ -1904,7 +1917,7 @@ def test_fmpz_mod_poly():
     R2 = fmpz_mod_poly_ctx(11)
     R3 = fmpz_mod_poly_ctx(13)
 
-    assert raises(lambda: fmpz_mod_ctx("AAA"), TypeError)
+    assert raises(lambda: fmpz_mod_ctx("AAA"), TypeError) # type: ignore
     assert raises(lambda: fmpz_mod_ctx(-1), ValueError)
     assert (R1 == R1) is True
     assert (R1 != R1) is False
@@ -2570,6 +2583,12 @@ def _all_polys() -> list[tuple[Any, Any, bool, flint.fmpz]]:
 
 
 def test_polys():
+    # To test type annotations, uncomment:
+    # P: type[flint.fmpq_poly]
+    # S: type[flint.fmpq]
+    # is_field: bool
+    # characteristic: flint.fmpz
+
     for P, S, is_field, characteristic in _all_polys():
 
         composite_characteristic = characteristic != 0 and not characteristic.is_prime()
@@ -2580,12 +2599,12 @@ def test_polys():
 
         assert P([S(1)]) == P([1]) == P(P([1])) == P(1)
 
-        assert raises(lambda: P([None]), TypeError)
-        assert raises(lambda: P(object()), TypeError)
-        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), TypeError)
+        assert raises(lambda: P([None]), TypeError) # type: ignore
+        assert raises(lambda: P(object()), TypeError) # type: ignore
+        assert raises(lambda: P(None), TypeError) # type: ignore
+        assert raises(lambda: P(None, None), TypeError) # type: ignore
+        assert raises(lambda: P([1,2], None), TypeError) # type: ignore
+        assert raises(lambda: P(1, None), TypeError) # type: ignore
 
         assert len(P([])) == P([]).length() == 0
         assert len(P([1])) == P([1]).length() == 1
@@ -2629,18 +2648,18 @@ def test_polys():
         assert (None == P([1])) is False
         assert (None != P([1])) is True
 
-        assert raises(lambda: P([1]) < P([1]), TypeError)
-        assert raises(lambda: P([1]) <= P([1]), TypeError)
-        assert raises(lambda: P([1]) > P([1]), TypeError)
-        assert raises(lambda: P([1]) >= P([1]), TypeError)
-        assert raises(lambda: P([1]) < None, TypeError)
-        assert raises(lambda: P([1]) <= None, TypeError)
-        assert raises(lambda: P([1]) > None, TypeError)
-        assert raises(lambda: P([1]) >= None, TypeError)
-        assert raises(lambda: None < P([1]), TypeError)
-        assert raises(lambda: None <= P([1]), TypeError)
-        assert raises(lambda: None > P([1]), TypeError)
-        assert raises(lambda: None >= P([1]), TypeError)
+        assert raises(lambda: P([1]) < P([1]), TypeError) # type: ignore
+        assert raises(lambda: P([1]) <= P([1]), TypeError) # type: ignore
+        assert raises(lambda: P([1]) > P([1]), TypeError) # type: ignore
+        assert raises(lambda: P([1]) >= P([1]), TypeError) # type: ignore
+        assert raises(lambda: P([1]) < None, TypeError) # type: ignore
+        assert raises(lambda: P([1]) <= None, TypeError) # type: ignore
+        assert raises(lambda: P([1]) > None, TypeError) # type: ignore
+        assert raises(lambda: P([1]) >= None, TypeError) # type: ignore
+        assert raises(lambda: None < P([1]), TypeError) # type: ignore
+        assert raises(lambda: None <= P([1]), TypeError) # type: ignore
+        assert raises(lambda: None > P([1]), TypeError) # type: ignore
+        assert raises(lambda: None >= P([1]), TypeError) # type: ignore
 
         assert P([1, 2, 3])[1] == S(2)
         assert P([1, 2, 3])[-1] == S(0)
@@ -2680,7 +2699,7 @@ def setbad(obj, i, val):
         assert p(0) == p(S(0)) == S(1) == 1
         assert p(1) == p(S(1)) == S(6) == 6
         assert p(p) == P([6, 16, 36, 36, 27])
-        assert raises(lambda: p(None), TypeError)
+        assert raises(lambda: p(None), TypeError) # type: ignore
 
         assert bool(P([])) is False
         assert bool(P([1])) is True
@@ -2700,8 +2719,8 @@ def setbad(obj, i, val):
             assert P([1, 2, 3]) + T(1) == P([2, 2, 3])
             assert T(1) + P([1, 2, 3]) == P([2, 2, 3])
 
-        assert raises(lambda: P([1, 2, 3]) + None, TypeError)
-        assert raises(lambda: None + P([1, 2, 3]), TypeError)
+        assert raises(lambda: P([1, 2, 3]) + None, TypeError) # type: ignore
+        assert raises(lambda: None + P([1, 2, 3]), TypeError) # type: ignore
 
         assert P([1, 2, 3]) - P([4, 5, 6]) == P([-3, -3, -3])
 
@@ -2709,8 +2728,8 @@ def setbad(obj, i, val):
             assert P([1, 2, 3]) - T(1) == P([0, 2, 3])
             assert T(1) - P([1, 2, 3]) == P([0, -2, -3])
 
-        assert raises(lambda: P([1, 2, 3]) - None, TypeError)
-        assert raises(lambda: None - P([1, 2, 3]), TypeError)
+        assert raises(lambda: P([1, 2, 3]) - None, TypeError) # type: ignore
+        assert raises(lambda: None - P([1, 2, 3]), TypeError) # type: ignore
 
         assert P([1, 2, 3]) * P([4, 5, 6]) == P([4, 13, 28, 27, 18])
 
@@ -2718,8 +2737,8 @@ def setbad(obj, i, val):
             assert P([1, 2, 3]) * T(2) == P([2, 4, 6])
             assert T(2) * P([1, 2, 3]) == P([2, 4, 6])
 
-        assert raises(lambda: P([1, 2, 3]) * None, TypeError)
-        assert raises(lambda: None * P([1, 2, 3]), TypeError)
+        assert raises(lambda: P([1, 2, 3]) * None, TypeError) # type: ignore
+        assert raises(lambda: None * P([1, 2, 3]), TypeError) # type: ignore
 
         assert P([1, 2, 1]) // P([1, 1]) == P([1, 1])
         assert P([1, 2, 1]) % P([1, 1]) == P([0])
@@ -2745,13 +2764,13 @@ def setbad(obj, i, val):
         assert 1 % P([1, 1]) == P([1])
         assert divmod(1, P([1, 1])) == (P([0]), P([1]))
 
-        assert raises(lambda: P([1, 2, 1]) // None, TypeError)
-        assert raises(lambda: P([1, 2, 1]) % None, TypeError)
-        assert raises(lambda: divmod(P([1, 2, 1]), None), TypeError)
+        assert raises(lambda: P([1, 2, 1]) // None, TypeError) # type: ignore
+        assert raises(lambda: P([1, 2, 1]) % None, TypeError) # type: ignore
+        assert raises(lambda: divmod(P([1, 2, 1]), None), TypeError) # type: ignore
 
-        assert raises(lambda: None // P([1, 1]), TypeError)
-        assert raises(lambda: None % P([1, 1]), TypeError)
-        assert raises(lambda: divmod(None, P([1, 1])), TypeError)
+        assert raises(lambda: None // P([1, 1]), TypeError) # type: ignore
+        assert raises(lambda: None % P([1, 1]), TypeError) # type: ignore
+        assert raises(lambda: divmod(None, P([1, 1])), TypeError) # type: ignore
 
         assert raises(lambda: P([1, 2, 1]) // 0, ZeroDivisionError)
         assert raises(lambda: P([1, 2, 1]) % 0, ZeroDivisionError)
@@ -2786,32 +2805,32 @@ def setbad(obj, i, val):
         assert P([1, 1]) ** 1 == P([1, 1])
         assert P([1, 1]) ** 2 == P([1, 2, 1])
         assert raises(lambda: P([1, 1]) ** -1, DomainError)
-        assert raises(lambda: P([1, 1]) ** None, TypeError)
+        assert raises(lambda: P([1, 1]) ** None, TypeError) # type: ignore
 
         # 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, flint.fq_default_poly]:
-            assert raises(lambda: pow(p, 2, mod), NotImplementedError)
+            assert raises(lambda: pow(p, 2, mod), NotImplementedError) # type: ignore
         else:
             assert p * p % mod == pow(p, 2, mod)
 
         if not composite_characteristic:
             assert P([1, 2, 1]).gcd(P([1, 1])) == P([1, 1])
-            assert raises(lambda: P([1, 2, 1]).gcd(None), TypeError)
+            assert raises(lambda: P([1, 2, 1]).gcd(None), TypeError) # type: ignore
         elif nmod_poly_will_crash:
             pass
         else:
             # Z/nZ for n not prime
             assert raises(lambda: P([1, 2, 1]).gcd(P([1, 1])), DomainError)
-            assert raises(lambda: P([1, 2, 1]).gcd(None), TypeError)
+            assert raises(lambda: P([1, 2, 1]).gcd(None), TypeError) # type: ignore
 
         if is_field:
             p1 = P([1, 0, 1])
             p2 = P([2, 1])
             g, s, t = P([1]), P([1])/5, P([2, -1])/5
             assert p1.xgcd(p2) == (g, s, t)
-            assert raises(lambda: p1.xgcd(None), TypeError)
+            assert raises(lambda: p1.xgcd(None), TypeError) # type: ignore
 
         if not composite_characteristic:
             assert P([1, 2, 1]).factor() == (S(1), [(P([1, 1]), 2)])
@@ -2882,29 +2901,6 @@ def test_poly_resultants():
                 assert a.resultant(b) == prime**tot
 
 
-#
-# It is not really possible to explain what _all_mpolys returns with type
-# annotations...
-#
-Tctx = TypeVar('Tctx', covariant=True)
-Tval = TypeVar('Tval')
-
-
-class _get_ctx(Protocol[Tctx]):
-    def __call__(self,
-        args: Iterable[str | tuple[str, int]] | tuple[str, int],
-        ordering: str | flint.Ordering = "lex"
-    ) -> Tctx:
-        ...
-# _get_ctx = Callable[[Iterable[str | tuple[str, int]] | tuple[str, int]], Tctx]
-_get_val = Callable[[int], Tval]
-_all_mpolys_type = tuple[
-    tuple[type[flint.fmpz_mpoly], _get_ctx[flint.fmpz_mpoly_ctx], _get_val[flint.fmpz], bool, flint.fmpz],
-    tuple[type[flint.fmpq_mpoly], _get_ctx[flint.fmpq_mpoly_ctx], _get_val[flint.fmpq], bool, flint.fmpz],
-    tuple[type[flint.nmod_mpoly], _get_ctx[flint.nmod_mpoly_ctx], _get_val[flint.nmod], bool, flint.fmpz],
-]
-
-
 def _all_mpolys(): # -> _all_mpolys_type:
     return (
         (flint.fmpz_mpoly, flint.fmpz_mpoly_ctx.get, flint.fmpz, False, flint.fmpz(0)),
@@ -2940,494 +2936,595 @@ def _all_mpolys(): # -> _all_mpolys_type:
     )
 
 
-def test_mpolys():
-    for P, get_context, S, is_field, characteristic in _all_mpolys():
-        # Division under modulo will raise a flint exception if something is
-        # not invertible, crashing the program. We can't tell before what is
-        # invertible and what is not before hand so we always raise an
-        # exception, except for fmpz_mpoly, that returns an bool noting if the
-        # division is exact or not.
-        composite_characteristic = characteristic != 0 and not characteristic.is_prime()
+class _GetMPolyCtx(Protocol[Tmpolyctx_co]):
+    def __call__(self,
+        names: Iterable[str | tuple[str, int]] | tuple[str, int],
+        ordering: str | flint.Ordering = "lex"
+    ) -> Tmpolyctx_co:
+        ...
 
-        ctx = get_context(("x", 2))
-
-        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 raises(lambda : ctx.__class__("x", flint.Ordering.lex), RuntimeError) # type: ignore
-        assert raises(lambda: get_context(("x", 2), ordering="bad"), ValueError) # type: ignore
-        assert raises(lambda: get_context(("x", -1)), ValueError)
-        assert raises(lambda: ctx.constant("bad"), TypeError) # type: ignore
-        assert raises(lambda: ctx.from_dict("bad"), ValueError) # type: ignore
-        assert raises(lambda: ctx.from_dict({(0, 0): "bad"}), TypeError) # type: ignore
-        assert raises(lambda: ctx.from_dict({(0, "bad"): 1}), TypeError) # type: ignore
-        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=get_context(("x",)).constant(0), ctx=ctx), IncompatibleContextError)
-        assert raises(lambda: P(val={}, ctx=None), ValueError)
-        assert raises(lambda: P(val={"bad": 1}, ctx=None), ValueError) # type: ignore
-        assert raises(lambda: P(val="1", ctx=None), ValueError)
-
-        ctx1 = get_context(("x", 4))
-        ctx2 = get_context(("x", 4), ordering="deglex")
-        assert ctx1.drop_gens(ctx1.names()).names() == tuple()
-        assert ctx1.drop_gens((ctx1.name(1), ctx1.name(2))).names() == (ctx1.name(0), ctx1.name(3))
-        assert ctx1.drop_gens(tuple()).names() == ctx1.names()
-        assert ctx1.drop_gens((-1,)).names() == ctx1.names()[:-1]
-
-        assert ctx.infer_generator_mapping(ctx) == {i: i for i in range(ctx.nvars())}
-        assert ctx1.infer_generator_mapping(ctx) == {0: 0, 1: 1}
-        assert ctx1.drop_gens(ctx.names()).infer_generator_mapping(ctx) == {}
-
-        assert quick_poly().project_to_context(ctx1) == \
-            ctx1.from_dict(
-                {(0, 0, 0, 0): 1, (0, 1, 0, 0): 2, (1, 0, 0, 0): 3, (2, 2, 0, 0): 4}
-            )
-        new_poly = quick_poly().project_to_context(ctx1)
-        assert ctx1.drop_gens(new_poly.unused_gens()) == ctx
-        assert new_poly.project_to_context(ctx) == quick_poly()
-
-        new_poly = quick_poly().project_to_context(ctx2)
-        new_ctx = ctx2.drop_gens(new_poly.unused_gens())
-        assert new_ctx != ctx
-        assert new_poly != quick_poly()
-
-        new_ctx = new_ctx.from_context(new_ctx, ordering=ctx.ordering())
-        assert new_ctx == ctx
-        assert new_poly.project_to_context(new_ctx) == quick_poly()
-
-        assert ctx.append_gens(*ctx1.names()[-2:]) == ctx1
-
-        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})
-
-        assert ctx.nvars() == 2
-        assert ctx.ordering() == flint.Ordering.lex
-
-        ctx1 = get_context(("x", 4))
-        assert [ctx1.name(i) for i in range(4)] == ['x0', 'x1', 'x2', 'x3']
-        for order in list(flint.Ordering):
-            ctx1 = get_context(("x", 4), ordering=order)
-            assert ctx1.ordering() == order
-
-        assert ctx.constant(1) == mpoly({(0, 0): 1}) == P(1, ctx=ctx)
 
-        assert raises(lambda: P([None]), TypeError) # type: ignore
-        assert raises(lambda: P(object()), TypeError) # type: ignore
-        assert raises(lambda: P(None), TypeError) # type: ignore
-        assert raises(lambda: P(None, None), TypeError) # type: ignore
-        assert raises(lambda: P([1,2], None), TypeError) # type: ignore
-        assert raises(lambda: P(1, None), ValueError) # type: ignore
-
-        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 True
-        assert (P(1, ctx=ctx) != 1) is False
-        assert (1 == P(1, ctx=ctx)) is True
-        assert (1 != P(1, ctx=ctx)) is False
-
-        assert (P(1, ctx=ctx) == S(1)) is True
-        assert (P(1, ctx=ctx) != S(1)) is False
-        assert (S(1) == P(1, ctx=ctx)) is True
-        assert (S(1) != P(1, ctx=ctx)) is False
-
-        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})
-
-        # Prevent pyright from getting confused in type inference:
-        P2 = P
-        if P2 is flint.fmpq_mpoly:
-            ctx_z = flint.fmpz_mpoly_ctx.get((("x", 2),))
-            assert quick_poly() == P2(ctx_z.from_dict({(0, 0): 1, (0, 1): 2, (1, 0): 3, (2, 2): 4}))
-            assert isinstance(ctx, flint.fmpq_mpoly_ctx)
-            assert P2(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) # type: ignore
-        assert raises(lambda: P(ctx=ctx) <= P(ctx=ctx), TypeError) # type: ignore
-        assert raises(lambda: P(ctx=ctx) > P(ctx=ctx), TypeError) # type: ignore
-        assert raises(lambda: P(ctx=ctx) >= P(ctx=ctx), TypeError) # type: ignore
-        assert raises(lambda: P(ctx=ctx) < None, TypeError) # type: ignore
-        assert raises(lambda: P(ctx=ctx) <= None, TypeError) # type: ignore
-        assert raises(lambda: P(ctx=ctx) > None, TypeError) # type: ignore
-        assert raises(lambda: P(ctx=ctx) >= None, TypeError) # type: ignore
-        assert raises(lambda: None < P(ctx=ctx), TypeError) # type: ignore
-        assert raises(lambda: None <= P(ctx=ctx), TypeError) # type: ignore
-        assert raises(lambda: None > P(ctx=ctx), TypeError) # type: ignore
-        assert raises(lambda: None >= P(ctx=ctx), TypeError) # type: ignore
+_MPolyTestCase = tuple[
+    type[Tmpoly],
+    _GetMPolyCtx['flint_base.flint_mpoly_context[Tmpoly, Tscalar, int]'],
+    Callable[[int], Tscalar],
+    bool,
+    flint.fmpz
+]
 
-        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) # type: ignore
-        assert raises(lambda: p[4], TypeError) # type: ignore
+class _Q(Protocol[Tscalar_co]):
+    def __call__(self, a: int, b: int | None = None, /) -> Tscalar_co:
+        ...
 
-        assert p[(2, 2)] == 4
-        assert p[(0, 0)] == 1
-        assert raises(lambda: p[(1,)], ValueError)
-        assert raises(lambda: p[(1, "bad")], TypeError) # type: ignore
-        assert raises(lambda: p["bad"], TypeError) # type: ignore
 
-        p = quick_poly()
-        p[(1, 0)] = S(10)
-        assert p == mpoly({(0, 0): 1, (0, 1): 2, (1, 0): 10, (2, 2): 4})
+def _is_Q(typ: object) -> TypeIs[_Q]:
+    return typ is flint.fmpq
+
+
+def _for_all_mpolys(test: Callable[[_MPolyTestCase], None]) -> None:
+    """Test all mpoly types with the given test function."""
+    # Spell it out like this so that a type checker can understand the types
+    # in the generics for each call of test().
+    test((
+        flint.fmpz_mpoly,
+        flint.fmpz_mpoly_ctx.get,
+        flint.fmpz,
+        False,
+        flint.fmpz(0)
+    ))
+    test((
+        flint.fmpq_mpoly,
+        flint.fmpq_mpoly_ctx.get,
+        flint.fmpq,
+        True,
+        flint.fmpz(0)
+    ))
+    test((
+        flint.fmpz_mod_mpoly,
+        lambda *args, **kwargs: flint.fmpz_mod_mpoly_ctx.get(*args, **kwargs, modulus=101),
+        lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(101)),
+        True,
+        flint.fmpz(101),
+    ))
+    test((
+        flint.fmpz_mod_mpoly,
+        lambda *args, **kwargs: flint.fmpz_mod_mpoly_ctx.get(*args, **kwargs, modulus=100),
+        lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(100)),
+        False,
+        flint.fmpz(100),
+    ))
+    test((
+        flint.nmod_mpoly,
+        lambda *args, **kwargs: flint.nmod_mpoly_ctx.get(*args, **kwargs, modulus=101),
+        lambda x: flint.nmod(x, 101),
+        True,
+        flint.fmpz(101),
+    ))
+    test((
+        flint.nmod_mpoly,
+        lambda *args, **kwargs: flint.nmod_mpoly_ctx.get(*args, **kwargs, modulus=100),
+        lambda x: flint.nmod(x, 100),
+        False,
+        flint.fmpz(100),
+    ))
+
+
+def all_mpolys(f: Callable[[_MPolyTestCase], None]) -> Callable[[], None]:
+    """Decorator to run a test function for all mpoly types."""
+    def wrapper():
+        _for_all_mpolys(f)
+    wrapper.__name__ = f.__name__
+    return wrapper
+
+
+@all_mpolys
+def test_mpolys_constructor(args: _MPolyTestCase[Tmpoly, Tscalar]) -> None:
+    P, get_context, S, _, _ = args
+
+    ctx = get_context(("x", 2))
+
+    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 raises(lambda : ctx.__class__("x", flint.Ordering.lex), RuntimeError) # type: ignore
+    assert raises(lambda: get_context(("x", 2), ordering="bad"), ValueError) # type: ignore
+    assert raises(lambda: get_context(("x", -1)), ValueError)
+    assert raises(lambda: ctx.constant("bad"), TypeError) # type: ignore
+    assert raises(lambda: ctx.from_dict("bad"), ValueError) # type: ignore
+    assert raises(lambda: ctx.from_dict({(0, 0): "bad"}), TypeError) # type: ignore
+    assert raises(lambda: ctx.from_dict({(0, "bad"): 1}), TypeError) # type: ignore
+    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=get_context(("x",)).constant(0), ctx=ctx), IncompatibleContextError) # type: ignore
+    assert raises(lambda: P(val={}, ctx=None), ValueError)
+    assert raises(lambda: P(val={"bad": 1}, ctx=None), ValueError) # type: ignore
+    assert raises(lambda: P(val="1", ctx=None), ValueError)
+
+    ctx1 = get_context(("x", 4))
+    ctx2 = get_context(("x", 4), ordering="deglex")
+    assert ctx1.drop_gens(ctx1.names()).names() == tuple()
+    assert ctx1.drop_gens((ctx1.name(1), ctx1.name(2))).names() == (ctx1.name(0), ctx1.name(3))
+    assert ctx1.drop_gens(tuple()).names() == ctx1.names()
+    assert ctx1.drop_gens((-1,)).names() == ctx1.names()[:-1]
+
+    assert ctx.infer_generator_mapping(ctx) == {i: i for i in range(ctx.nvars())}
+    assert ctx1.infer_generator_mapping(ctx) == {0: 0, 1: 1}
+    assert ctx1.drop_gens(ctx.names()).infer_generator_mapping(ctx) == {}
+
+    assert quick_poly().project_to_context(ctx1) == \
+        ctx1.from_dict(
+            {(0, 0, 0, 0): 1, (0, 1, 0, 0): 2, (1, 0, 0, 0): 3, (2, 2, 0, 0): 4}
+        )
+    new_poly = quick_poly().project_to_context(ctx1)
+    assert ctx1.drop_gens(new_poly.unused_gens()) == ctx
+    assert new_poly.project_to_context(ctx) == quick_poly()
+
+    new_poly = quick_poly().project_to_context(ctx2)
+    new_ctx = ctx2.drop_gens(new_poly.unused_gens())
+    assert new_ctx != ctx
+    assert new_poly != quick_poly()
+
+    new_ctx = new_ctx.from_context(new_ctx, ordering=ctx.ordering())
+    assert new_ctx == ctx
+    assert new_poly.project_to_context(new_ctx) == quick_poly()
+
+    assert ctx.append_gens(*ctx1.names()[-2:]) == ctx1
+
+    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() / 2
+    assert ctx.gen(1) == ctx.from_dict({(0, 1): 1})
+
+    assert ctx.nvars() == 2
+    assert ctx.ordering() == flint.Ordering.lex
+
+    ctx1 = get_context(("x", 4))
+    assert [ctx1.name(i) for i in range(4)] == ['x0', 'x1', 'x2', 'x3']
+    for order in list(flint.Ordering):
+        ctx1 = get_context(("x", 4), ordering=order)
+        assert ctx1.ordering() == order
+
+    assert ctx.constant(1) == mpoly({(0, 0): 1}) == P(1, ctx=ctx)
+
+    assert raises(lambda: P([None]), TypeError) # type: ignore
+    assert raises(lambda: P(object()), TypeError) # type: ignore
+    assert raises(lambda: P(None), TypeError) # type: ignore
+    assert raises(lambda: P(None, None), TypeError) # type: ignore
+    assert raises(lambda: P([1,2], None), TypeError) # type: ignore
+    assert raises(lambda: P(1, None), ValueError) # type: ignore
+
+    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 True
+    assert (P(1, ctx=ctx) != 1) is False
+    assert (1 == P(1, ctx=ctx)) is True
+    assert (1 != P(1, ctx=ctx)) is False
+
+    assert (P(1, ctx=ctx) == S(1)) is True
+    assert (P(1, ctx=ctx) != S(1)) is False
+    assert (S(1) == P(1, ctx=ctx)) is True
+    assert (S(1) != P(1, ctx=ctx)) is False
+
+    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 and isinstance(ctx, flint.fmpq_mpoly_ctx):
+        ctx_z = flint.fmpz_mpoly_ctx.get((("x", 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)
+
+
+@all_mpolys
+def test_mpolys_properties(args: _MPolyTestCase[Tmpoly, Tscalar]) -> None:
+
+    P, get_context, S, is_field, characteristic = args
+
+    # Division under modulo will raise a flint exception if something is
+    # not invertible, crashing the program. We can't tell before what is
+    # invertible and what is not before hand so we always raise an
+    # exception, except for fmpz_mpoly, that returns an bool noting if the
+    # division is exact or not.
+    composite_characteristic = characteristic != 0 and not characteristic.is_prime()
+
+    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})
+
+    ctx = get_context(("x", 2))
+    ctx1 = get_context(("x", 4))
+
+    assert raises(lambda: P(ctx=ctx) < P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) <= P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) > P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) >= P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) < None, TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) <= None, TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) > None, TypeError) # type: ignore
+    assert raises(lambda: P(ctx=ctx) >= None, TypeError) # type: ignore
+    assert raises(lambda: None < P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: None <= P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: None > P(ctx=ctx), TypeError) # type: ignore
+    assert raises(lambda: None >= P(ctx=ctx), TypeError) # type: ignore
+
+    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) # type: ignore
+    assert raises(lambda: p[4], TypeError) # type: ignore
+
+    assert p[(2, 2)] == 4
+    assert p[(0, 0)] == 1
+    assert raises(lambda: p[(1,)], ValueError)
+    assert raises(lambda: p[(1, "bad")], TypeError) # type: ignore
+    assert raises(lambda: p["bad"], TypeError) # type: ignore
+
+    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) # type: ignore
+    assert raises(lambda: p.__setitem__(("bad", 1), 1), TypeError) # type: ignore
+
+    assert raises(lambda: p.__setitem__((2, 1), None), TypeError) # type: ignore
+
+    assert P(ctx=ctx).repr() == f"{ctx.__class__.__name__}(2, '<Ordering.lex: 'lex'>', ('x0', 'x1')).from_dict({{}})"
+    assert P(1, ctx=ctx).repr() == f"{ctx.__class__.__name__}(2, '<Ordering.lex: 'lex'>', ('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) # type: ignore
+    assert raises(lambda: p(1), ValueError) # type: ignore
+    assert raises(lambda: p(0, 1, 2), ValueError) # type: ignore
+
+    assert raises(lambda: p.subs({"x0": None}), TypeError) # type: ignore
+    assert raises(lambda: p.subs({"x0": None, "x1": None}), TypeError) # type: ignore
+    assert raises(lambda: p.subs({"a": 1}), ValueError)
+    assert raises(lambda: p.subs({"x0": 0, "x1": 1, "x2": 2}), ValueError)
+
+    no_gens_ctx = get_context(tuple())
+    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[(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) # type: ignore
-        assert raises(lambda: p.__setitem__(("bad", 1), 1), TypeError) # type: ignore
-
-        assert raises(lambda: p.__setitem__((2, 1), None), TypeError) # type: ignore
-
-        assert P(ctx=ctx).repr() == f"{ctx.__class__.__name__}(2, '<Ordering.lex: 'lex'>', ('x0', 'x1')).from_dict({{}})"
-        assert P(1, ctx=ctx).repr() == f"{ctx.__class__.__name__}(2, '<Ordering.lex: 'lex'>', ('x0', 'x1')).from_dict({{(0, 0): 1}})"
-        assert str(quick_poly()) == repr(quick_poly()) == '4*x0^2*x1^2 + 3*x0 + 2*x1 + 1'
+        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) # type: ignore
+    assert raises(lambda: None + mpoly({(0, 0): 2, (0, 1): 2, (1, 0): 3, (2, 2): 4}), TypeError) # type: ignore
+    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) # type: ignore
+
+    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) # type: ignore
+    assert raises(lambda: None - quick_poly(), TypeError) # type: ignore
+    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) # type: ignore
+
+    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
+        })
 
-        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)
+    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) # type: ignore
+    assert raises(lambda: None * quick_poly(), TypeError) # type: ignore
+    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) # type: ignore
+
+    if composite_characteristic:
+        assert raises(lambda: quick_poly() // mpoly({(1, 1): 1}), DomainError)
+        assert raises(lambda: quick_poly() % mpoly({(1, 1): 1}), DomainError)
+        assert raises(lambda: divmod(quick_poly(), mpoly({(1, 1): 1})), DomainError)
+    else:
+        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))
+
+        assert S(1) / P(1, ctx=ctx) == P(1, ctx=ctx)
+        assert quick_poly() / S(1) == quick_poly()
+        assert quick_poly() // S(1) == quick_poly()
+        assert quick_poly() % S(1) == P(ctx=ctx)
+        assert divmod(quick_poly(), S(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 p.total_degree() == 4
-        assert P(ctx=ctx).total_degree() == -1
-        assert P(1, ctx=ctx).total_degree() == 0
+    f = mpoly({(1, 1): 4, (0, 0): 1})
+    g = mpoly({(0, 1): 2, (1, 0): 2})
+    if not composite_characteristic:
+        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))
 
-        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
+        assert S(1) // quick_poly() == P(ctx=ctx)
+        assert S(1) % quick_poly() == P(1, ctx=ctx)
+        assert divmod(S(1), quick_poly()) == (P(ctx=ctx), P(1, ctx=ctx))
 
-        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) # type: ignore
-        assert raises(lambda: p(1), ValueError) # type: ignore
-        assert raises(lambda: p(0, 1, 2), ValueError) # type: ignore
-
-        assert raises(lambda: p.subs({"x0": None}), TypeError) # type: ignore
-        assert raises(lambda: p.subs({"x0": None, "x1": None}), TypeError) # type: ignore
-        assert raises(lambda: p.subs({"a": 1}), ValueError)
-        assert raises(lambda: p.subs({"x0": 0, "x1": 1, "x2": 2}), ValueError)
-
-        no_gens_ctx = get_context(tuple())
-        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) # type: ignore
-        assert raises(lambda: None + mpoly({(0, 0): 2, (0, 1): 2, (1, 0): 3, (2, 2): 4}), TypeError) # type: ignore
-        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) # type: ignore
-
-        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) # type: ignore
-        assert raises(lambda: None - quick_poly(), TypeError) # type: ignore
-        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) # type: ignore
-
-        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) # type: ignore
-        assert raises(lambda: None * quick_poly(), TypeError) # type: ignore
-        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) # type: ignore
-
-        if composite_characteristic:
-            assert raises(lambda: quick_poly() // mpoly({(1, 1): 1}), DomainError)
-            assert raises(lambda: quick_poly() % mpoly({(1, 1): 1}), DomainError)
-            assert raises(lambda: divmod(quick_poly(), mpoly({(1, 1): 1})), DomainError)
-        else:
-            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))
-
-            assert S(1) / P(1, ctx=ctx) == P(1, ctx=ctx)
-            assert quick_poly() / S(1) == quick_poly()
-            assert quick_poly() // S(1) == quick_poly()
-            assert quick_poly() % S(1) == P(ctx=ctx)
-            assert divmod(quick_poly(), S(1)) == (quick_poly(), P(ctx=ctx))
+        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)
+
+    # We prefer various other errors to the "division not supported" domain error so these are safe.
+    assert raises(lambda: quick_poly() / None, TypeError) # type: ignore
+    assert raises(lambda: quick_poly() // None, TypeError) # type: ignore
+    assert raises(lambda: quick_poly() % None, TypeError) # type: ignore
+    assert raises(lambda: divmod(quick_poly(), None), TypeError) # type: ignore
+
+    assert raises(lambda: None / quick_poly(), TypeError) # type: ignore
+    assert raises(lambda: None // quick_poly(), TypeError) # type: ignore
+    assert raises(lambda: None % quick_poly(), TypeError) # type: ignore
+    assert raises(lambda: divmod(None, quick_poly()), TypeError) # type: ignore
+
+    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)
+
+    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, ZeroDivisionError)
+    assert raises(lambda: P(ctx=ctx) ** None, TypeError) # type: ignore
+
+    # # XXX: Not sure what this should do in general:
+    assert raises(lambda: pow(P(1, ctx=ctx), 2, 3), NotImplementedError) # type: ignore
+
+    if composite_characteristic:
+        assert raises(lambda: (f * g).gcd(f), DomainError)
+    else:
         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})
+            assert (f * g).gcd(f) == f / 4
         else:
-            assert raises(lambda: quick_poly() / 3, DomainError)
+            assert (f * g).gcd(f) == f
+        assert raises(lambda: quick_poly().gcd(None), TypeError) # type: ignore
+        assert raises(lambda: quick_poly().gcd(P(ctx=ctx1)), IncompatibleContextError)
 
-        f = mpoly({(1, 1): 4, (0, 0): 1})
-        g = mpoly({(0, 1): 2, (1, 0): 2})
-        if not composite_characteristic:
-            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))
+    x0, x1 = ctx.gens()
+    assert (2*x0**3 + 4*x0**2 + 6*x0).term_content() == x0 if is_field else 2*x0
+    assert (3*x0**2*x1 + 6*x0*x1**2 + 9*x1**3).term_content() == x1 if is_field else 3*x1
 
-            assert S(1) // quick_poly() == P(ctx=ctx)
-            assert S(1) % quick_poly() == P(1, ctx=ctx)
-            assert divmod(S(1), quick_poly()) == (P(ctx=ctx), P(1, ctx=ctx))
+    assert (x0**2 - 1).resultant(x0**2 - 2*x0 + 1, 'x0') == 0
+    assert (x0**2 + x1**2 - 1).resultant(x0 - x1, 'x0') == 2*x1**2 - 1
 
-            assert f * g / mpoly({(0, 1): 1, (1, 0): 1}) \
-                == mpoly({(1, 1): 8, (0, 0): 2})
+    assert (x0**3 - 6*x0**2 + 11*x0 - 6).discriminant('x0') == 4
+    assert (x0**2 + 4*x0*x1 + 4*x1**2 - 1).discriminant('x0') == 4
 
-            if not is_field:
-                assert raises(lambda: 1 / quick_poly(), DomainError)
-                assert raises(lambda: quick_poly() / P(2, ctx=ctx), DomainError)
-
-        # We prefer various other errors to the "division not supported" domain error so these are safe.
-        assert raises(lambda: quick_poly() / None, TypeError) # type: ignore
-        assert raises(lambda: quick_poly() // None, TypeError) # type: ignore
-        assert raises(lambda: quick_poly() % None, TypeError) # type: ignore
-        assert raises(lambda: divmod(quick_poly(), None), TypeError) # type: ignore
-
-        assert raises(lambda: None / quick_poly(), TypeError) # type: ignore
-        assert raises(lambda: None // quick_poly(), TypeError) # type: ignore
-        assert raises(lambda: None % quick_poly(), TypeError) # type: ignore
-        assert raises(lambda: divmod(None, quick_poly()), TypeError) # type: ignore
-
-        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)
-
-        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, ZeroDivisionError)
-        assert raises(lambda: P(ctx=ctx) ** None, TypeError) # type: ignore
+    f1 = 3*x0**2*x1**2 + 6*x0*x1**2 + 9*x1**2
+    res, stride = f1.deflation()
+    assert res == 3*x0**2*x1 + 6*x0*x1 + 9*x1
+    assert tuple(stride) == (1, 2)
 
-        # # XXX: Not sure what this should do in general:
-        assert raises(lambda: pow(P(1, ctx=ctx), 2, 3), NotImplementedError) # type: ignore
+    g1 = ((x0**2 + x1**2)**3 + (x0**2 + x1**2)**2 + 1)
+    res, stride = g1.deflation()
+    assert res == x0**3 + 3*x0**2*x1 + x0**2 + 3*x0*x1**2 + 2*x0*x1 + x1**3 + x1**2 + 1
+    assert tuple(stride) == (2, 2)
 
-        if composite_characteristic:
-            assert raises(lambda: (f * g).gcd(f), DomainError)
-        else:
-            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) # type: ignore
-            assert raises(lambda: quick_poly().gcd(P(ctx=ctx1)), IncompatibleContextError)
-
-        x0, x1 = ctx.gens()
-        assert (2*x0**3 + 4*x0**2 + 6*x0).term_content() == x0 if is_field else 2*x0
-        assert (3*x0**2*x1 + 6*x0*x1**2 + 9*x1**3).term_content() == x1 if is_field else 3*x1
-
-        assert (x0**2 - 1).resultant(x0**2 - 2*x0 + 1, 'x0') == 0
-        assert (x0**2 + x1**2 - 1).resultant(x0 - x1, 'x0') == 2*x1**2 - 1
-
-        assert (x0**3 - 6*x0**2 + 11*x0 - 6).discriminant('x0') == 4
-        assert (x0**2 + 4*x0*x1 + 4*x1**2 - 1).discriminant('x0') == 4
-
-        f1 = 3*x0**2*x1**2 + 6*x0*x1**2 + 9*x1**2
-        res, stride = f1.deflation()
-        assert res == 3*x0**2*x1 + 6*x0*x1 + 9*x1
-        assert tuple(stride) == (1, 2)
-
-        g1 = ((x0**2 + x1**2)**3 + (x0**2 + x1**2)**2 + 1)
-        res, stride = g1.deflation()
-        assert res == x0**3 + 3*x0**2*x1 + x0**2 + 3*x0*x1**2 + 2*x0*x1 + x1**3 + x1**2 + 1
-        assert tuple(stride) == (2, 2)
-
-        for p in [f1, g1]:
-            pd, n = p.deflation()
-            assert pd.inflate(n) == p
-            assert p.deflate(n).inflate(n) == p
-
-            pd, n, m = p.deflation_monom()
-            assert m * pd.inflate(n) == p
-
-            if not composite_characteristic:
-                n, i = p.deflation_index()
-                m = ctx.term(exp_vec=i)
-                assert (p / m).deflate(n).inflate(n) * m == p
-
-        if P is flint.fmpz_mpoly:
-            assert (x0**2 * x1 + x0 * x1).primitive() == (1, x0**2*x1 + x0*x1)
-            assert (4 * x0 + 2 * x0 * x1).primitive() == (2, x0 * x1 + 2 * x0)
-
-        if composite_characteristic:
-            # Factorisation not allowed over Z/nZ for n not prime.
-            # Flint would abort so we raise an exception instead:
-            assert raises(lambda: (f * g).factor(), DomainError)
-        elif characteristic == 0:
-            # Primitive factors over Z for fmpz_mpoly and fmpq_mpoly
-            assert (f * g).factor() == (S(2), [(g / 2, 1), (f, 1)])
-        elif is_field:
-            # Monic polynomials over Z/pZ for nmod_mpoly and fmpz_mod_mpoly
-            assert (f * g).factor() == (S(8), [(g / 2, 1), (f / 4, 1)])
-
-        if composite_characteristic:
-            assert raises(lambda: (f * g).sqrt(), DomainError)
-        else:
-            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(), DomainError)
+    for p in [f1, g1]:
+        pd, n = p.deflation()
+        assert pd.inflate(n) == p
+        assert p.deflate(n).inflate(n) == p
 
-        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})
+        pd, n, m = p.deflation_monom()
+        assert m * pd.inflate(n) == p
 
-        assert raises(lambda: p.derivative("x3"), ValueError)
-        assert raises(lambda: p.derivative(3), IndexError)
-        assert raises(lambda: p.derivative(None), TypeError) # type: ignore
+        if not composite_characteristic:
+            n, i = p.deflation_index()
+            m = ctx.term(exp_vec=i)
+            assert (p / m).deflate(n).inflate(n) * m == p
+
+    if isinstance(x0, flint.fmpz_mpoly) and isinstance(x1, flint.fmpz_mpoly):
+        assert (x0**2 * x1 + x0 * x1).primitive() == (1, x0**2*x1 + x0*x1)
+        assert (4 * x0 + 2 * x0 * x1).primitive() == (2, x0 * x1 + 2 * x0)
+
+    if composite_characteristic:
+        # Factorisation not allowed over Z/nZ for n not prime.
+        # Flint would abort so we raise an exception instead:
+        assert raises(lambda: (f * g).factor(), DomainError)
+    elif characteristic == 0:
+        # Primitive factors over Z for fmpz_mpoly and fmpq_mpoly
+        assert (f * g).factor() == (S(2), [(g / 2, 1), (f, 1)])
+    elif is_field:
+        # Monic polynomials over Z/pZ for nmod_mpoly and fmpz_mod_mpoly
+        assert (f * g).factor() == (S(8), [(g / 2, 1), (f / 4, 1)])
+
+    if composite_characteristic:
+        assert raises(lambda: (f * g).sqrt(), DomainError)
+    else:
+        assert (f * f).sqrt() == f
+        if isinstance(f, flint.fmpz_mpoly):
+            assert (f * f).sqrt(assume_perfect_square=True) == f
+        assert raises(lambda: quick_poly().sqrt(), DomainError)
+
+    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) # type: ignore
+
+    if isinstance(p, (flint.fmpz_mpoly, flint.fmpq_mpoly)):
+        if isinstance(p, flint.fmpq_mpoly) and _is_Q(S):
+            assert p.integral(0) == p.integral("x0") == \
+                mpoly({(3, 2): S(4, 3), (2, 0): S(3, 2), (1, 1): S(2), (1, 0): S(1)})
+            assert p.integral(1) == p.integral("x1") == \
+                mpoly({(2, 3): S(4, 3), (1, 1): S(3), (0, 2): S(1), (0, 1): S(1)})
+        else:
+            assert p.integral(0) == p.integral("x0") == \
+                (6, mpoly({(3, 2): 8, (2, 0): 9, (1, 1): 12, (1, 0): 6}))
+            assert p.integral(1) == p.integral("x1") == \
+                (3, mpoly({(2, 3): 4, (1, 1): 9, (0, 2): 3, (0, 1): 3}))
 
-        if (P is not flint.fmpz_mod_mpoly and P is not flint.nmod_mpoly):
-            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) # type: ignore
 
-            assert raises(lambda: p.integral("x3"), ValueError)
-            assert raises(lambda: p.integral(3), IndexError)
-            assert raises(lambda: p.integral(None), TypeError) # type: ignore
+    elif isinstance(p, (flint.nmod_mpoly, flint.fmpz_mod_mpoly)):
+        pass # No integral method
+    else:
+        assert False # New _mpoly types?
 
 
 def _all_mpoly_vecs():
@@ -3812,7 +3909,7 @@ def test_division_poly_mpoly():
             assert two ** -1 == two ** Z(-1) == half
         elif characteristic == 0:
             assert raises(lambda: x / 2, DomainError)
-            assert raises(lambda: x / two, DomainError), characteristic
+            assert raises(lambda: x / two, DomainError)
             assert raises(lambda: two ** -1, DomainError)
             assert raises(lambda: two ** Z(-1), DomainError)
         else:
@@ -4830,7 +4927,7 @@ def test_python_threads():
     from threading import Thread
 
     iterations = 10**5
-    threads = 3 + 1
+    nthreads = 3 + 1
     size = 3
     M = flint.fmpz_mat([[0]*size for _ in range(size)])
 
@@ -4849,7 +4946,7 @@ def get_dets():
         for _ in range(iterations):
             M.det()
 
-    threads = [Thread(target=set_values) for _ in range(threads-1)]
+    threads = [Thread(target=set_values) for _ in range(nthreads-1)]
     threads.append(Thread(target=get_dets))
 
     for t in threads:
@@ -4902,7 +4999,8 @@ def test_all_tests():
     test_division_poly_mpoly,
 
     test_polys,
-    test_mpolys,
+    test_mpolys_constructor,
+    test_mpolys_properties,
 
     test_poly_resultants,
 
diff --git a/src/flint/types/fmpq.pyi b/src/flint/types/fmpq.pyi
index 47a4fb45..9232f93a 100644
--- a/src/flint/types/fmpq.pyi
+++ b/src/flint/types/fmpq.pyi
@@ -3,19 +3,18 @@ from typing import overload
 from .fmpz import fmpz
 
 
+_str = str
 ifmpz = int | fmpz
 ifmpq = int | fmpz | fmpq
 
 
 class fmpq(flint_scalar):
     @overload
-    def __init__(self): ...
+    def __init__(self) -> None: ...
     @overload
-    def __init__(self, arg: ifmpq | str, /): ...
+    def __init__(self, arg: ifmpq | _str, /) -> None: ...
     @overload
-    def __init__(self, num: ifmpz, den: ifmpz, /): ...
-
-    def __init__(self, arg1: ifmpq | str = 0, arg2: ifmpz = 1, /): ...
+    def __init__(self, num: ifmpz, den: ifmpz, /) -> None: ...
 
     @property
     def p(self) -> fmpz: ...
@@ -34,11 +33,11 @@ class fmpq(flint_scalar):
 
     def next(self, signed: bool = True, minimal: bool = True) -> fmpq: ...
 
-    def str(self, base: int = 10, condense: int = 0) -> str: ...
-    def repr(self) -> str: ...
+    def str(self, base: int = 10, condense: int = 0) -> _str: ...
+    def repr(self) -> _str: ...
 
-    def __str__(self) -> str: ...
-    def __repr__(self) -> str: ...
+    def __str__(self) -> _str: ...
+    def __repr__(self) -> _str: ...
 
     def __int__(self) -> int: ...
 
@@ -54,7 +53,8 @@ class fmpq(flint_scalar):
     @overload
     def __round__(self) -> fmpz: ...
 
-    def __round__(self, ndigits: int | None = None) -> fmpq | fmpz: ...
+    def __eq__(self, other: object, /) -> bool: ...
+    def __hash__(self) -> int: ...
 
     def __lt__(self, other: ifmpq, /) -> bool: ...
     def __le__(self, other: ifmpq, /) -> bool: ...
diff --git a/src/flint/types/fmpq_mpoly.pyi b/src/flint/types/fmpq_mpoly.pyi
index 99c0857d..833ba9c7 100644
--- a/src/flint/types/fmpq_mpoly.pyi
+++ b/src/flint/types/fmpq_mpoly.pyi
@@ -5,6 +5,7 @@ from .fmpq import fmpq
 from .fmpz_mpoly import fmpz_mpoly
 
 
+_str = str
 ifmpz = int | fmpz
 ifmpq = int | fmpz | fmpq
 
@@ -12,8 +13,8 @@ ifmpq = int | fmpz | fmpq
 class fmpq_mpoly_ctx(flint_mpoly_context[fmpq_mpoly, fmpq, ifmpq]):
     @classmethod
     def get(cls,
-        names: str | Iterable[str | tuple[str, int]] | tuple[str, int],
-        ordering: Ordering | str = Ordering.lex
+        names: _str | Iterable[_str | tuple[_str, int]] | tuple[_str, int],
+        ordering: Ordering | _str = Ordering.lex
     ) -> fmpq_mpoly_ctx:
         ...
 
@@ -27,12 +28,12 @@ class fmpq_mpoly_ctx(flint_mpoly_context[fmpq_mpoly, fmpq, ifmpq]):
 
 class fmpq_mpoly(flint_mpoly[fmpq_mpoly_ctx, fmpq, ifmpq]):
     def __init__(self,
-        val: fmpq_mpoly | fmpz_mpoly | ifmpq | dict[tuple[int, ...], ifmpq] | str = 0,
+        val: fmpq_mpoly | fmpz_mpoly | ifmpq | dict[tuple[int, ...], ifmpq] | _str = 0,
         ctx: fmpq_mpoly_ctx | None = None
     ) -> None: ...
 
-    def str(self) -> str: ...
-    def repr(self) -> str: ...
+    def str(self) -> _str: ...
+    def repr(self) -> _str: ...
 
     def context(self) -> fmpq_mpoly_ctx: ...
 
@@ -52,7 +53,7 @@ class fmpq_mpoly(flint_mpoly[fmpq_mpoly_ctx, fmpq, ifmpq]):
     def __getitem__(self, index: tuple[int, ...]) -> fmpq: ...
     def __setitem__(self, index: tuple[int, ...], coeff: ifmpq) -> None: ...
 
-    def subs(self, mapping: dict[str | int, ifmpq]) -> fmpq_mpoly: ...
+    def subs(self, mapping: dict[_str | int, ifmpq]) -> fmpq_mpoly: ...
     def compose(self, *args: fmpq_mpoly, ctx: fmpq_mpoly_ctx | None = None) -> fmpq_mpoly: ...
 
     def __call__(self, *args: ifmpq) -> fmpq: ...
@@ -89,8 +90,8 @@ class fmpq_mpoly(flint_mpoly[fmpq_mpoly_ctx, fmpq, ifmpq]):
 
     def sqrt(self, assume_perfect_square: bool = False) -> fmpq_mpoly: ...
 
-    def resultant(self, other: fmpq_mpoly, var: str | int) -> fmpq_mpoly: ...
-    def discriminant(self, var: str | int) -> fmpq_mpoly: ...
+    def resultant(self, other: fmpq_mpoly, var: _str | int) -> fmpq_mpoly: ...
+    def discriminant(self, var: _str | int) -> fmpq_mpoly: ...
 
     def deflation(self) -> tuple[fmpq_mpoly, list[int]]: ...
     def inflate(self, N: list[int]) -> fmpq_mpoly: ...
@@ -98,8 +99,8 @@ class fmpq_mpoly(flint_mpoly[fmpq_mpoly_ctx, fmpq, ifmpq]):
     def deflation_monom(self) -> tuple[fmpq_mpoly, list[int], fmpq_mpoly]: ...
     def deflation_index(self) -> tuple[list[int], list[int]]: ...
 
-    def derivative(self, var: str | int) -> fmpq_mpoly: ...
-    def integral(self, var: str | int) -> fmpq_mpoly: ...
+    def derivative(self, var: _str | int) -> fmpq_mpoly: ...
+    def integral(self, var: _str | int) -> fmpq_mpoly: ...
 
 
 class fmpq_mpoly_vec:
diff --git a/src/flint/types/fmpq_poly.pyi b/src/flint/types/fmpq_poly.pyi
new file mode 100644
index 00000000..8ff1d09d
--- /dev/null
+++ b/src/flint/types/fmpq_poly.pyi
@@ -0,0 +1,94 @@
+from typing import overload, Any
+from flint.flint_base.flint_base import flint_poly
+from flint.types.fmpz import fmpz, ifmpz
+from flint.types.fmpq import fmpq, ifmpq
+from flint.types.fmpz_poly import fmpz_poly, ifmpz_poly
+
+
+ifmpq_poly = fmpq_poly | ifmpq | ifmpz_poly
+
+
+class fmpq_poly(flint_poly[fmpq]):
+    """
+    The *fmpq_poly* type represents dense univariate polynomials
+    over the rational numbers. For efficiency reasons, an *fmpq_poly* is
+    structurally an integer polynomial with a single common denominator.
+    """
+
+    @overload
+    def __init__(self) -> None: ...
+    @overload
+    def __init__(self, other: fmpq_poly | fmpz_poly | list[int] | list[ifmpq] | ifmpq, /) -> None: ...
+    @overload
+    def __init__(self, other: fmpq_poly | fmpz_poly | list[int] | list[ifmpq] | ifmpq,
+                       den: ifmpz, /) -> None: ...
+
+    def __len__(self) -> int: ...
+    def length(self) -> int: ...
+    def degree(self) -> int: ...
+
+    def numer(self) -> fmpz_poly: ...
+    def denom(self) -> fmpz: ...
+
+    @property
+    def p(self) -> fmpz_poly: ...
+    @property
+    def q(self) -> fmpz: ...
+
+    def __getitem__(self, i: int, /) -> fmpq: ...
+    def __setitem__(self, i: int, x: ifmpq, /) -> None: ...
+
+    def repr(self) -> str: ...
+
+    def __bool__(self) -> bool: ...
+    def is_zero(self) -> bool: ...
+    def is_one(self) -> bool: ...
+    def is_constant(self) -> bool: ...
+
+    def leading_coefficient(self) -> fmpq: ...
+
+    @overload
+    def __call__(self, other: ifmpq, /) -> fmpq: ...
+    @overload
+    def __call__(self, other: fmpq_poly | fmpz_poly, /) -> fmpq_poly: ...
+
+    def derivative(self) -> fmpq_poly: ...
+    def integral(self) -> fmpq_poly: ...
+
+    def __pos__(self) -> fmpq_poly: ...
+    def __neg__(self) -> fmpq_poly: ...
+    def __add__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __radd__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __sub__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __rsub__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __mul__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __rmul__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __floordiv__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __rfloordiv__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __mod__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __rmod__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __truediv__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __rtruediv__(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def __divmod__(self, other: ifmpq_poly, /) -> tuple[fmpq_poly, fmpq_poly]: ...
+    def __rdivmod__(self, other: ifmpq_poly, /) -> tuple[fmpq_poly, fmpq_poly]: ...
+    def __pow__(self, exp: int | ifmpz, /) -> fmpq_poly: ...
+
+    def gcd(self, other: ifmpq_poly, /) -> fmpq_poly: ...
+    def resultant(self, other: ifmpq_poly, /) -> fmpq: ...
+    def xgcd(self, other: ifmpq_poly, /) -> tuple[fmpq_poly, fmpq_poly, fmpq_poly]: ...
+
+    def factor(self, *, monic: bool = False) -> tuple[fmpq, list[tuple[fmpq_poly, int]]]: ...
+    def factor_squarefree(self) -> tuple[fmpq, list[tuple[fmpq_poly, int]]]: ...
+
+    def sqrt(self) -> fmpq_poly: ...
+
+    def deflation(self) -> tuple[fmpq_poly, int]: ...
+    # XXX: acb
+    def complex_roots(self, **kwargs) -> list[Any]: ...
+
+    @staticmethod
+    def bernoulli_poly(n: int, /) -> fmpq_poly: ...
+    @staticmethod
+    def euler_poly(n: int, /) -> fmpq_poly: ...
+    @staticmethod
+    def legendre_p(n: int, /) -> fmpq_poly: ...
diff --git a/src/flint/types/fmpq_series.pyi b/src/flint/types/fmpq_series.pyi
new file mode 100644
index 00000000..678c81b9
--- /dev/null
+++ b/src/flint/types/fmpq_series.pyi
@@ -0,0 +1,78 @@
+from typing import Any
+from flint.flint_base.flint_base import flint_series
+from flint.types.fmpz import fmpz, ifmpz
+from flint.types.fmpq import fmpq, ifmpq
+from flint.types.fmpq_poly import fmpq_poly
+from flint.types.fmpz_series import fmpz_series, ifmpz_series
+
+
+ifmpq_series = fmpq_series | fmpq_poly | ifmpz_series
+
+
+class fmpq_series(flint_series[fmpq]):
+    """Approximate truncated power series with rational coefficients."""
+
+    def __init__(self,
+                 val: list[int] | list[ifmpq] | ifmpq | ifmpq_series | None = None,
+                 den: ifmpz | None = None,
+                 prec: int | None = None): ...
+
+    @property
+    def prec(self) -> int: ...
+
+    def _equal_repr(self, other: object) -> bool: ...
+
+    def __len__(self) -> int: ...
+    def length(self) -> int: ...
+
+    def numer(self) -> fmpz_series: ...
+    def denom(self) -> fmpz: ...
+
+    @property
+    def p(self) -> fmpq_poly: ...
+    @property
+    def q(self) -> fmpz: ...
+
+    def __getitem__(self, i: int, /) -> fmpq: ...
+    def __setitem__(self, i: int, x: ifmpq, /) -> None: ...
+
+    def repr(self, **kwargs: Any) -> str: ...
+    def str(self, **kwargs: Any) -> str: ...
+
+    def __pos__(self) -> fmpq_series: ...
+    def __neg__(self) -> fmpq_series: ...
+    def __add__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __radd__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __sub__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __rsub__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __mul__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __rmul__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __truediv__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __rtruediv__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def __pow__(self, exp: int, /) -> fmpq_series: ...
+
+    def valuation(self) -> int: ...
+
+    def __call__(self, other: ifmpq_series, /) -> fmpq_series: ...
+    def reversion(self) -> fmpq_series: ...
+    def inv(self) -> fmpq_series: ...
+
+    def zero_constant_term(self) -> bool: ...
+    def one_constant_term(self) -> bool: ...
+
+    def derivative(self) -> fmpq_series: ...
+    def integral(self) -> fmpq_series: ...
+    def sqrt(self) -> fmpq_series: ...
+    def rsqrt(self) -> fmpq_series: ...
+    def exp(self) -> fmpq_series: ...
+    def log(self) -> fmpq_series: ...
+    def atan(self) -> fmpq_series: ...
+    def atanh(self) -> fmpq_series: ...
+    def asin(self) -> fmpq_series: ...
+    def asinh(self) -> fmpq_series: ...
+    def sin(self) -> fmpq_series: ...
+    def cos(self) -> fmpq_series: ...
+    def tan(self) -> fmpq_series: ...
+    def sinh(self) -> fmpq_series: ...
+    def cosh(self) -> fmpq_series: ...
+    def tanh(self) -> fmpq_series: ...
diff --git a/src/flint/types/fmpz.pyi b/src/flint/types/fmpz.pyi
index 2d4b76a5..a94df7b8 100644
--- a/src/flint/types/fmpz.pyi
+++ b/src/flint/types/fmpz.pyi
@@ -1,11 +1,14 @@
 #
 # It would be better if many methods accepted SupportsIndex rather than just
-# int | fmpz but for now these annotations are accurate.
+# ifmpz but for now these annotations are accurate.
 #
 from ..flint_base.flint_base import flint_scalar
 
+_str = str
+ifmpz = int | fmpz
+
 class fmpz(flint_scalar):
-    def __init__(self, arg: int | fmpz | str = 0, /): ...
+    def __init__(self, arg: ifmpz | _str = 0, /): ...
 
     @property
     def numerator(self) -> fmpz: ...
@@ -15,11 +18,11 @@ class fmpz(flint_scalar):
     def bit_length(self) -> int: ...
     def height_bits(self, signed: bool = False) -> int: ...
 
-    def str(self, base: int = 10, condense: int = 0) -> str: ...
-    def repr(self) -> str: ...
+    def str(self, base: int = 10, condense: int = 0) -> _str: ...
+    def repr(self) -> _str: ...
 
-    def __str__(self) -> str: ...
-    def __repr__(self) -> str: ...
+    def __str__(self) -> _str: ...
+    def __repr__(self) -> _str: ...
 
     def __int__(self) -> int: ...
     def __index__(self) -> int: ...
@@ -29,47 +32,50 @@ class fmpz(flint_scalar):
     def __trunc__(self) -> fmpz: ...
     def __round__(self, ndigits: int | None = None) -> fmpz: ...
 
-    def __lt__(self, other: int | fmpz, /) -> bool: ...
-    def __le__(self, other: int | fmpz, /) -> bool: ...
-    def __gt__(self, other: int | fmpz, /) -> bool: ...
-    def __ge__(self, other: int | fmpz, /) -> bool: ...
+    def __eq__(self, other: object, /) -> bool: ...
+    def __hash__(self) -> int: ...
+
+    def __lt__(self, other: ifmpz, /) -> bool: ...
+    def __le__(self, other: ifmpz, /) -> bool: ...
+    def __gt__(self, other: ifmpz, /) -> bool: ...
+    def __ge__(self, other: ifmpz, /) -> bool: ...
 
     def __neg__(self) -> fmpz: ...
     def __pos__(self) -> fmpz: ...
-    def __add__(self, other: int | fmpz, /) -> fmpz: ...
+    def __add__(self, other: ifmpz, /) -> fmpz: ...
     def __radd__(self, other: int, /) -> fmpz: ...
-    def __sub__(self, other: int | fmpz, /) -> fmpz: ...
+    def __sub__(self, other: ifmpz, /) -> fmpz: ...
     def __rsub__(self, other: int, /) -> fmpz: ...
-    def __mul__(self, other: int | fmpz, /) -> fmpz: ...
+    def __mul__(self, other: ifmpz, /) -> fmpz: ...
     def __rmul__(self, other: int, /) -> fmpz: ...
-    def __truediv__(self, other: int | fmpz, /) -> fmpz: ...
+    def __truediv__(self, other: ifmpz, /) -> fmpz: ...
     def __rtruediv__(self, other: int, /) -> fmpz: ...
-    def __floordiv__(self, other: int | fmpz, /) -> fmpz: ...
+    def __floordiv__(self, other: ifmpz, /) -> fmpz: ...
     def __rfloordiv__(self, other: int, /) -> fmpz: ...
-    def __mod__(self, other: int | fmpz, /) -> fmpz: ...
+    def __mod__(self, other: ifmpz, /) -> fmpz: ...
     def __rmod__(self, other: int, /) -> fmpz: ...
-    def __divmod__(self, other: int | fmpz, /) -> tuple[fmpz, fmpz]: ...
+    def __divmod__(self, other: ifmpz, /) -> tuple[fmpz, fmpz]: ...
     def __rdivmod__(self, other: int, /) -> tuple[fmpz, fmpz]: ...
-    def __pow__(self, other: int | fmpz, modulo: int | fmpz | None = None, /) -> fmpz: ...
-    def __rpow__(self, other: int, modulo: int | fmpz | None = None, /) -> fmpz: ...
+    def __pow__(self, other: ifmpz, modulo: int | fmpz | None = None, /) -> fmpz: ...
+    def __rpow__(self, other: int, modulo: ifmpz | None = None, /) -> fmpz: ...
     def __abs__(self) -> fmpz: ...
 
-    def __lshift__(self, other: int | fmpz, /) -> fmpz: ...
-    def __rlshift__(self, other: int | fmpz, /) -> fmpz: ...
-    def __rshift__(self, other: int | fmpz, /) -> fmpz: ...
-    def __rrshift__(self, other: int | fmpz, /) -> fmpz: ...
-    def __and__(self, other: int | fmpz, /) -> fmpz: ...
+    def __lshift__(self, other: ifmpz, /) -> fmpz: ...
+    def __rlshift__(self, other: ifmpz, /) -> fmpz: ...
+    def __rshift__(self, other: ifmpz, /) -> fmpz: ...
+    def __rrshift__(self, other: ifmpz, /) -> fmpz: ...
+    def __and__(self, other: ifmpz, /) -> fmpz: ...
     def __rand__(self, other: int, /) -> fmpz: ...
-    def __or__(self, other: int | fmpz, /) -> fmpz: ...
+    def __or__(self, other: ifmpz, /) -> fmpz: ...
     def __ror__(self, other: int, /) -> fmpz: ...
-    def __xor__(self, other: int | fmpz, /) -> fmpz: ...
+    def __xor__(self, other: ifmpz, /) -> fmpz: ...
     def __rxor__(self, other: int, /) -> fmpz: ...
     def __invert__(self) -> fmpz: ...
 
-    def gcd(self, other: int | fmpz, /) -> fmpz: ...
-    def lcm(self, other: int | fmpz, /) -> fmpz: ...
+    def gcd(self, other: ifmpz, /) -> fmpz: ...
+    def lcm(self, other: ifmpz, /) -> fmpz: ...
     def factor(self, trial_limit: int | None = None) -> list[tuple[fmpz, int]]: ...
-    def factor_smooth(self, bits: int = 15, proved: int = -1): ...
+    def factor_smooth(self, bits: int = 15, proved: int = -1) -> list[tuple[fmpz, int]]: ...
 
     def is_prime(self) -> bool: ...
     def is_probable_prime(self) -> bool: ...
@@ -100,7 +106,7 @@ class fmpz(flint_scalar):
     def isqrt(self) -> fmpz: ...
     def sqrt(self) -> fmpz: ...
     def sqrtrem(self) -> tuple[fmpz, fmpz]: ...
-    def sqrtmod(self, p: int | fmpz) -> fmpz: ...
+    def sqrtmod(self, p: ifmpz) -> fmpz: ...
     def root(self, n: int) -> fmpz: ...
 
-    def jacobi(self, other: int | fmpz, /) -> int: ...
+    def jacobi(self, other: ifmpz, /) -> int: ...
diff --git a/src/flint/types/fmpz_mod.pyi b/src/flint/types/fmpz_mod.pyi
new file mode 100644
index 00000000..23f07e6e
--- /dev/null
+++ b/src/flint/types/fmpz_mod.pyi
@@ -0,0 +1,66 @@
+from ..flint_base.flint_base import flint_scalar
+from .fmpz import fmpz
+from .fmpq import fmpq
+
+
+ifmpz = int | fmpz
+ifmpq = int | fmpz | fmpq
+ifmpz_mod = ifmpz | fmpz_mod
+
+
+class fmpz_mod_ctx:
+
+    def __init__(self, mod: ifmpz) -> None: ...
+    def modulus(self) -> fmpz: ...
+    def is_prime(self) -> bool: ...
+    def zero(self) -> fmpz_mod: ...
+    def one(self) -> fmpz_mod: ...
+    def random_element(self) -> fmpz_mod: ...
+    def __eq__(self, other: object) -> bool: ...
+    def __hash__(self) -> int: ...
+    def __str__(self) -> str: ...
+    def __repr__(self) -> str: ...
+    def __call__(self, val: ifmpz_mod) -> fmpz_mod: ...
+
+
+class fmpz_mod(flint_scalar):
+
+    def __init__(self, val: ifmpz_mod, mod: fmpz_mod_ctx) -> None: ...
+
+    def is_zero(self) -> bool: ...
+
+
+    def is_one(self) -> bool: ...
+    def is_unit(self) -> bool: ...
+
+    def inverse(self, check: bool = True) -> fmpz_mod: ...
+
+    def discrete_log(self, a: ifmpz_mod) -> fmpz: ...
+
+    def __eq__(self, other: object) -> bool: ...
+    def __ne__(self, other: object) -> bool: ...
+
+    def __bool__(self) -> bool: ...
+
+    def repr(self) -> str: ...
+    def __repr__(self) -> str: ...
+    def __str__(self) -> str: ...
+    def __int__(self) -> int: ...
+    def str(self) -> str: ...
+
+    def __pos__(self) -> fmpz_mod: ...
+
+    def __neg__(self) -> fmpz_mod: ...
+
+    def __add__(self, other: ifmpz_mod, /) -> fmpz_mod: ...
+    def __radd__(self, other: ifmpq, /) -> fmpz_mod: ...
+    def __sub__(self, other: ifmpz_mod, /) -> fmpz_mod: ...
+    def __rsub__(self, other: ifmpq, /) -> fmpz_mod: ...
+    def __mul__(self, other: ifmpz_mod, /) -> fmpz_mod: ...
+    def __rmul__(self, other: ifmpq, /) -> fmpz_mod: ...
+    def __truediv__(self, other: ifmpz_mod, /) -> fmpz_mod: ...
+    def __rtruediv__(self, other: ifmpq, /) -> fmpz_mod: ...
+    def __pow__(self, other: ifmpz, /) -> fmpz_mod: ...
+    def __invert__(self) -> fmpz_mod: ...
+
+    def sqrt(self) -> fmpz_mod: ...
diff --git a/src/flint/types/fmpz_mod.pyx b/src/flint/types/fmpz_mod.pyx
index 5559e20e..2700bf65 100644
--- a/src/flint/types/fmpz_mod.pyx
+++ b/src/flint/types/fmpz_mod.pyx
@@ -563,7 +563,7 @@ cdef class fmpz_mod(flint_scalar):
         # Attempt to convert exponent to fmpz
         e = any_as_fmpz(e)
         if e is NotImplemented:
-            raise NotImplementedError
+            raise TypeError("fmpz_mod.__pow__: exponent must be an integer")
 
         check = fmpz_mod_pow_fmpz(
             res.val, self.val, (<fmpz>e).val, self.ctx.val
diff --git a/src/flint/types/fmpz_mod_mpoly.pyi b/src/flint/types/fmpz_mod_mpoly.pyi
new file mode 100644
index 00000000..02db61fd
--- /dev/null
+++ b/src/flint/types/fmpz_mod_mpoly.pyi
@@ -0,0 +1,111 @@
+from typing import Iterable, Mapping
+from ..flint_base.flint_base import flint_mpoly, flint_mod_mpoly_context, Ordering
+from .fmpz import fmpz
+from .fmpz_mod import fmpz_mod
+from .fmpz_mpoly import fmpz_mpoly
+
+
+_str = str
+ifmpz = int | fmpz
+ifmpz_mod = int | fmpz | fmpz_mod
+
+
+class fmpz_mod_mpoly_ctx(flint_mod_mpoly_context[fmpz_mod_mpoly, fmpz_mod, ifmpz_mod]):
+
+    @classmethod
+    def get(cls,
+        names: _str | Iterable[_str | tuple[_str, int]] | tuple[_str, int],
+        ordering: Ordering | _str = Ordering.lex,
+        *,
+        modulus: ifmpz,
+    ) -> fmpz_mod_mpoly_ctx:
+        ...
+
+    def nvars(self) -> int: ...
+    def ordering(self) -> Ordering: ...
+    def modulus(self) -> int: ...
+    def gen(self, i: int, /) -> fmpz_mod_mpoly: ...
+    def constant(self, z: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def from_dict(self, d: Mapping[tuple[int, ...], ifmpz_mod]) -> fmpz_mod_mpoly: ...
+
+
+class fmpz_mod_mpoly(flint_mpoly[fmpz_mod_mpoly_ctx, fmpz_mod, ifmpz_mod]):
+    def __init__(self,
+        val: fmpz_mod_mpoly | fmpz_mpoly | ifmpz_mod | dict[tuple[int, ...], ifmpz_mod] | _str = 0,
+        ctx: fmpz_mod_mpoly_ctx | None = None
+    ) -> None: ...
+
+    def str(self) -> _str: ...
+    def repr(self) -> _str: ...
+
+    def context(self) -> fmpz_mod_mpoly_ctx: ...
+
+    def degrees(self) -> tuple[int, ...]: ...
+    def total_degree(self) -> int: ...
+
+    def is_one(self) -> bool: ...
+    def is_zero(self) -> bool: ...
+    def is_constant(self) -> bool: ...
+
+    def __len__(self) -> int: ...
+    def coefficient(self, i: int) -> fmpz_mod: ...
+    def monomial(self, i: int) -> tuple[int, ...]: ...
+
+    def monoms(self) -> list[tuple[int, ...]]: ...
+    def coeffs(self) -> list[fmpz_mod]: ...
+    def __getitem__(self, index: tuple[int, ...]) -> fmpz_mod: ...
+    def __setitem__(self, index: tuple[int, ...], coeff: ifmpz_mod) -> None: ...
+
+    def subs(self, mapping: dict[_str | int, ifmpz_mod]) -> fmpz_mod_mpoly: ...
+    def compose(self, *args: fmpz_mod_mpoly, ctx: fmpz_mod_mpoly_ctx | None = None) -> fmpz_mod_mpoly: ...
+
+    def __call__(self, *args: ifmpz_mod) -> fmpz_mod: ...
+
+    def __pos__(self) -> fmpz_mod_mpoly: ...
+    def __neg__(self) -> fmpz_mod_mpoly: ...
+    def __add__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __radd__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __sub__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __rsub__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __mul__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __rmul__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __truediv__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __rtruediv__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __floordiv__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __rfloordiv__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __mod__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __rmod__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __divmod__(self, other: fmpz_mod_mpoly | ifmpz_mod) -> tuple[fmpz_mod_mpoly, fmpz_mod_mpoly]: ...
+    def __rdivmod__(self, other: ifmpz_mod) -> tuple[fmpz_mod_mpoly, fmpz_mod_mpoly]: ...
+    def __pow__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+    def __rpow__(self, other: ifmpz_mod) -> fmpz_mod_mpoly: ...
+
+    def iadd(self, other: fmpz_mod_mpoly | ifmpz_mod) -> None: ...
+    def isub(self, other: fmpz_mod_mpoly | ifmpz_mod) -> None: ...
+    def imul(self, other: fmpz_mod_mpoly | ifmpz_mod) -> None: ...
+
+    def gcd(self, other: fmpz_mod_mpoly) -> fmpz_mod_mpoly: ...
+    def xgcd(self, other: fmpz_mod_mpoly, /) -> fmpz_mod: ...
+    def term_content(self) -> fmpz_mod_mpoly: ...
+    def primitive(self) -> tuple[fmpz_mod, fmpz_mod_mpoly]: ...
+
+    def factor(self) -> tuple[fmpz_mod, list[tuple[fmpz_mod_mpoly, int]]]: ...
+    def factor_squarefree(self) -> tuple[fmpz_mod, list[tuple[fmpz_mod_mpoly, int]]]: ...
+
+    def sqrt(self, assume_perfect_square: bool = False) -> fmpz_mod_mpoly: ...
+
+    def resultant(self, other: fmpz_mod_mpoly, var: _str | int) -> fmpz_mod_mpoly: ...
+    def discriminant(self, var: _str | int) -> fmpz_mod_mpoly: ...
+
+    def deflation(self) -> tuple[fmpz_mod_mpoly, list[int]]: ...
+    def inflate(self, N: list[int]) -> fmpz_mod_mpoly: ...
+    def deflate(self, N: list[int]) -> fmpz_mod_mpoly: ...
+    def deflation_monom(self) -> tuple[fmpz_mod_mpoly, list[int], fmpz_mod_mpoly]: ...
+    def deflation_index(self) -> tuple[list[int], list[int]]: ...
+
+    def derivative(self, var: _str | int) -> fmpz_mod_mpoly: ...
+
+
+class fmpz_mod_mpoly_vec:
+    ...
+
diff --git a/src/flint/types/fmpz_mpoly.pyi b/src/flint/types/fmpz_mpoly.pyi
index fe24f210..1ec15c31 100644
--- a/src/flint/types/fmpz_mpoly.pyi
+++ b/src/flint/types/fmpz_mpoly.pyi
@@ -4,14 +4,15 @@ from ..flint_base.flint_base import flint_mpoly, flint_mpoly_context, Ordering
 from .fmpz import fmpz
 
 
+_str = str
 ifmpz = int | fmpz
 
 
 class fmpz_mpoly_ctx(flint_mpoly_context[fmpz_mpoly, fmpz, ifmpz]):
     @classmethod
     def get(cls,
-        names: str | Iterable[str | tuple[str, int]] | tuple[str, int],
-        ordering: Ordering | str = Ordering.lex
+        names: _str | Iterable[_str | tuple[_str, int]] | tuple[_str, int],
+        ordering: Ordering | _str = Ordering.lex
     ) -> fmpz_mpoly_ctx:
         ...
 
@@ -26,12 +27,12 @@ class fmpz_mpoly_ctx(flint_mpoly_context[fmpz_mpoly, fmpz, ifmpz]):
 class fmpz_mpoly(flint_mpoly[fmpz_mpoly_ctx, fmpz, ifmpz]):
 
     def __init__(self,
-        val: fmpz_mpoly | ifmpz | dict[tuple[int, ...], ifmpz] | str = 0,
+        val: fmpz_mpoly | ifmpz | dict[tuple[int, ...], ifmpz] | _str = 0,
         ctx: fmpz_mpoly_ctx | None = None
     ) -> None: ...
 
-    def str(self) -> str: ...
-    def repr(self) -> str: ...
+    def str(self) -> _str: ...
+    def repr(self) -> _str: ...
 
     def context(self) -> fmpz_mpoly_ctx: ...
 
@@ -51,7 +52,7 @@ class fmpz_mpoly(flint_mpoly[fmpz_mpoly_ctx, fmpz, ifmpz]):
     def __getitem__(self, index: tuple[int, ...]) -> fmpz: ...
     def __setitem__(self, index: tuple[int, ...], coeff: ifmpz) -> None: ...
 
-    def subs(self, mapping: dict[str | int, ifmpz]) -> fmpz_mpoly: ...
+    def subs(self, mapping: dict[_str | int, ifmpz]) -> fmpz_mpoly: ...
     def compose(self, *args: fmpz_mpoly, ctx: fmpz_mpoly_ctx | None = None) -> fmpz_mpoly: ...
 
     def __call__(self, *args: ifmpz) -> fmpz: ...
@@ -88,8 +89,8 @@ class fmpz_mpoly(flint_mpoly[fmpz_mpoly_ctx, fmpz, ifmpz]):
 
     def sqrt(self, assume_perfect_square: bool = False) -> fmpz_mpoly: ...
 
-    def resultant(self, other: fmpz_mpoly, var: str | int) -> fmpz_mpoly: ...
-    def discriminant(self, var: str | int) -> fmpz_mpoly: ...
+    def resultant(self, other: fmpz_mpoly, var: _str | int) -> fmpz_mpoly: ...
+    def discriminant(self, var: _str | int) -> fmpz_mpoly: ...
 
     def deflation(self) -> tuple[fmpz_mpoly, list[int]]: ...
     def inflate(self, N: list[int]) -> fmpz_mpoly: ...
@@ -97,7 +98,8 @@ class fmpz_mpoly(flint_mpoly[fmpz_mpoly_ctx, fmpz, ifmpz]):
     def deflation_monom(self) -> tuple[fmpz_mpoly, list[int], fmpz_mpoly]: ...
     def deflation_index(self) -> tuple[list[int], list[int]]: ...
 
-    def derivative(self, var: str | int) -> fmpz_mpoly: ...
+    def derivative(self, var: _str | int) -> fmpz_mpoly: ...
+    def integral(self, var: _str | int) -> tuple[int, fmpz_mpoly]: ...
 
 
 class fmpz_mpoly_vec:
diff --git a/src/flint/types/fmpz_poly.pyi b/src/flint/types/fmpz_poly.pyi
new file mode 100644
index 00000000..32b7d296
--- /dev/null
+++ b/src/flint/types/fmpz_poly.pyi
@@ -0,0 +1,109 @@
+from typing import overload, Any
+from flint.flint_base.flint_base import flint_poly
+from flint.types.fmpz import fmpz, ifmpz
+from flint.types.fmpq import fmpq
+from flint.types.fmpq_poly import fmpq_poly
+
+
+ifmpz_poly = fmpz_poly | ifmpz
+
+
+class fmpz_poly(flint_poly[fmpz]):
+    """
+    The *fmpz_poly* type represents dense univariate polynomials over
+    the integers.
+    """
+
+    @overload
+    def __init__(self) -> None: ...
+    @overload
+    def __init__(self, arg: list[int] | list[ifmpz] | ifmpz | fmpz_poly, /) -> None: ...
+
+    def __len__(self) -> int: ...
+    def length(self) -> int: ...
+    def degree(self) -> int: ...
+
+    def __getitem__(self, i: int, /) -> fmpz: ...
+    def __setitem__(self, i: int, x: ifmpz, /) -> None: ...
+
+    def repr(self) -> str: ...
+
+    def __bool__(self) -> bool: ...
+    def is_zero(self) -> bool: ...
+    def is_one(self) -> bool: ...
+    def is_constant(self) -> bool: ...
+
+    def leading_coefficient(self) -> fmpz: ...
+
+    @overload
+    def __call__(self, other: ifmpz, /) -> fmpz: ...
+    @overload
+    def __call__(self, other: fmpq, /) -> fmpq: ...
+    @overload
+    def __call__(self, other: fmpz_poly, /) -> fmpz_poly: ...
+    @overload
+    def __call__(self, other: fmpq_poly, /) -> fmpq_poly: ...
+    #@overload
+    #def __call__(self, other: arb, /) -> acb: ...
+    #@overload
+    #def __call__(self, other: acb, /) -> acb: ...
+
+
+    def derivative(self) -> fmpz_poly: ...
+
+    def __pos__(self) -> fmpz_poly: ...
+    def __neg__(self) -> fmpz_poly: ...
+    def __add__(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def __radd__(self, other: ifmpz, /) -> fmpz_poly: ...
+    def __sub__(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def __rsub__(self, other: ifmpz, /) -> fmpz_poly: ...
+    def __mul__(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def __rmul__(self, other: ifmpz, /) -> fmpz_poly: ...
+    def __truediv__(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def __rtruediv__(self, other: ifmpz, /) -> fmpz_poly: ...
+    def __floordiv__(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def __rfloordiv__(self, other: ifmpz, /) -> fmpz_poly: ...
+    def __mod__(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def __rmod__(self, other: ifmpz, /) -> fmpz_poly: ...
+    def __divmod__(self, other: ifmpz_poly, /) -> tuple[fmpz_poly, fmpz_poly]: ...
+    def __rdivmod__(self, other: ifmpz, /) -> tuple[fmpz_poly, fmpz_poly]: ...
+    def __pow__(self, other: int, /) -> fmpz_poly: ...
+
+    def gcd(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def content(self) -> fmpz: ...
+
+    def resultant(self, other: ifmpz_poly, /) -> fmpz_poly: ...
+    def factor(self) -> tuple[fmpz_poly, list[tuple[fmpz_poly, int]]]: ...
+    def factor_squarefree(self) -> tuple[fmpz_poly, list[tuple[fmpz_poly, int]]]: ...
+    def sqrt(self) -> fmpz_poly: ...
+
+    def deflation(self) -> tuple[fmpz_poly, int]: ...
+    def deflation_monom(self) -> tuple[fmpz_poly, int, fmpz_poly]: ...
+    def inflate(self, n: int, /) -> fmpz_poly: ...
+    def deflate(self, n: int, /) -> fmpz_poly: ...
+    def deflation_index(self) -> tuple[int, int]: ...
+
+    # XXX acb
+    def complex_roots(self, verbose: bool = False) -> list[Any]: ...
+
+    @staticmethod
+    def cyclotomic(n: int, /) -> fmpz_poly: ...
+
+    def is_cyclotomic(self) -> int: ...
+
+    @staticmethod
+    def cos_minpoly(n: int, /) -> fmpz_poly: ...
+
+    @staticmethod
+    def chebyshev_t(n: int, /) -> fmpz_poly: ...
+
+    @staticmethod
+    def chebyshev_u(n: int, /) -> fmpz_poly: ...
+
+    @staticmethod
+    def swinnerton_dyer(n: int, /, use_arb: bool = True) -> fmpz_poly: ...
+
+    @staticmethod
+    def hilbert_class_poly(D: int, /) -> fmpz_poly: ...
+
+    def height_bits(self, signed: bool = False) -> int: ...
diff --git a/src/flint/types/fmpz_series.pyi b/src/flint/types/fmpz_series.pyi
new file mode 100644
index 00000000..a34d5bd5
--- /dev/null
+++ b/src/flint/types/fmpz_series.pyi
@@ -0,0 +1,46 @@
+from typing import Any, Sequence
+from flint.flint_base.flint_base import flint_series
+from flint.types.fmpz import fmpz, ifmpz
+from flint.types.fmpz_poly import fmpz_poly, ifmpz_poly
+
+
+ifmpz_series = fmpz_series | ifmpz_poly
+
+
+class fmpz_series(flint_series[fmpz]):
+
+    def __init__(self,
+                 val: list[int] | list[ifmpz] | fmpz_series | fmpz_poly | ifmpz | None = None,
+                 prec: int | None = None): ...
+
+    @property
+    def prec(self) -> int: ...
+
+    def _equal_repr(self, other: object) -> bool: ...
+
+    def __len__(self) -> int: ...
+    def length(self) -> int: ...
+
+    def __getitem__(self, i: int, /) -> fmpz: ...
+    def __setitem__(self, i: int, x: ifmpz, /) -> None: ...
+
+    def repr(self, **kwargs: Any) -> str: ...
+    def str(self, **kwargs: Any) -> str: ...
+
+    def __pos__(self) -> fmpz_series: ...
+    def __neg__(self) -> fmpz_series: ...
+    def __add__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __radd__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __sub__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __rsub__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __mul__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __rmul__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __truediv__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __rtruediv__(self, other: ifmpz_series, /) -> fmpz_series: ...
+    def __pow__(self, exp: int, /) -> fmpz_series: ...
+
+    def valuation(self) -> int: ...
+
+    def __call__(self, other: ifmpz_series, /) -> fmpz_series: ...
+
+    def reversion(self) -> fmpz_series: ...
diff --git a/src/flint/types/nmod.pyi b/src/flint/types/nmod.pyi
index 449f89d4..8dbe8a4e 100644
--- a/src/flint/types/nmod.pyi
+++ b/src/flint/types/nmod.pyi
@@ -2,6 +2,7 @@ from ..flint_base.flint_base import flint_scalar
 from .fmpz import fmpz
 from .fmpq import fmpq
 
+_str = str
 ifmpz = int | fmpz
 ifmpq = int | fmpz | fmpq
 inmod = int | fmpz | fmpq | nmod
@@ -11,11 +12,11 @@ class nmod(flint_scalar):
 
     def modulus(self) -> int: ...
 
-    def str(self, base: int = 10, condense: int = 0) -> str: ...
-    def repr(self) -> str: ...
+    def str(self, base: int = 10, condense: int = 0) -> _str: ...
+    def repr(self) -> _str: ...
 
-    def __str__(self) -> str: ...
-    def __repr__(self) -> str: ...
+    def __str__(self) -> _str: ...
+    def __repr__(self) -> _str: ...
 
     def __int__(self) -> int: ...
 
diff --git a/src/flint/types/nmod_mpoly.pyi b/src/flint/types/nmod_mpoly.pyi
index 4dca7750..b2eed49b 100644
--- a/src/flint/types/nmod_mpoly.pyi
+++ b/src/flint/types/nmod_mpoly.pyi
@@ -1,4 +1,4 @@
-from typing import Iterable, Mapping
+from typing import Iterable, Mapping, Self
 from ..flint_base.flint_base import flint_mpoly, flint_mpoly_context, Ordering
 from .fmpz import fmpz
 from .fmpq import fmpq
@@ -6,6 +6,7 @@ from .nmod import nmod
 from .fmpz_mpoly import fmpz_mpoly
 
 
+_str = str
 ifmpz = int | fmpz
 inmod = int | fmpz | fmpq | nmod
 
@@ -14,8 +15,8 @@ class nmod_mpoly_ctx(flint_mpoly_context[nmod_mpoly, nmod, inmod]):
 
     @classmethod
     def get(cls,
-        names: str | Iterable[str | tuple[str, int]] | tuple[str, int],
-        ordering: Ordering | str = Ordering.lex,
+        names: _str | Iterable[_str | tuple[_str, int]] | tuple[_str, int],
+        ordering: Ordering | _str = Ordering.lex,
         *,
         modulus: int,
     ) -> nmod_mpoly_ctx:
@@ -31,12 +32,12 @@ class nmod_mpoly_ctx(flint_mpoly_context[nmod_mpoly, nmod, inmod]):
 
 class nmod_mpoly(flint_mpoly[nmod_mpoly_ctx, nmod, inmod]):
     def __init__(self,
-        val: nmod_mpoly | fmpz_mpoly | inmod | dict[tuple[int, ...], inmod] | str = 0,
+        val: nmod_mpoly | fmpz_mpoly | inmod | dict[tuple[int, ...], inmod] | _str = 0,
         ctx: nmod_mpoly_ctx | None = None
     ) -> None: ...
 
-    def str(self) -> str: ...
-    def repr(self) -> str: ...
+    def str(self) -> _str: ...
+    def repr(self) -> _str: ...
 
     def context(self) -> nmod_mpoly_ctx: ...
 
@@ -56,7 +57,7 @@ class nmod_mpoly(flint_mpoly[nmod_mpoly_ctx, nmod, inmod]):
     def __getitem__(self, index: tuple[int, ...]) -> nmod: ...
     def __setitem__(self, index: tuple[int, ...], coeff: inmod) -> None: ...
 
-    def subs(self, mapping: dict[str | int, inmod]) -> nmod_mpoly: ...
+    def subs(self, mapping: dict[_str | int, inmod]) -> nmod_mpoly: ...
     def compose(self, *args: nmod_mpoly, ctx: nmod_mpoly_ctx | None = None) -> nmod_mpoly: ...
 
     def __call__(self, *args: inmod) -> nmod: ...
@@ -94,8 +95,8 @@ class nmod_mpoly(flint_mpoly[nmod_mpoly_ctx, nmod, inmod]):
 
     def sqrt(self, assume_perfect_square: bool = False) -> nmod_mpoly: ...
 
-    def resultant(self, other: nmod_mpoly, var: str | int) -> nmod_mpoly: ...
-    def discriminant(self, var: str | int) -> nmod_mpoly: ...
+    def resultant(self, other: Self, var: _str | int) -> Self: ...
+    def discriminant(self, var: _str | int) -> nmod_mpoly: ...
 
     def deflation(self) -> tuple[nmod_mpoly, list[int]]: ...
     def inflate(self, N: list[int]) -> nmod_mpoly: ...
@@ -103,9 +104,9 @@ class nmod_mpoly(flint_mpoly[nmod_mpoly_ctx, nmod, inmod]):
     def deflation_monom(self) -> tuple[nmod_mpoly, list[int], nmod_mpoly]: ...
     def deflation_index(self) -> tuple[list[int], list[int]]: ...
 
-    def derivative(self, var: str | int) -> nmod_mpoly: ...
+    def derivative(self, var: _str | int) -> nmod_mpoly: ...
 
 
-class fmpq_mpoly_vec:
+class nmod_mpoly_vec:
     ...