diff --git a/setup.py b/setup.py
index 2b4a1ace..2cb837bd 100644
--- a/setup.py
+++ b/setup.py
@@ -97,9 +97,11 @@
     ("flint.types.nmod_poly", ["src/flint/types/nmod_poly.pyx"]),
     ("flint.types.nmod_mat", ["src/flint/types/nmod_mat.pyx"]),
     ("flint.types.nmod_series", ["src/flint/types/nmod_series.pyx"]),
+    ("flint.types.nmod_mpoly", ["src/flint/types/nmod_mpoly.pyx"]),
 
     ("flint.types.fmpz_mod", ["src/flint/types/fmpz_mod.pyx"]),
     ("flint.types.fmpz_mod_poly", ["src/flint/types/fmpz_mod_poly.pyx"]),
+    ("flint.types.fmpz_mod_mpoly", ["src/flint/types/fmpz_mod_mpoly.pyx"]),
     ("flint.types.fmpz_mod_mat", ["src/flint/types/fmpz_mod_mat.pyx"]),
 
     ("flint.types.fmpq_mpoly", ["src/flint/types/fmpq_mpoly.pyx"]),
diff --git a/src/flint/__init__.py b/src/flint/__init__.py
index 4786804b..1b26c1ff 100644
--- a/src/flint/__init__.py
+++ b/src/flint/__init__.py
@@ -14,12 +14,14 @@
 
 from .types.nmod import *
 from .types.nmod_poly import *
+from .types.nmod_mpoly import nmod_mpoly_ctx, nmod_mpoly, nmod_mpoly_vec
 from .types.nmod_mat import *
 from .types.nmod_series import *
 
 from .types.fmpz_mpoly import fmpz_mpoly_ctx, fmpz_mpoly, fmpz_mpoly_vec
 from .types.fmpz_mod import *
 from .types.fmpz_mod_poly import *
+from .types.fmpz_mod_mpoly import fmpz_mod_mpoly_ctx, fmpz_mod_mpoly, fmpz_mod_mpoly_vec
 from .types.fmpz_mod_mat import fmpz_mod_mat
 
 from .types.fmpq_mpoly import fmpq_mpoly_ctx, fmpq_mpoly, fmpq_mpoly_vec
diff --git a/src/flint/flint_base/flint_base.pyx b/src/flint/flint_base/flint_base.pyx
index e6151845..e0b02fcf 100644
--- a/src/flint/flint_base/flint_base.pyx
+++ b/src/flint/flint_base/flint_base.pyx
@@ -11,6 +11,9 @@ from flint.utils.typecheck cimport typecheck
 cimport libc.stdlib
 
 from typing import Optional
+from flint.utils.flint_exceptions import IncompatibleContextError
+
+from flint.types.fmpz cimport fmpz, any_as_fmpz
 
 
 FLINT_BITS = _FLINT_BITS
@@ -38,7 +41,7 @@ cdef class flint_scalar(flint_elem):
     def is_zero(self):
         return False
 
-    def _any_as_self(self):
+    def _any_as_self(self, other):
         return NotImplemented
 
     def _neg_(self):
@@ -329,14 +332,13 @@ cdef class flint_mpoly_context(flint_elem):
         return nametup
 
     @classmethod
-    def get_context(cls, slong nvars=1, ordering=Ordering.lex, names: Optional[str] = "x", nametup: Optional[tuple] = None):
+    def create_context_key(cls, slong nvars=1, ordering=Ordering.lex, names: Optional[str] = "x", nametup: Optional[tuple] = None):
         """
-        Retrieve a context via the number of variables, `nvars`, the ordering, `ordering`, and either a variable
-        name string, `names`, or a tuple of variable names, `nametup`.
+        Create a key for the context cache via the number of variables, the ordering, and
+        either a variable name string, or a tuple of variable names.
         """
-
         # A type hint of `ordering: Ordering` results in the error "TypeError: an integer is required" if a Ordering
-        # object is not provided. This is pretty obtuse so we check it's type ourselves
+        # object is not provided. This is pretty obtuse so we check its type ourselves
         if not isinstance(ordering, Ordering):
             raise TypeError(f"`ordering` ('{ordering}') is not an instance of flint.Ordering")
 
@@ -346,6 +348,15 @@ cdef class flint_mpoly_context(flint_elem):
             key = nvars, ordering, cls.create_variable_names(nvars, names)
         else:
             raise ValueError("must provide either `names` or `nametup`")
+        return key
+
+    @classmethod
+    def get_context(cls, *args, **kwargs):
+        """
+        Retrieve a context via the number of variables, `nvars`, the ordering, `ordering`, and either a variable
+        name string, `names`, or a tuple of variable names, `nametup`.
+        """
+        key = cls.create_context_key(*args, **kwargs)
 
         ctx = cls._ctx_cache.get(key)
         if ctx is None:
@@ -361,6 +372,18 @@ cdef class flint_mpoly_context(flint_elem):
             nametup=ctx.names()
         )
 
+    def any_as_scalar(self, other):
+        raise NotImplementedError("abstract method")
+
+    def scalar_as_mpoly(self, other):
+        raise NotImplementedError("abstract method")
+
+    def compatible_context_check(self, other):
+        if not typecheck(other, type(self)):
+            raise TypeError(f"type {type(other)} is not {type(self)}")
+        elif other is not self:
+            raise IncompatibleContextError(f"{other} is not {self}")
+
 
 cdef class flint_mpoly(flint_elem):
     """
@@ -373,6 +396,281 @@ cdef class flint_mpoly(flint_elem):
     def to_dict(self):
         return {self.monomial(i): self.coefficient(i) for i in range(len(self))}
 
+    def _division_check(self, other):
+        if not other:
+            raise ZeroDivisionError("nmod_mpoly division by zero")
+
+    def _add_scalar_(self, other):
+        return NotImplemented
+
+    def _add_mpoly_(self, other):
+        return NotImplemented
+
+    def _iadd_scalar_(self, other):
+        return NotImplemented
+
+    def _iadd_mpoly_(self, other):
+        return NotImplemented
+
+    def _sub_scalar_(self, other):
+        return NotImplemented
+
+    def _sub_mpoly_(self, other):
+        return NotImplemented
+
+    def _isub_scalar_(self, other):
+        return NotImplemented
+
+    def _isub_mpoly_(self, other):
+        return NotImplemented
+
+    def _mul_scalar_(self, other):
+        return NotImplemented
+
+    def _imul_mpoly_(self, other):
+        return NotImplemented
+
+    def _imul_scalar_(self, other):
+        return NotImplemented
+
+    def _mul_mpoly_(self, other):
+        return NotImplemented
+
+    def _pow_(self, other):
+        return NotImplemented
+
+    def _divmod_mpoly_(self, other):
+        return NotImplemented
+
+    def _floordiv_mpoly_(self, other):
+        return NotImplemented
+
+    def _truediv_mpoly_(self, other):
+        return NotImplemented
+
+    def __add__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            return self._add_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        return self._add_scalar_(other)
+
+    def __radd__(self, other):
+        return self.__add__(other)
+
+    def iadd(self, other):
+        """
+        In-place addition, mutates self.
+
+            >>> from flint import Ordering, fmpz_mpoly_ctx
+            >>> ctx = fmpz_mpoly_ctx.get_context(2, Ordering.lex, 'x')
+            >>> f = ctx.from_dict({(1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f
+            4*x0*x1 + 2*x0 + 3*x1
+            >>> f.iadd(5)
+            >>> f
+            4*x0*x1 + 2*x0 + 3*x1 + 5
+
+        """
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._iadd_mpoly_(other)
+            return
+
+        other_scalar = self.context().any_as_scalar(other)
+        if other_scalar is NotImplemented:
+            raise NotImplementedError(f"cannot add {type(self)} and {type(other)}")
+
+        self._iadd_scalar_(other_scalar)
+
+    def __sub__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            return self._sub_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        return self._sub_scalar_(other)
+
+    def __rsub__(self, other):
+        return -self.__sub__(other)
+
+    def isub(self, other):
+        """
+        In-place subtraction, mutates self.
+
+            >>> from flint import Ordering, fmpz_mpoly_ctx
+            >>> ctx = fmpz_mpoly_ctx.get_context(2, Ordering.lex, 'x')
+            >>> f = ctx.from_dict({(1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f
+            4*x0*x1 + 2*x0 + 3*x1
+            >>> f.isub(5)
+            >>> f
+            4*x0*x1 + 2*x0 + 3*x1 - 5
+
+        """
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._isub_mpoly_(other)
+            return
+
+        other_scalar = self.context().any_as_scalar(other)
+        if other_scalar is NotImplemented:
+            raise NotImplementedError(f"cannot subtract {type(self)} and {type(other)}")
+
+        self._isub_scalar_(other_scalar)
+
+    def __mul__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            return self._mul_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        return self._mul_scalar_(other)
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
+    def imul(self, other):
+        """
+        In-place multiplication, mutates self.
+
+            >>> from flint import Ordering, fmpz_mpoly_ctx
+            >>> ctx = fmpz_mpoly_ctx.get_context(2, Ordering.lex, 'x')
+            >>> f = ctx.from_dict({(1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f
+            4*x0*x1 + 2*x0 + 3*x1
+            >>> f.imul(2)
+            >>> f
+            8*x0*x1 + 4*x0 + 6*x1
+
+        """
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._imul_mpoly_(other)
+            return
+
+        other_scalar = self.context().any_as_scalar(other)
+        if other_scalar is NotImplemented:
+            raise NotImplementedError(f"cannot multiply {type(self)} and {type(other)}")
+
+        self._imul_scalar_(other_scalar)
+
+    def __pow__(self, other, modulus):
+        if modulus is not None:
+            raise NotImplementedError("cannot specify modulus outside of the context")
+        elif typecheck(other, fmpz):
+            return self._pow_(other)
+
+        other = any_as_fmpz(other)
+        if other is NotImplemented:
+            return NotImplemented
+        elif other < 0:
+            raise ValueError("cannot raise to a negative power")
+
+        return self._pow_(other)
+
+    def __divmod__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._division_check(other)
+            return self._divmod_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        self._division_check(other)
+        return self._divmod_mpoly_(other)
+
+    def __rdivmod__(self, other):
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        other._division_check(self)
+        return other._divmod_mpoly_(self)
+
+    def __truediv__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._division_check(other)
+            return self._truediv_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        self._division_check(other)
+        return self._truediv_mpoly_(other)
+
+    def __rtruediv__(self, other):
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        other._division_check(self)
+        return other._truediv_mpoly_(self)
+
+    def __floordiv__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._division_check(other)
+            return self._floordiv_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        self._division_check(other)
+        return self._floordiv_mpoly_(other)
+
+    def __rfloordiv__(self, other):
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        other._division_check(self)
+        return other._floordiv_mpoly_(self)
+
+    def __mod__(self, other):
+        if typecheck(other, type(self)):
+            self.context().compatible_context_check(other.context())
+            self._division_check(other)
+            return self._mod_mpoly_(other)
+
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        self._division_check(other)
+        return self._mod_mpoly_(other)
+
+    def __rmod__(self, other):
+        other = self.context().any_as_scalar(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        other = self.context().scalar_as_mpoly(other)
+        other._division_check(self)
+        return other._mod_mpoly_(self)
+
     def __contains__(self, x):
         """
         Returns True if `self` contains a term with exponent vector `x` and a non-zero coefficient.
diff --git a/src/flint/flintlib/flint.pxd b/src/flint/flintlib/flint.pxd
index 65ed5d84..60be1da3 100644
--- a/src/flint/flintlib/flint.pxd
+++ b/src/flint/flintlib/flint.pxd
@@ -50,6 +50,10 @@ cdef extern from *:
     """
 
 cdef extern from "flint/flint.h":
+    """
+    #define SIZEOF_ULONG sizeof(ulong)
+    """
+    int SIZEOF_ULONG
     const char * FLINT_VERSION
     const int __FLINT_RELEASE
     const int FLINT_BITS
@@ -60,6 +64,11 @@ cdef extern from "flint/flint.h":
     long flint_get_num_threads()
     void flint_cleanup()
 
+    ctypedef struct nmod_t:
+        mp_limb_t n
+        mp_limb_t ninv
+        flint_bitcnt_t norm
+
 cdef extern from *:
     """
     /* FLINT_BITS is not known until C compile time. We need to check if long
diff --git a/src/flint/flintlib/fmpz_mod_mpoly.pxd b/src/flint/flintlib/fmpz_mod_mpoly.pxd
new file mode 100644
index 00000000..bb00391b
--- /dev/null
+++ b/src/flint/flintlib/fmpz_mod_mpoly.pxd
@@ -0,0 +1,176 @@
+from flint.flintlib.fmpz cimport fmpz_t
+from flint.flintlib.flint cimport flint_bitcnt_t, ulong, mp_limb_t, slong, fmpz_struct, flint_rand_t
+from flint.flintlib.fmpz_poly cimport fmpz_poly_struct, fmpz_poly_t
+from flint.flintlib.mpoly cimport mpoly_ctx_t, ordering_t
+from flint.flintlib.fmpz_mod cimport fmpz_mod_ctx_t
+
+cdef extern from "flint/fmpz_mod_mpoly.h":
+    ctypedef struct fmpz_mod_mpoly_ctx_struct:
+        mpoly_ctx_t minfo
+        fmpz_mod_ctx_t ffinfo
+
+    ctypedef fmpz_mod_mpoly_ctx_struct fmpz_mod_mpoly_ctx_t[1]
+
+    ctypedef struct fmpz_mod_mpoly_struct:
+        fmpz_struct * coeffs
+        ulong * exps
+        slong length
+        flint_bitcnt_t bits
+        slong coeffs_alloc
+        slong exps_alloc
+
+    ctypedef fmpz_mod_mpoly_struct fmpz_mod_mpoly_t[1]
+
+    ctypedef struct fmpz_mod_mpoly_univar_struct:
+        fmpz_mod_mpoly_struct * coeffs
+        fmpz_struct * exps
+        slong alloc
+        slong length
+
+    ctypedef fmpz_mod_mpoly_univar_struct fmpz_mod_mpoly_univar_t[1]
+
+# from here on is parsed
+    void fmpz_mod_mpoly_ctx_init(fmpz_mod_mpoly_ctx_t ctx, slong nvars, const ordering_t ord, const fmpz_t p)
+    slong fmpz_mod_mpoly_ctx_nvars(const fmpz_mod_mpoly_ctx_t ctx)
+    ordering_t fmpz_mod_mpoly_ctx_ord(const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_ctx_get_modulus(fmpz_t n, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_ctx_clear(fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_init(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_init2(fmpz_mod_mpoly_t A, slong alloc, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_init3(fmpz_mod_mpoly_t A, slong alloc, flint_bitcnt_t bits, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_clear(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    char * fmpz_mod_mpoly_get_str_pretty(const fmpz_mod_mpoly_t A, const char ** x, const fmpz_mod_mpoly_ctx_t ctx)
+    # int fmpz_mod_mpoly_fprint_pretty(FILE * file, const fmpz_mod_mpoly_t A, const char ** x, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_print_pretty(const fmpz_mod_mpoly_t A, const char ** x, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_set_str_pretty(fmpz_mod_mpoly_t A, const char * str, const char ** x, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_gen(fmpz_mod_mpoly_t A, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_is_gen(const fmpz_mod_mpoly_t A, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_equal(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_swap(fmpz_mod_mpoly_t poly1, fmpz_mod_mpoly_t poly2, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_is_fmpz(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_fmpz(fmpz_t c, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_fmpz(fmpz_mod_mpoly_t A, const fmpz_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_ui(fmpz_mod_mpoly_t A, ulong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_si(fmpz_mod_mpoly_t A, slong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_zero(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_one(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_equal_fmpz(const fmpz_mod_mpoly_t A, const fmpz_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_equal_ui(const fmpz_mod_mpoly_t A, ulong c, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_equal_si(const fmpz_mod_mpoly_t A, slong c, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_is_zero(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_is_one(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_degrees_fit_si(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_degrees_fmpz(fmpz_struct ** degs, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_degrees_si(slong * degs, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_degree_fmpz(fmpz_t deg, const fmpz_mod_mpoly_t A, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_degree_si(const fmpz_mod_mpoly_t A, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_total_degree_fits_si(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_total_degree_fmpz(fmpz_t tdeg, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_total_degree_si(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_used_vars(int * used, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_coeff_fmpz_monomial(fmpz_t c, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t M, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_fmpz_monomial(fmpz_mod_mpoly_t A, const fmpz_t c, const fmpz_mod_mpoly_t M, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_coeff_fmpz_fmpz(fmpz_t c, const fmpz_mod_mpoly_t A, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_coeff_fmpz_ui(fmpz_t c, const fmpz_mod_mpoly_t A, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_fmpz_fmpz(fmpz_mod_mpoly_t A, const fmpz_t c, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_ui_fmpz(fmpz_mod_mpoly_t A, ulong c, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_si_fmpz(fmpz_mod_mpoly_t A, slong c, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_fmpz_ui(fmpz_mod_mpoly_t A, const fmpz_t c, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_ui_ui(fmpz_mod_mpoly_t A, ulong c, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_coeff_si_ui(fmpz_mod_mpoly_t A, slong c, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_coeff_vars_ui(fmpz_mod_mpoly_t C, const fmpz_mod_mpoly_t A, const slong * vars, const ulong * exps, slong length, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_cmp(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_is_canonical(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_length(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_resize(fmpz_mod_mpoly_t A, slong new_length, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_term_coeff_fmpz(fmpz_t c, const fmpz_mod_mpoly_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_term_coeff_fmpz(fmpz_mod_mpoly_t A, slong i, const fmpz_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_term_coeff_ui(fmpz_mod_mpoly_t A, slong i, ulong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_term_coeff_si(fmpz_mod_mpoly_t A, slong i, slong c, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_term_exp_fits_si(const fmpz_mod_mpoly_t poly, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_term_exp_fits_ui(const fmpz_mod_mpoly_t poly, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_term_exp_fmpz(fmpz_struct ** exp, const fmpz_mod_mpoly_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_term_exp_ui(ulong * exp, const fmpz_mod_mpoly_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_term_exp_si(slong * exp, const fmpz_mod_mpoly_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    ulong fmpz_mod_mpoly_get_term_var_exp_ui(const fmpz_mod_mpoly_t A, slong i, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_get_term_var_exp_si(const fmpz_mod_mpoly_t A, slong i, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_term_exp_fmpz(fmpz_mod_mpoly_t A, slong i, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_set_term_exp_ui(fmpz_mod_mpoly_t A, slong i, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_term(fmpz_mod_mpoly_t M, const fmpz_mod_mpoly_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_get_term_monomial(fmpz_mod_mpoly_t M, const fmpz_mod_mpoly_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_fmpz_fmpz(fmpz_mod_mpoly_t A, const fmpz_t c, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_fmpz_ffmpz(fmpz_mod_mpoly_t A, const fmpz_t c, const fmpz_struct * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_ui_fmpz(fmpz_mod_mpoly_t A, ulong c, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_ui_ffmpz(fmpz_mod_mpoly_t A, ulong c, const fmpz_struct * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_si_fmpz(fmpz_mod_mpoly_t A, slong c, fmpz_struct * const * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_si_ffmpz(fmpz_mod_mpoly_t A, slong c, const fmpz_struct * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_fmpz_ui(fmpz_mod_mpoly_t A, const fmpz_t c, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_ui_ui(fmpz_mod_mpoly_t A, ulong c, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_push_term_si_ui(fmpz_mod_mpoly_t A, slong c, const ulong * exp, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_sort_terms(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_combine_like_terms(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_reverse(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_randtest_bound(fmpz_mod_mpoly_t A, flint_rand_t state, slong length, ulong exp_bound, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_randtest_bounds(fmpz_mod_mpoly_t A, flint_rand_t state, slong length, ulong * exp_bounds, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_randtest_bits(fmpz_mod_mpoly_t A, flint_rand_t state, slong length, mp_limb_t exp_bits, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_add_fmpz(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_add_ui(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, ulong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_add_si(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, slong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_sub_fmpz(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_sub_ui(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, ulong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_sub_si(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, slong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_add(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_t C, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_sub(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_t C, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_neg(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_scalar_mul_fmpz(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_scalar_mul_ui(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, ulong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_scalar_mul_si(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, slong c, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_scalar_addmul_fmpz(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_t C, const fmpz_t d, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_make_monic(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_derivative(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_evaluate_all_fmpz(fmpz_t eval, const fmpz_mod_mpoly_t A, fmpz_struct * const * vals, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_evaluate_one_fmpz(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, slong var, const fmpz_t val, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_compose_fmpz_poly(fmpz_poly_t A, const fmpz_mod_mpoly_t B, fmpz_poly_struct * const * C, const fmpz_mod_mpoly_ctx_t ctxB)
+    int fmpz_mod_mpoly_compose_fmpz_mod_mpoly_geobucket(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C, const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
+    int fmpz_mod_mpoly_compose_fmpz_mod_mpoly(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, fmpz_mod_mpoly_struct * const * C, const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
+    void fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const slong * c, const fmpz_mod_mpoly_ctx_t ctxB, const fmpz_mod_mpoly_ctx_t ctxAC)
+    void fmpz_mod_mpoly_mul(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_t C, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_mul_johnson(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_t C, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_mul_dense(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_t C, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_pow_fmpz(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_t k, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_pow_ui(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, ulong k, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_divides(fmpz_mod_mpoly_t Q, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_div(fmpz_mod_mpoly_t Q, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_divrem(fmpz_mod_mpoly_t Q, fmpz_mod_mpoly_t R, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_divrem_ideal(fmpz_mod_mpoly_struct ** Q, fmpz_mod_mpoly_t R, const fmpz_mod_mpoly_t A, fmpz_mod_mpoly_struct * const * B, slong len, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_term_content(fmpz_mod_mpoly_t M, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_content_vars(fmpz_mod_mpoly_t g, const fmpz_mod_mpoly_t A, slong * vars, slong vars_length, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd(fmpz_mod_mpoly_t G, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd_cofactors(fmpz_mod_mpoly_t G, fmpz_mod_mpoly_t Abar, fmpz_mod_mpoly_t Bbar, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd_brown(fmpz_mod_mpoly_t G, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd_hensel(fmpz_mod_mpoly_t G, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd_subresultant(fmpz_mod_mpoly_t G, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd_zippel(fmpz_mod_mpoly_t G, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_gcd_zippel2(fmpz_mod_mpoly_t G, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_resultant(fmpz_mod_mpoly_t R, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_discriminant(fmpz_mod_mpoly_t D, const fmpz_mod_mpoly_t A, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_sqrt(fmpz_mod_mpoly_t Q, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_is_square(const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_quadratic_root(fmpz_mod_mpoly_t Q, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_univar_init(fmpz_mod_mpoly_univar_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_univar_clear(fmpz_mod_mpoly_univar_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_univar_swap(fmpz_mod_mpoly_univar_t A, fmpz_mod_mpoly_univar_t B, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_to_univar(fmpz_mod_mpoly_univar_t A, const fmpz_mod_mpoly_t B, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_from_univar(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_univar_t B, slong var, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_univar_degree_fits_si(const fmpz_mod_mpoly_univar_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_univar_length(const fmpz_mod_mpoly_univar_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_univar_get_term_exp_si(fmpz_mod_mpoly_univar_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_univar_get_term_coeff(fmpz_mod_mpoly_t c, const fmpz_mod_mpoly_univar_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_univar_swap_term_coeff(fmpz_mod_mpoly_t c, fmpz_mod_mpoly_univar_t A, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_univar_set_coeff_ui(fmpz_mod_mpoly_univar_t Ax, ulong e, const fmpz_mod_mpoly_t c, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_univar_resultant(fmpz_mod_mpoly_t R, const fmpz_mod_mpoly_univar_t Ax, const fmpz_mod_mpoly_univar_t Bx, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_univar_discriminant(fmpz_mod_mpoly_t D, const fmpz_mod_mpoly_univar_t Ax, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_inflate(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_struct * shift, const fmpz_struct * stride, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_deflate(fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_t B, const fmpz_struct * shift, const fmpz_struct * stride, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_deflation(fmpz_struct * shift, fmpz_struct * stride, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
diff --git a/src/flint/flintlib/fmpz_mod_mpoly_factor.pxd b/src/flint/flintlib/fmpz_mod_mpoly_factor.pxd
new file mode 100644
index 00000000..aa9b9cf2
--- /dev/null
+++ b/src/flint/flintlib/fmpz_mod_mpoly_factor.pxd
@@ -0,0 +1,26 @@
+from flint.flintlib.flint cimport slong
+from flint.flintlib.fmpz cimport fmpz_t, fmpz_struct
+from flint.flintlib.fmpz_mod_mpoly cimport fmpz_mod_mpoly_ctx_t, fmpz_mod_mpoly_t, fmpz_mod_mpoly_struct
+
+
+cdef extern from "flint/fmpz_mod_mpoly_factor.h":
+    ctypedef struct fmpz_mod_mpoly_factor_struct:
+        fmpz_t constant
+        fmpz_mod_mpoly_struct * poly
+        fmpz_struct * exp
+        slong num
+        slong alloc
+
+    ctypedef fmpz_mod_mpoly_factor_struct fmpz_mod_mpoly_factor_t[1]
+
+    void fmpz_mod_mpoly_factor_init(fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_factor_clear(fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_factor_swap(fmpz_mod_mpoly_factor_t f, fmpz_mod_mpoly_factor_t g, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_factor_length(const fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_factor_get_constant_fmpz(fmpz_t c, const fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_factor_get_base(fmpz_mod_mpoly_t B, const fmpz_mod_mpoly_factor_t f, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_factor_swap_base(fmpz_mod_mpoly_t B, fmpz_mod_mpoly_factor_t f, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    slong fmpz_mod_mpoly_factor_get_exp_si(fmpz_mod_mpoly_factor_t f, slong i, const fmpz_mod_mpoly_ctx_t ctx)
+    void fmpz_mod_mpoly_factor_sort(fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_factor_squarefree(fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
+    int fmpz_mod_mpoly_factor(fmpz_mod_mpoly_factor_t f, const fmpz_mod_mpoly_t A, const fmpz_mod_mpoly_ctx_t ctx)
diff --git a/src/flint/flintlib/fmpz_mpoly_factor.pxd b/src/flint/flintlib/fmpz_mpoly_factor.pxd
index bba20a9b..62335351 100644
--- a/src/flint/flintlib/fmpz_mpoly_factor.pxd
+++ b/src/flint/flintlib/fmpz_mpoly_factor.pxd
@@ -4,8 +4,6 @@ from flint.flintlib.flint cimport slong, fmpz_struct
 from flint.flintlib.fmpq cimport fmpq_t
 
 
-# unimported types  set()
-
 cdef extern from "flint/fmpz_mpoly_factor.h":
 
     ctypedef struct fmpz_mpoly_factor_struct:
diff --git a/src/flint/flintlib/nmod.pxd b/src/flint/flintlib/nmod.pxd
index b748f90f..c88b876b 100644
--- a/src/flint/flintlib/nmod.pxd
+++ b/src/flint/flintlib/nmod.pxd
@@ -1,11 +1,7 @@
-from flint.flintlib.flint cimport mp_limb_t, mp_bitcnt_t, ulong
+from flint.flintlib.flint cimport mp_limb_t, mp_bitcnt_t, ulong, nmod_t
 from flint.flintlib.fmpz cimport fmpz_t
 
 cdef extern from "flint/nmod.h":
-    ctypedef struct nmod_t:
-        mp_limb_t n
-        mp_limb_t ninv
-        mp_bitcnt_t norm
 # TODO add macros
 
 # from here on is parsed
diff --git a/src/flint/flintlib/nmod_mpoly.pxd b/src/flint/flintlib/nmod_mpoly.pxd
new file mode 100644
index 00000000..4185c3b1
--- /dev/null
+++ b/src/flint/flintlib/nmod_mpoly.pxd
@@ -0,0 +1,159 @@
+from flint.flintlib.flint cimport fmpz_struct, flint_rand_t, mp_limb_t, slong, ulong, flint_bitcnt_t, nmod_t
+from flint.flintlib.fmpz cimport fmpz_t
+from flint.flintlib.mpoly cimport mpoly_ctx_t, ordering_t
+from flint.flintlib.nmod_poly cimport nmod_poly_struct, nmod_poly_t
+
+
+cdef extern from "flint/nmod_mpoly.h":
+    ctypedef struct nmod_mpoly_struct:
+        mp_limb_t * coeffs
+        ulong * exps
+        slong length
+        flint_bitcnt_t bits
+        slong coeffs_alloc
+        slong exps_alloc
+
+    ctypedef nmod_mpoly_struct nmod_mpoly_t[1]
+
+    ctypedef struct nmod_mpoly_ctx_struct:
+        mpoly_ctx_t minfo
+        nmod_t mod
+
+    ctypedef nmod_mpoly_ctx_struct nmod_mpoly_ctx_t[1]
+
+    ctypedef struct nmod_mpoly_univar_struct:
+        nmod_mpoly_struct * coeffs
+        fmpz_struct * exps
+        slong alloc
+        slong length
+
+    ctypedef nmod_mpoly_univar_struct nmod_mpoly_univar_t[1]
+
+# from here on is parsed
+    void nmod_mpoly_ctx_init(nmod_mpoly_ctx_t ctx, slong nvars, const ordering_t ord, mp_limb_t n)
+    slong nmod_mpoly_ctx_nvars(const nmod_mpoly_ctx_t ctx)
+    ordering_t nmod_mpoly_ctx_ord(const nmod_mpoly_ctx_t ctx)
+    mp_limb_t nmod_mpoly_ctx_modulus(const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_ctx_clear(nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_init(nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_init2(nmod_mpoly_t A, slong alloc, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_init3(nmod_mpoly_t A, slong alloc, flint_bitcnt_t bits, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_fit_length(nmod_mpoly_t A, slong len, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_realloc(nmod_mpoly_t A, slong alloc, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_clear(nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    char * nmod_mpoly_get_str_pretty(const nmod_mpoly_t A, const char ** x, const nmod_mpoly_ctx_t ctx)
+    # int nmod_mpoly_fprint_pretty(FILE * file, const nmod_mpoly_t A, const char ** x, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_print_pretty(const nmod_mpoly_t A, const char ** x, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_set_str_pretty(nmod_mpoly_t A, const char * str, const char ** x, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_gen(nmod_mpoly_t A, slong var, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_is_gen(const nmod_mpoly_t A, slong var, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_equal(const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_swap(nmod_mpoly_t A, nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_is_ui(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_get_ui(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_ui(nmod_mpoly_t A, ulong c, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_zero(nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_one(nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_equal_ui(const nmod_mpoly_t A, ulong c, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_is_zero(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_is_one(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_degrees_fit_si(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_degrees_fmpz(fmpz_struct ** degs, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_degrees_si(slong * degs, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_degree_fmpz(fmpz_t deg, const nmod_mpoly_t A, slong var, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_degree_si(const nmod_mpoly_t A, slong var, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_total_degree_fits_si(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_total_degree_fmpz(fmpz_t tdeg, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_total_degree_si(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_used_vars(int * used, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_get_coeff_ui_monomial(const nmod_mpoly_t A, const nmod_mpoly_t M, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_coeff_ui_monomial(nmod_mpoly_t A, ulong c, const nmod_mpoly_t M, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_get_coeff_ui_fmpz(const nmod_mpoly_t A, fmpz_struct * const * exp, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_get_coeff_ui_ui(const nmod_mpoly_t A, const ulong * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_coeff_ui_fmpz(nmod_mpoly_t A, ulong c, fmpz_struct * const * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_coeff_ui_ui(nmod_mpoly_t A, ulong c, const ulong * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_get_coeff_vars_ui(nmod_mpoly_t C, const nmod_mpoly_t A, const slong * vars, const ulong * exps, slong length, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_cmp(const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    mp_limb_t * nmod_mpoly_term_coeff_ref(nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_is_canonical(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_length(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_resize(nmod_mpoly_t A, slong new_length, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_get_term_coeff_ui(const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_term_coeff_ui(nmod_mpoly_t A, slong i, ulong c, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_term_exp_fits_si(const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_term_exp_fits_ui(const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_get_term_exp_fmpz(fmpz_struct ** exp, const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_get_term_exp_ui(ulong * exp, const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_get_term_exp_si(slong * exp, const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_get_term_var_exp_ui(const nmod_mpoly_t A, slong i, slong var, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_get_term_var_exp_si(const nmod_mpoly_t A, slong i, slong var, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_term_exp_fmpz(nmod_mpoly_t A, slong i, fmpz_struct * const * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_set_term_exp_ui(nmod_mpoly_t A, slong i, const ulong * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_get_term(nmod_mpoly_t M, const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_get_term_monomial(nmod_mpoly_t M, const nmod_mpoly_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_push_term_ui_fmpz(nmod_mpoly_t A, ulong c, fmpz_struct * const * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_push_term_ui_ffmpz(nmod_mpoly_t A, ulong c, const fmpz_struct * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_push_term_ui_ui(nmod_mpoly_t A, ulong c, const ulong * exp, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_sort_terms(nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_combine_like_terms(nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_reverse(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_randtest_bound(nmod_mpoly_t A, flint_rand_t state, slong length, ulong exp_bound, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_randtest_bounds(nmod_mpoly_t A, flint_rand_t state, slong length, ulong * exp_bounds, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_randtest_bits(nmod_mpoly_t A, flint_rand_t state, slong length, mp_limb_t exp_bits, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_add_ui(nmod_mpoly_t A, const nmod_mpoly_t B, ulong c, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_sub_ui(nmod_mpoly_t A, const nmod_mpoly_t B, ulong c, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_add(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_sub(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_neg(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_scalar_mul_ui(nmod_mpoly_t A, const nmod_mpoly_t B, ulong c, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_make_monic(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_derivative(nmod_mpoly_t A, const nmod_mpoly_t B, slong var, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_evaluate_all_ui(const nmod_mpoly_t A, const ulong * vals, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_evaluate_one_ui(nmod_mpoly_t A, const nmod_mpoly_t B, slong var, ulong val, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_compose_nmod_poly(nmod_poly_t A, const nmod_mpoly_t B, nmod_poly_struct * const * C, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_compose_nmod_mpoly_geobucket(nmod_mpoly_t A, const nmod_mpoly_t B, nmod_mpoly_struct * const * C, const nmod_mpoly_ctx_t ctxB, const nmod_mpoly_ctx_t ctxAC)
+    int nmod_mpoly_compose_nmod_mpoly_horner(nmod_mpoly_t A, const nmod_mpoly_t B, nmod_mpoly_struct * const * C, const nmod_mpoly_ctx_t ctxB, const nmod_mpoly_ctx_t ctxAC)
+    int nmod_mpoly_compose_nmod_mpoly(nmod_mpoly_t A, const nmod_mpoly_t B, nmod_mpoly_struct * const * C, const nmod_mpoly_ctx_t ctxB, const nmod_mpoly_ctx_t ctxAC)
+    void nmod_mpoly_compose_nmod_mpoly_gen(nmod_mpoly_t A, const nmod_mpoly_t B, const slong * c, const nmod_mpoly_ctx_t ctxB, const nmod_mpoly_ctx_t ctxAC)
+    void nmod_mpoly_mul(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_mul_johnson(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_mul_heap_threaded(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_mul_array(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_mul_array_threaded(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_mul_dense(nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_pow_fmpz(nmod_mpoly_t A, const nmod_mpoly_t B, const fmpz_t k, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_pow_ui(nmod_mpoly_t A, const nmod_mpoly_t B, ulong k, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_divides(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_div(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_divrem(nmod_mpoly_t Q, nmod_mpoly_t R, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_divrem_ideal(nmod_mpoly_struct ** Q, nmod_mpoly_t R, const nmod_mpoly_t A, nmod_mpoly_struct * const * B, slong len, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_divides_dense(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_divides_monagan_pearce(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_divides_heap_threaded(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_term_content(nmod_mpoly_t M, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_content_vars(nmod_mpoly_t g, const nmod_mpoly_t A, slong * vars, slong vars_length, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_gcd(nmod_mpoly_t G, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_gcd_cofactors(nmod_mpoly_t G, nmod_mpoly_t Abar, nmod_mpoly_t Bbar, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_gcd_brown(nmod_mpoly_t G, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_gcd_hensel(nmod_mpoly_t G, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_gcd_zippel(nmod_mpoly_t G, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_resultant(nmod_mpoly_t R, const nmod_mpoly_t A, const nmod_mpoly_t B, slong var, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_discriminant(nmod_mpoly_t D, const nmod_mpoly_t A, slong var, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_sqrt(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_is_square(const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_quadratic_root(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_univar_init(nmod_mpoly_univar_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_univar_clear(nmod_mpoly_univar_t A, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_univar_swap(nmod_mpoly_univar_t A, nmod_mpoly_univar_t B, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_to_univar(nmod_mpoly_univar_t A, const nmod_mpoly_t B, slong var, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_from_univar(nmod_mpoly_t A, const nmod_mpoly_univar_t B, slong var, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_univar_degree_fits_si(const nmod_mpoly_univar_t A, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_univar_length(const nmod_mpoly_univar_t A, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_univar_get_term_exp_si(nmod_mpoly_univar_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_univar_get_term_coeff(nmod_mpoly_t c, const nmod_mpoly_univar_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_univar_swap_term_coeff(nmod_mpoly_t c, nmod_mpoly_univar_t A, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_pow_rmul(nmod_mpoly_t A, const nmod_mpoly_t B, ulong k, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_div_monagan_pearce(nmod_mpoly_t polyq, const nmod_mpoly_t poly2, const nmod_mpoly_t poly3, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_divrem_monagan_pearce(nmod_mpoly_t q, nmod_mpoly_t r, const nmod_mpoly_t poly2, const nmod_mpoly_t poly3, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_divrem_ideal_monagan_pearce(nmod_mpoly_struct ** q, nmod_mpoly_t r, const nmod_mpoly_t poly2, nmod_mpoly_struct * const * poly3, slong len, const nmod_mpoly_ctx_t ctx)
diff --git a/src/flint/flintlib/nmod_mpoly_factor.pxd b/src/flint/flintlib/nmod_mpoly_factor.pxd
new file mode 100644
index 00000000..5e6fefcd
--- /dev/null
+++ b/src/flint/flintlib/nmod_mpoly_factor.pxd
@@ -0,0 +1,26 @@
+from flint.flintlib.nmod_mpoly cimport nmod_mpoly_ctx_t, nmod_mpoly_struct, nmod_mpoly_t
+from flint.flintlib.flint cimport slong, ulong, mp_limb_t
+from flint.flintlib.fmpz cimport fmpz_struct
+
+
+cdef extern from "flint/nmod_mpoly_factor.h":
+    ctypedef struct nmod_mpoly_factor_struct:
+        mp_limb_t constant
+        nmod_mpoly_struct * poly
+        fmpz_struct * exp
+        slong num
+        slong alloc
+
+    ctypedef nmod_mpoly_factor_struct nmod_mpoly_factor_t[1]
+
+    void nmod_mpoly_factor_init(nmod_mpoly_factor_t f, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_factor_clear(nmod_mpoly_factor_t f, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_factor_swap(nmod_mpoly_factor_t f, nmod_mpoly_factor_t g, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_factor_length(const nmod_mpoly_factor_t f, const nmod_mpoly_ctx_t ctx)
+    ulong nmod_mpoly_factor_get_constant_ui(const nmod_mpoly_factor_t f, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_factor_get_base(nmod_mpoly_t p, const nmod_mpoly_factor_t f, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_factor_swap_base(nmod_mpoly_t p, nmod_mpoly_factor_t f, slong i, const nmod_mpoly_ctx_t ctx)
+    slong nmod_mpoly_factor_get_exp_si(nmod_mpoly_factor_t f, slong i, const nmod_mpoly_ctx_t ctx)
+    void nmod_mpoly_factor_sort(nmod_mpoly_factor_t f, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_factor_squarefree(nmod_mpoly_factor_t f, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
+    int nmod_mpoly_factor(nmod_mpoly_factor_t f, const nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx)
diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py
index 323d4af5..783a4084 100644
--- a/src/flint/test/test_all.py
+++ b/src/flint/test/test_all.py
@@ -2750,19 +2750,52 @@ def setbad(obj, i, val):
 
 def _all_mpolys():
     return [
-        (flint.fmpz_mpoly, flint.fmpz_mpoly_ctx, flint.fmpz, False),
-        (flint.fmpq_mpoly, flint.fmpq_mpoly_ctx, flint.fmpq, True),
+        (flint.fmpz_mpoly, flint.fmpz_mpoly_ctx.get_context, flint.fmpz, False),
+        (flint.fmpq_mpoly, flint.fmpq_mpoly_ctx.get_context, flint.fmpq, True),
+        (
+            flint.fmpz_mod_mpoly,
+            lambda *args, **kwargs: flint.fmpz_mod_mpoly_ctx.get_context(*args, **kwargs, modulus=101),
+            lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(101)),
+            True,
+        ),
+        (
+            flint.fmpz_mod_mpoly,
+            lambda *args, **kwargs: flint.fmpz_mod_mpoly_ctx.get_context(*args, **kwargs, modulus=100),
+            lambda x: flint.fmpz_mod(x, flint.fmpz_mod_ctx(100)),
+            False,
+        ),
+        (
+            flint.nmod_mpoly,
+            lambda *args, **kwargs: flint.nmod_mpoly_ctx.get_context(*args, **kwargs, modulus=101),
+            lambda x: flint.nmod(x, 101),
+            True,
+        ),
+        (
+            flint.nmod_mpoly,
+            lambda *args, **kwargs: flint.nmod_mpoly_ctx.get_context(*args, **kwargs, modulus=100),
+            lambda x: flint.nmod(x, 100),
+            False,
+        ),
     ]
 
 
 def test_mpolys():
-    for P, C, S, is_field in _all_mpolys():
+    for P, get_context, S, is_field in _all_mpolys():
 
-        ctx = C.get_context(nvars=2)
+        # 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.
+        division_not_supported = P is not flint.fmpz_mpoly and not is_field
+        characteristic_zero = not (P is flint.fmpz_mod_mpoly or P is flint.nmod_mpoly)
 
-        assert raises(lambda: C.get_context(nvars=2, ordering="bad"), TypeError)
-        assert raises(lambda: C.get_context(nvars=-1), ValueError)
-        assert raises(lambda: C(-1, flint.Ordering.lex, []), ValueError)
+        ctx = get_context(nvars=2)
+
+        assert raises(lambda: get_context(nvars=2, ordering="bad"), TypeError)
+        assert raises(lambda: get_context(nvars=-1), ValueError)
+        if ctx.__class__ is flint.fmpz_mod_mpoly_ctx or ctx.__class__ is flint.nmod_mpoly_ctx:
+            assert raises(lambda: ctx.__class__(-1, flint.Ordering.lex, [], 4), ValueError)
+        else:
+            assert raises(lambda: ctx.__class__(-1, flint.Ordering.lex, []), ValueError)
         assert raises(lambda: ctx.constant("bad"), TypeError)
         assert raises(lambda: ctx.from_dict("bad"), ValueError)
         assert raises(lambda: ctx.from_dict({(0, 0): "bad"}), TypeError)
@@ -2771,7 +2804,7 @@ def test_mpolys():
         assert raises(lambda: ctx.gen(-1), IndexError)
         assert raises(lambda: ctx.gen(10), IndexError)
 
-        assert raises(lambda: P(val=C.get_context(nvars=1).constant(0), ctx=ctx), IncompatibleContextError)
+        assert raises(lambda: P(val=get_context(nvars=1).constant(0), ctx=ctx), IncompatibleContextError)
         assert raises(lambda: P(val={}, ctx=None), ValueError)
         assert raises(lambda: P(val={"bad": 1}, ctx=None), ValueError)
         assert raises(lambda: P(val="1", ctx=None), ValueError)
@@ -2790,10 +2823,10 @@ def quick_poly():
         assert ctx.nvars() == 2
         assert ctx.ordering() == flint.Ordering.lex
 
-        ctx1 = C.get_context(4)
+        ctx1 = get_context(4)
         assert [ctx1.name(i) for i in range(4)] == ['x0', 'x1', 'x2', 'x3']
         for order in list(flint.Ordering):
-            ctx1 = C.get_context(4, order)
+            ctx1 = get_context(4, order)
             assert ctx1.ordering() == order
 
         assert ctx.constant(1) == mpoly({(0, 0): 1}) == P(1, ctx=ctx)
@@ -2822,10 +2855,15 @@ def quick_poly():
         assert (P(1, ctx=ctx) == P(2, ctx=ctx)) is False
         assert (P(1, ctx=ctx) != P(2, ctx=ctx)) is True
 
-        assert (P(1, ctx=ctx) == 1) is False
-        assert (P(1, ctx=ctx) != 1) is True
-        assert (1 == P(1, ctx=ctx)) is False
-        assert (1 != P(1, ctx=ctx)) is True
+        assert (P(1, ctx=ctx) == 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
@@ -2930,7 +2968,7 @@ def quick_poly():
         assert raises(lambda: p.subs({"a": 1}), ValueError)
         assert raises(lambda: p.subs({"x0": 0, "x1": 1, "x2": 2}), ValueError)
 
-        no_gens_ctx = C.get_context(0)
+        no_gens_ctx = get_context(0)
         no_gens_p = P("2", no_gens_ctx)
         assert no_gens_p.compose(ctx=ctx1).context() is ctx1
         assert raises(lambda: no_gens_p.compose(), ValueError)
@@ -2948,8 +2986,8 @@ def quick_poly():
 
         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): -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}) \
@@ -3013,27 +3051,53 @@ def quick_poly():
         assert raises(lambda: quick_poly().imul(P(ctx=ctx1)), IncompatibleContextError)
         assert raises(lambda: quick_poly().imul(None), NotImplementedError)
 
-        assert quick_poly() // mpoly({(1, 1): 1}) == mpoly({(1, 1): 4})
-        assert quick_poly() % mpoly({(1, 1): 1}) \
-            == mpoly({(1, 0): 3, (0, 1): 2, (0, 0): 1})
-        assert divmod(quick_poly(), mpoly({(1, 1): 1})) \
-            == (mpoly({(1, 1): 4}), mpoly({(1, 0): 3, (0, 1): 2, (0, 0): 1}))
-
-        assert 1 / P(1, ctx=ctx) == P(1, ctx=ctx)
-        assert quick_poly() / 1 == quick_poly()
-        assert quick_poly() // 1 == quick_poly()
-        assert quick_poly() % 1 == P(ctx=ctx)
-        assert divmod(quick_poly(), 1) == (quick_poly(), P(ctx=ctx))
+        if division_not_supported:
+            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)})
+            assert quick_poly() / 3 == mpoly({(0, 0): S(1) / 3, (0, 1): S(2) / 3, (1, 0): S(1), (2, 2): S(4) / 3})
         else:
             assert raises(lambda: quick_poly() / 3, DomainError)
 
-        assert 1 // quick_poly() == P(ctx=ctx)
-        assert 1 % quick_poly() == P(1, ctx=ctx)
-        assert divmod(1, quick_poly()) == (P(ctx=ctx), P(1, ctx=ctx))
+        f = mpoly({(1, 1): 4, (0, 0): 1})
+        g = mpoly({(0, 1): 2, (1, 0): 2})
+        if not division_not_supported:
+            assert 1 // quick_poly() == P(ctx=ctx)
+            assert 1 % quick_poly() == P(1, ctx=ctx)
+            assert divmod(1, quick_poly()) == (P(ctx=ctx), P(1, ctx=ctx))
+
+            assert 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 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)
         assert raises(lambda: quick_poly() // None, TypeError)
         assert raises(lambda: quick_poly() % None, TypeError)
@@ -3064,15 +3128,6 @@ def quick_poly():
         assert raises(lambda: quick_poly() % P(1, ctx=ctx1), IncompatibleContextError)
         assert raises(lambda: divmod(quick_poly(), P(1, ctx=ctx1)), IncompatibleContextError)
 
-        f = mpoly({(1, 1): 4, (0, 0): 1})
-        g = mpoly({(0, 1): 2, (1, 0): 2})
-        assert f * g / mpoly({(0, 1): 1, (1, 0): 1}) \
-            == mpoly({(1, 1): 8, (0, 0): 2})
-
-        if not is_field:
-            assert raises(lambda: 1 / quick_poly(), DomainError)
-            assert raises(lambda: quick_poly() / P(2, ctx=ctx), DomainError)
-
         assert quick_poly() ** 0 == P(1, ctx=ctx)
         assert quick_poly() ** 1 == quick_poly()
         assert quick_poly() ** 2 == mpoly({
@@ -3093,19 +3148,34 @@ def quick_poly():
         # # XXX: Not sure what this should do in general:
         assert raises(lambda: pow(P(1, ctx=ctx), 2, 3), NotImplementedError)
 
-        if is_field:
-            assert (f * g).gcd(f) == f / 4
+        if division_not_supported:
+            assert raises(lambda: (f * g).gcd(f), DomainError)
         else:
-            assert (f * g).gcd(f) == f
-        assert raises(lambda: quick_poly().gcd(None), TypeError)
-        assert raises(lambda: quick_poly().gcd(P(ctx=ctx1)), IncompatibleContextError)
-
-        assert (f * g).factor() == (S(2), [(mpoly({(0, 1): 1, (1, 0): 1}), 1), (f, 1)])
-
-        assert (f * f).sqrt() == f
-        if P is flint.fmpz_mpoly:
-            assert (f * f).sqrt(assume_perfect_square=True) == f
-        assert raises(lambda: quick_poly().sqrt(), ValueError)
+            if is_field:
+                assert (f * g).gcd(f) == f / 4
+            else:
+                assert (f * g).gcd(f) == f
+            assert raises(lambda: quick_poly().gcd(None), TypeError)
+            assert raises(lambda: quick_poly().gcd(P(ctx=ctx1)), IncompatibleContextError)
+
+        if division_not_supported:
+            # 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_zero:
+            # 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 division_not_supported:
+            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(), ValueError)
 
         p = quick_poly()
         assert p.derivative(0) == p.derivative("x0") == mpoly({(0, 0): 3, (1, 2): 8})
@@ -3115,20 +3185,21 @@ def quick_poly():
         assert raises(lambda: p.derivative(3), IndexError)
         assert raises(lambda: p.derivative(None), TypeError)
 
-        if is_field:
-            assert quick_poly().integral(0) == quick_poly().integral("x0") == \
-                mpoly({(3, 2): S(4, 3), (2, 0): S(3, 2), (1, 1): S(2), (1, 0): S(1)})
-            assert quick_poly().integral(1) == quick_poly().integral("x1") == \
-                mpoly({(2, 3): S(4, 3), (1, 1): S(3), (0, 2): S(1), (0, 1): S(1)})
-        else:
-            assert quick_poly().integral(0) == quick_poly().integral("x0") == \
-                (6, mpoly({(3, 2): 8, (2, 0): 9, (1, 1): 12, (1, 0): 6}))
-            assert quick_poly().integral(1) == quick_poly().integral("x1") == \
-                (3, mpoly({(2, 3): 4, (1, 1): 9, (0, 2): 3, (0, 1): 3}))
-
-        assert raises(lambda: p.integral("x3"), ValueError)
-        assert raises(lambda: p.integral(3), IndexError)
-        assert raises(lambda: p.integral(None), TypeError)
+        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)
 
 
 def _all_mpoly_vecs():
diff --git a/src/flint/types/fmpq_mpoly.pyx b/src/flint/types/fmpq_mpoly.pyx
index bcb3ea20..88219b54 100644
--- a/src/flint/types/fmpq_mpoly.pyx
+++ b/src/flint/types/fmpq_mpoly.pyx
@@ -28,6 +28,8 @@ from flint.flintlib.fmpq_mpoly cimport (
     fmpq_mpoly_divides,
     fmpq_mpoly_divrem,
     fmpq_mpoly_equal,
+    fmpq_mpoly_equal_fmpq,
+    fmpq_mpoly_equal_fmpz,
     fmpq_mpoly_evaluate_all_fmpq,
     fmpq_mpoly_evaluate_one_fmpq,
     fmpq_mpoly_gcd,
@@ -256,11 +258,18 @@ cdef class fmpq_mpoly(flint_mpoly):
             return op == Py_NE
         elif typecheck(self, fmpq_mpoly) and typecheck(other, fmpq_mpoly):
             if (<fmpq_mpoly>self).ctx is (<fmpq_mpoly>other).ctx:
-                return (op == Py_NE) ^ bool(
-                    fmpq_mpoly_equal((<fmpq_mpoly>self).val, (<fmpq_mpoly>other).val, (<fmpq_mpoly>self).ctx.val)
-                )
+                return (op == Py_NE) ^ <bint>fmpq_mpoly_equal(self.val, (<fmpq_mpoly>other).val, self.ctx.val)
             else:
                 return op == Py_NE
+        elif typecheck(other, fmpq):
+            return (op == Py_NE) ^ <bint>fmpq_mpoly_equal_fmpq(self.val, (<fmpq>other).val, self.ctx.val)
+        elif typecheck(other, fmpz):
+            return (op == Py_NE) ^ <bint>fmpq_mpoly_equal_fmpz(self.val, (<fmpz>other).val, self.ctx.val)
+        elif isinstance(other, int):
+            other = any_as_fmpz(other)
+            if other is NotImplemented:
+                return NotImplemented
+            return (op == Py_NE) ^ <bint>fmpq_mpoly_equal_fmpz(self.val, (<fmpz>other).val, self.ctx.val)
         else:
             return NotImplemented
 
@@ -533,17 +542,10 @@ cdef class fmpq_mpoly(flint_mpoly):
             fmpq vres
             slong nvars = self.ctx.nvars(), nargs = len(args)
 
-        if nargs < nvars:
-            raise ValueError("not enough arguments provided")
-        elif nargs > nvars:
-            raise ValueError("more arguments provided than variables")
-
-        args_fmpq = tuple(any_as_fmpq(v) for v in args)
-        for arg in args_fmpq:
-            if arg is NotImplemented:
-                raise TypeError(f"cannot coerce argument ('{arg}') to fmpq")
+        if nargs != nvars:
+            raise ValueError("number of generators does not match number of arguments")
 
-        V = fmpq_vec(args_fmpq, double_indirect=True)
+        V = fmpq_vec(args, double_indirect=True)
         vres = fmpq.__new__(fmpq)
         if fmpq_mpoly_evaluate_all_fmpq(vres.val, self.val, V.double_indirect, self.ctx.val) == 0:
             raise ValueError("unreasonably large polynomial")  # pragma: no cover
@@ -916,7 +918,7 @@ cdef class fmpq_mpoly(flint_mpoly):
             >>> Zm = fmpq_mpoly
             >>> ctx = fmpq_mpoly_ctx.get_context(3, Ordering.lex, 'x,y,z')
             >>> p1 = Zm("2*x + 4", ctx)
-            >>> p2 = Zm("3*x*z +  + 3*x + 3*z + 3", ctx)
+            >>> p2 = Zm("3*x*z + 3*x + 3*z + 3", ctx)
             >>> (p1 * p2).factor()
             (6, [(z + 1, 1), (x + 2, 1), (x + 1, 1)])
             >>> (p2 * p1 * p2).factor()
@@ -927,7 +929,8 @@ cdef class fmpq_mpoly(flint_mpoly):
             fmpq_mpoly u
 
         fmpq_mpoly_factor_init(fac, self.ctx.val)
-        fmpq_mpoly_factor(fac, self.val, self.ctx.val)
+        if not fmpq_mpoly_factor(fac, self.val, self.ctx.val):
+            raise RuntimeError("factorisation failed")
         res = [0] * fac.num
 
         for i in range(fac.num):
diff --git a/src/flint/types/fmpz_mod_mpoly.pxd b/src/flint/types/fmpz_mod_mpoly.pxd
new file mode 100644
index 00000000..3f259c29
--- /dev/null
+++ b/src/flint/types/fmpz_mod_mpoly.pxd
@@ -0,0 +1,37 @@
+from flint.flint_base.flint_base cimport flint_mpoly, flint_mpoly_context
+
+from flint.flintlib.fmpz_mod_mpoly cimport (
+    fmpz_mod_mpoly_ctx_t,
+    fmpz_mod_mpoly_t,
+    fmpz_mod_mpoly_init,
+    fmpz_mod_mpoly_struct
+)
+from flint.flintlib.flint cimport slong
+
+cdef inline init_fmpz_mod_mpoly(fmpz_mod_mpoly var, fmpz_mod_mpoly_ctx ctx):
+    var.ctx = ctx
+    fmpz_mod_mpoly_init(var.val, ctx.val)
+    var._init = True
+
+cdef inline fmpz_mod_mpoly create_fmpz_mod_mpoly(fmpz_mod_mpoly_ctx ctx):
+    cdef fmpz_mod_mpoly var
+    var = fmpz_mod_mpoly.__new__(fmpz_mod_mpoly)
+    var.ctx = ctx
+    fmpz_mod_mpoly_init(var.val, ctx.val)
+    var._init = True
+    return var
+
+cdef class fmpz_mod_mpoly_ctx(flint_mpoly_context):
+    cdef fmpz_mod_mpoly_ctx_t val
+    cdef readonly object __prime_modulus
+
+cdef class fmpz_mod_mpoly(flint_mpoly):
+    cdef fmpz_mod_mpoly_t val
+    cdef fmpz_mod_mpoly_ctx ctx
+    cdef bint _init
+
+cdef class fmpz_mod_mpoly_vec:
+    cdef fmpz_mod_mpoly_struct *val
+    cdef slong length
+    cdef fmpz_mod_mpoly_ctx ctx
+    cdef fmpz_mod_mpoly_struct **double_indirect
diff --git a/src/flint/types/fmpz_mod_mpoly.pyx b/src/flint/types/fmpz_mod_mpoly.pyx
new file mode 100644
index 00000000..43f5b328
--- /dev/null
+++ b/src/flint/types/fmpz_mod_mpoly.pyx
@@ -0,0 +1,1029 @@
+from flint.flint_base.flint_base cimport (
+    flint_mpoly,
+    flint_mpoly_context,
+    ordering_py_to_c,
+    ordering_c_to_py,
+)
+from flint.flint_base.flint_base import Ordering
+
+from flint.utils.typecheck cimport typecheck
+from flint.utils.flint_exceptions import DomainError, IncompatibleContextError
+
+from flint.types.fmpz cimport any_as_fmpz, fmpz
+from flint.types.fmpz_vec cimport fmpz_vec
+from flint.types.fmpz_mod cimport fmpz_mod
+from flint.types.nmod cimport nmod
+
+from flint.flintlib.fmpz cimport fmpz_set
+from flint.flintlib.fmpz_mod_mpoly cimport (
+    fmpz_mod_mpoly_add,
+    fmpz_mod_mpoly_add_fmpz,
+    fmpz_mod_mpoly_clear,
+    fmpz_mod_mpoly_compose_fmpz_mod_mpoly,
+    fmpz_mod_mpoly_ctx_get_modulus,
+    fmpz_mod_mpoly_ctx_init,
+    fmpz_mod_mpoly_degrees_fmpz,
+    fmpz_mod_mpoly_derivative,
+    fmpz_mod_mpoly_div,
+    fmpz_mod_mpoly_divides,
+    fmpz_mod_mpoly_divrem,
+    fmpz_mod_mpoly_equal,
+    fmpz_mod_mpoly_equal_fmpz,
+    fmpz_mod_mpoly_evaluate_all_fmpz,
+    fmpz_mod_mpoly_evaluate_one_fmpz,
+    fmpz_mod_mpoly_gcd,
+    fmpz_mod_mpoly_gen,
+    fmpz_mod_mpoly_get_coeff_fmpz_fmpz,
+    fmpz_mod_mpoly_get_str_pretty,
+    fmpz_mod_mpoly_get_term,
+    fmpz_mod_mpoly_get_term_coeff_fmpz,
+    fmpz_mod_mpoly_get_term_exp_fmpz,
+    fmpz_mod_mpoly_is_one,
+    fmpz_mod_mpoly_is_zero,
+    fmpz_mod_mpoly_length,
+    fmpz_mod_mpoly_mul,
+    fmpz_mod_mpoly_neg,
+    fmpz_mod_mpoly_pow_fmpz,
+    fmpz_mod_mpoly_push_term_fmpz_ffmpz,
+    fmpz_mod_mpoly_scalar_mul_fmpz,
+    fmpz_mod_mpoly_set,
+    fmpz_mod_mpoly_set_coeff_fmpz_fmpz,
+    fmpz_mod_mpoly_set_fmpz,
+    fmpz_mod_mpoly_set_str_pretty,
+    fmpz_mod_mpoly_sort_terms,
+    fmpz_mod_mpoly_sub,
+    fmpz_mod_mpoly_sub_fmpz,
+    fmpz_mod_mpoly_total_degree_fmpz,
+    fmpz_mod_mpoly_sqrt,
+)
+from flint.flintlib.fmpz_mod_mpoly_factor cimport (
+    fmpz_mod_mpoly_factor,
+    fmpz_mod_mpoly_factor_clear,
+    fmpz_mod_mpoly_factor_init,
+    fmpz_mod_mpoly_factor_squarefree,
+    fmpz_mod_mpoly_factor_t,
+)
+
+from cpython.object cimport Py_EQ, Py_NE
+cimport libc.stdlib
+
+cdef dict _fmpz_mod_mpoly_ctx_cache = {}
+
+
+cdef class fmpz_mod_mpoly_ctx(flint_mpoly_context):
+    """
+    A class for storing the polynomial context
+
+    :param nvars: The number of variables in the ring
+    :param ordering:  The term order for the ring
+    :param names:  A tuple containing the names of the variables of the ring.
+
+    Do not construct one of these directly, use `fmpz_mod_mpoly_ctx.get_context`.
+    """
+
+    _ctx_cache = _fmpz_mod_mpoly_ctx_cache
+
+    def __init__(self, slong nvars, ordering, names, modulus):
+        cdef fmpz m
+        if not typecheck(modulus, fmpz):
+            m = any_as_fmpz(modulus)
+            if m is NotImplemented:
+                raise TypeError(f"modulus ({modulus}) is not coercible to fmpz")
+        else:
+            m = modulus
+        fmpz_mod_mpoly_ctx_init(self.val, nvars, ordering_py_to_c(ordering), m.val)
+        self.__prime_modulus = None
+        super().__init__(nvars, names)
+
+    @classmethod
+    def create_context_key(
+            cls,
+            slong nvars=1,
+            ordering=Ordering.lex,
+            modulus = None,
+            names: Optional[str] = "x",
+            nametup: Optional[tuple] = None,
+    ):
+        """
+        Create a key for the context cache via the number of variables, the ordering, the modulus, and either a
+        variable name string, or a tuple of variable names.
+        """
+        # A type hint of `ordering: Ordering` results in the error "TypeError: an integer is required" if a Ordering
+        # object is not provided. This is pretty obtuse so we check its type ourselves
+        if not isinstance(ordering, Ordering):
+            raise TypeError(f"`ordering` ('{ordering}') is not an instance of flint.Ordering")
+        elif not typecheck(modulus, fmpz):
+            m = any_as_fmpz(modulus)
+            if m is NotImplemented:
+                raise TypeError(f"`modulus` ('{modulus}') is not coercible to fmpz")
+            else:
+                modulus = m
+
+        if nametup is not None:
+            key = nvars, ordering, nametup, modulus
+        elif nametup is None and names is not None:
+            key = nvars, ordering, cls.create_variable_names(nvars, names), modulus
+        else:
+            raise ValueError("must provide either `names` or `nametup`")
+        return key
+
+    def any_as_scalar(self, other):
+        if isinstance(other, int):
+            return any_as_fmpz(other)
+        elif typecheck(other, nmod):
+            if (<nmod>other).modulus() != self.modulus():
+                raise DomainError(
+                    f"modulus does not match, got {(<nmod>other).modulus()}, expected {self.modulus()}"
+                )
+            return any_as_fmpz((<nmod>other).val)
+        elif typecheck(other, fmpz):
+            return fmpz(other)
+        elif typecheck(other, fmpz_mod):
+            if (<fmpz_mod>other).ctx.modulus() != self.modulus():
+                raise DomainError(
+                    f"modulus does not match, got {(<fmpz_mod>other).ctx.modulus()}, expected {self.modulus()}"
+                )
+            res = fmpz.__new__(fmpz)
+            fmpz_set((<fmpz>res).val, (<fmpz_mod>other).val)
+            return res
+        else:
+            return NotImplemented
+
+    def scalar_as_mpoly(self, other: fmpz):
+        # non-fmpz scalars should first be converted via cls.any_as_scalar
+        return self.constant(<fmpz>other)
+
+    def nvars(self):
+        """
+        Return the number of variables in the context
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.lex, 11, 'x')
+            >>> ctx.nvars()
+            4
+        """
+        return self.val.minfo.nvars
+
+    def ordering(self):
+        """
+        Return the term order of the context object.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.deglex, 11, 'w')
+            >>> ctx.ordering()
+            <Ordering.deglex: 1>
+        """
+        return ordering_c_to_py(self.val.minfo.ord)
+
+    def modulus(self):
+        """
+        Return the modulus of the context object.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.deglex, 2, 'w')
+            >>> ctx.modulus()
+            2
+
+        """
+        cdef fmpz m = fmpz.__new__(fmpz)
+        fmpz_mod_mpoly_ctx_get_modulus(m.val, self.val)
+        return m
+
+    def is_prime(self):
+        """
+        Return whether the modulus is prime
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.degrevlex, 2**127, 'z')
+            >>> ctx.is_prime()
+            False
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.degrevlex, 2**127 - 1, 'z')
+            >>> ctx.is_prime()
+            True
+        """
+        if self.__prime_modulus is None:
+            self.__prime_modulus = self.modulus().is_prime()
+        return self.__prime_modulus
+
+    def gen(self, slong i):
+        """
+        Return the `i`th generator of the polynomial ring
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(3, Ordering.degrevlex, 11, 'z')
+            >>> ctx.gen(1)
+            z1
+        """
+        cdef fmpz_mod_mpoly res
+        if not 0 <= i < self.val.minfo.nvars:
+            raise IndexError("generator index out of range")
+        res = create_fmpz_mod_mpoly(self)
+        fmpz_mod_mpoly_gen(res.val, i, res.ctx.val)
+        return res
+
+    def constant(self, z):
+        """
+        Create a constant polynomial in this context
+        """
+        cdef fmpz_mod_mpoly res
+        z = any_as_fmpz(z)
+        if z is NotImplemented:
+            raise TypeError("cannot coerce argument to fmpz")
+        res = create_fmpz_mod_mpoly(self)
+        fmpz_mod_mpoly_set_fmpz(res.val, (<fmpz>z).val, res.ctx.val)
+        return res
+
+    def from_dict(self, d):
+        """
+        Create a fmpz_mod_mpoly from a dictionary in this context.
+
+        The dictionary's keys are tuples of ints (or anything that implicitly converts
+        to fmpz) representing exponents, and corresponding values of fmpz.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x,y')
+            >>> ctx.from_dict({(1,0):2, (1,1):3, (0,1):1})
+            3*x*y + 2*x + y
+        """
+        cdef:
+            fmpz_vec exp_vec
+            slong i, nvars = self.nvars()
+            fmpz_mod_mpoly res
+
+        if not isinstance(d, dict):
+            raise ValueError("expected a dictionary")
+
+        res = create_fmpz_mod_mpoly(self)
+
+        for i, (exps, coeff) in enumerate(d.items()):
+            if len(exps) != nvars:
+                raise ValueError(f"expected {nvars} exponents, got {len(exps)}")
+            elif not coeff:
+                continue
+
+            exp_vec = fmpz_vec(exps)
+            coeff_scalar = self.any_as_scalar(coeff)
+            if coeff_scalar is NotImplemented:
+                raise TypeError(f"cannot coerce {repr(coeff)} to nmod_mpoly coefficient")
+
+            fmpz_mod_mpoly_push_term_fmpz_ffmpz(
+                res.val,
+                (<fmpz>coeff_scalar).val,
+                exp_vec.val,
+                self.val
+            )
+
+        fmpz_mod_mpoly_sort_terms(res.val, self.val)
+        return res
+
+
+cdef class fmpz_mod_mpoly(flint_mpoly):
+    """
+    The *fmpz_mod_mpoly* type represents sparse multivariate polynomials over
+    the integers modulo `n`, for large `n`.
+    """
+
+    def __cinit__(self):
+        self._init = False
+
+    def __dealloc__(self):
+        if self._init:
+            fmpz_mod_mpoly_clear(self.val, self.ctx.val)
+            self._init = False
+
+    def __init__(self, val=0, ctx=None):
+        if typecheck(val, fmpz_mod_mpoly):
+            if ctx is None or ctx == (<fmpz_mod_mpoly>val).ctx:
+                init_fmpz_mod_mpoly(self, (<fmpz_mod_mpoly>val).ctx)
+                fmpz_mod_mpoly_set(self.val, (<fmpz_mod_mpoly>val).val, self.ctx.val)
+            else:
+                raise IncompatibleContextError(f"{ctx} is not {(<fmpz_mod_mpoly>val).ctx}")
+        elif isinstance(val, dict):
+            if ctx is None:
+                raise ValueError("a context is required to create a fmpz_mod_mpoly from a dict")
+            x = ctx.from_dict(val)
+            # XXX: this copy is silly, have a ctx function that assigns an fmpz_mod_mpoly_t
+            init_fmpz_mod_mpoly(self, ctx)
+            fmpz_mod_mpoly_set(self.val, (<fmpz_mod_mpoly>x).val, self.ctx.val)
+        elif isinstance(val, str):
+            if ctx is None:
+                raise ValueError("cannot parse a polynomial without context")
+            val = bytes(val, 'utf-8')
+            init_fmpz_mod_mpoly(self, ctx)
+            if fmpz_mod_mpoly_set_str_pretty(self.val, val, self.ctx.c_names, self.ctx.val) == -1:
+                raise ValueError("unable to parse fmpz_mod_mpoly from string")
+            fmpz_mod_mpoly_sort_terms(self.val, self.ctx.val)
+        elif typecheck(val, fmpz_mod):
+            if ctx is None:
+                raise ValueError("need context to convert fmpz_mod to fmpz_mod_mpoly")
+            init_fmpz_mod_mpoly(self, ctx)
+            fmpz_mod_mpoly_set_fmpz(self.val, (<fmpz_mod>val).val, self.ctx.val)
+        else:
+            v = any_as_fmpz(val)
+            if v is NotImplemented:
+                raise TypeError("cannot create fmpz_mod_mpoly from type %s" % type(val))
+            elif ctx is None:
+                raise ValueError("need context to convert fmpz to fmpz_mod_mpoly")
+            init_fmpz_mod_mpoly(self, ctx)
+            fmpz_mod_mpoly_set_fmpz(self.val, (<fmpz>v).val, self.ctx.val)
+
+    def _division_check(self, other):
+        super()._division_check(other)
+        if not self.ctx.is_prime():
+            raise DomainError("division with non-prime modulus is not supported")
+
+    def __bool__(self):
+        return not fmpz_mod_mpoly_is_zero(self.val, self.ctx.val)
+
+    def is_zero(self):
+        return <bint>fmpz_mod_mpoly_is_zero(self.val, self.ctx.val)
+
+    def is_one(self):
+        return <bint>fmpz_mod_mpoly_is_one(self.val, self.ctx.val)
+
+    def __richcmp__(self, other, int op):
+        if not (op == Py_EQ or op == Py_NE):
+            return NotImplemented
+        elif other is None:
+            return op == Py_NE
+        elif typecheck(self, fmpz_mod_mpoly) and typecheck(other, fmpz_mod_mpoly):
+            if (<fmpz_mod_mpoly>self).ctx is (<fmpz_mod_mpoly>other).ctx:
+                return (op == Py_NE) ^ bool(
+                    fmpz_mod_mpoly_equal((<fmpz_mod_mpoly>self).val, (<fmpz_mod_mpoly>other).val, (<fmpz_mod_mpoly>self).ctx.val)
+                )
+            else:
+                return op == Py_NE
+        elif typecheck(other, fmpz):
+            return (op == Py_NE) ^ <bint>fmpz_mod_mpoly_equal_fmpz(self.val, (<fmpz>other).val, self.ctx.val)
+        elif typecheck(other, fmpz_mod):
+            return (op == Py_NE) ^ <bint>fmpz_mod_mpoly_equal_fmpz(self.val, (<fmpz_mod>other).val, self.ctx.val)
+        elif isinstance(other, int):
+            other = any_as_fmpz(other)
+            if other is NotImplemented:
+                return NotImplemented
+            return (op == Py_NE) ^ <bint>fmpz_mod_mpoly_equal_fmpz(self.val, (<fmpz>other).val, self.ctx.val)
+        else:
+            return NotImplemented
+
+    def __len__(self):
+        return fmpz_mod_mpoly_length(self.val, self.ctx.val)
+
+    def __getitem__(self, x):
+        """
+        Return the coefficient of the term with the exponent vector `x`.
+        Always returns a value, missing keys will return `0`.
+        Negative exponents are made positive.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p[1, 1]
+            3
+
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        if not isinstance(x, tuple):
+            raise TypeError("exponent vector index is not a tuple")
+        elif len(x) != nvars:
+            raise ValueError("exponent vector provided does not match number of variables")
+
+        res = fmpz()
+        exp_vec = fmpz_vec(x, double_indirect=True)
+        fmpz_mod_mpoly_get_coeff_fmpz_fmpz((<fmpz>res).val, self.val, exp_vec.double_indirect, self.ctx.val)
+        return res
+
+    def __setitem__(self, x, y):
+        """
+        Set the coefficient of the term with the exponent vector `x` to `y`.
+        Will always set a value, missing keys will create a new term.
+        Negative exponents are made positive.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p[1, 1] = 20
+            >>> p
+            20*x0*x1 + 2*x1
+
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        if not isinstance(x, tuple):
+            raise TypeError("exponent vector index is not a tuple")
+        elif len(x) != nvars:
+            raise ValueError("exponent vector provided does not match number of variables")
+        exp_vec = fmpz_vec(x, double_indirect=True)
+
+        coeff = self.ctx.any_as_scalar(y)
+        if coeff is NotImplemented:
+            raise TypeError("provided coefficient not coercible to fmpz")
+        fmpz_mod_mpoly_set_coeff_fmpz_fmpz(self.val, (<fmpz>coeff).val, exp_vec.double_indirect, self.ctx.val)
+
+    def __neg__(self):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_neg(res.val, (<fmpz_mod_mpoly>self).val, res.ctx.val)
+        return res
+
+    def _add_scalar_(self, other: fmpz):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_add_fmpz(res.val, self.val, other.val, self.ctx.val)
+        return res
+
+    def _add_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_add(res.val, self.val, other.val, res.ctx.val)
+        return res
+
+    def _sub_scalar_(self, other: fmpz):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_sub_fmpz(res.val, self.val, other.val, self.ctx.val)
+        return res
+
+    def _sub_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_sub(res.val, self.val, other.val, res.ctx.val)
+        return res
+
+    def _mul_scalar_(self, other: fmpz):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_scalar_mul_fmpz(res.val, self.val, other.val, self.ctx.val)
+        return res
+
+    def _mul_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_mul(res.val, self.val, other.val, res.ctx.val)
+        return res
+
+    def _pow_(self, other: fmpz):
+        cdef fmpz_mod_mpoly res
+        res = create_fmpz_mod_mpoly(self.ctx)
+        if fmpz_mod_mpoly_pow_fmpz(res.val, self.val, other.val, res.ctx.val) == 0:
+            raise ValueError("unreasonably large polynomial")  # pragma: no cover
+        return res
+
+    def _divmod_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly quotient, remainder
+        quotient = create_fmpz_mod_mpoly(self.ctx)
+        remainder = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_divrem(quotient.val, remainder.val, self.val, other.val, self.ctx.val)
+        return (quotient, remainder)
+
+    def _floordiv_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly quotient
+        quotient = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_div(quotient.val, self.val, other.val, self.ctx.val)
+        return quotient
+
+    def _truediv_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly quotient
+        quotient = create_fmpz_mod_mpoly(self.ctx)
+        if fmpz_mod_mpoly_divides(quotient.val, self.val, other.val, self.ctx.val):
+            return quotient
+        else:
+            raise DomainError("fmpz_mod_mpoly division is not exact")
+
+    def _mod_mpoly_(self, other: fmpz_mod_mpoly):
+        cdef fmpz_mod_mpoly quotient, remainder
+        quotient = create_fmpz_mod_mpoly(self.ctx)
+        remainder = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_divrem(quotient.val, remainder.val, self.val, other.val, self.ctx.val)
+        return remainder
+
+    def __call__(self, *args) -> fmpz:
+        cdef:
+            fmpz_vec V
+            fmpz vres
+            slong nvars = self.ctx.nvars(), nargs = len(args)
+
+        if nargs != nvars:
+            raise ValueError("number of generators does not match number of arguments")
+        args = [self.ctx.any_as_scalar(x) for x in args]
+        V = fmpz_vec(args, double_indirect=True)
+        vres = fmpz.__new__(fmpz)
+        fmpz_mod_mpoly_evaluate_all_fmpz(vres.val, self.val, V.double_indirect, self.ctx.val)
+        return vres
+
+    def _iadd_scalar_(self, other: fmpz):
+        fmpz_mod_mpoly_add_fmpz(self.val, self.val, other.val, self.ctx.val)
+
+    def _iadd_mpoly_(self, other: fmpz_mod_mpoly):
+        fmpz_mod_mpoly_add(self.val, self.val, other.val, self.ctx.val)
+
+    def _isub_scalar_(self, other: fmpz):
+        fmpz_mod_mpoly_sub_fmpz(self.val, self.val, other.val, self.ctx.val)
+
+    def _isub_mpoly_(self, other: fmpz_mod_mpoly):
+        fmpz_mod_mpoly_sub(self.val, self.val, other.val, self.ctx.val)
+
+    def _imul_scalar_(self, other: fmpz):
+        fmpz_mod_mpoly_scalar_mul_fmpz(self.val, self.val, other.val, self.ctx.val)
+
+    def _imul_mpoly_(self, other: fmpz_mod_mpoly):
+        fmpz_mod_mpoly_mul(self.val, self.val, other.val, self.ctx.val)
+
+    def monoms(self):
+        """
+        Return the exponent vectors of each term as a tuple of fmpz.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f.monoms()
+            [(1, 1), (1, 0), (0, 1), (0, 0)]
+
+        """
+        cdef:
+            slong i, nvars = self.ctx.nvars()
+            fmpz_vec vec = fmpz_vec(nvars, double_indirect=True)
+
+        res = []
+        for i in range(len(self)):
+            fmpz_mod_mpoly_get_term_exp_fmpz(vec.double_indirect, self.val, i, self.ctx.val)
+            res.append(vec.to_tuple())
+
+        return res
+
+    def coeffs(self):
+        """
+        Return the coefficients of each term as a fmpz
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f.coeffs()
+            [4, 2, 3, 1]
+
+        """
+        cdef:
+            fmpz coeff
+            slong i
+
+        res = []
+        for i in range(len(self)):
+            coeff = fmpz.__new__(fmpz)
+            fmpz_mod_mpoly_get_term_coeff_fmpz(coeff.val, self.val, i, self.ctx.val)
+            res.append(coeff)
+
+        return res
+
+    # def terms(self):
+    #     """
+    #     Return the terms of this polynomial as a list of fmpz_mod_mpolys.
+
+    #         >>> from flint import Ordering
+    #         >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+    #         >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+    #         >>> f.terms()
+    #         [4*x0*x1, 2*x0, 3*x1, 1]
+
+    #     """
+    #     cdef:
+    #         fmpz_mod_mpoly term
+    #         slong i
+
+    #     res = []
+    #     for i in range(len(self)):
+    #         term = create_fmpz_mod_mpoly(self.ctx)
+    #         fmpz_mod_mpoly_get_term(term.val, self.val, i, self.ctx.val)
+    #         res.append(term)
+
+    #     return res
+
+    def subs(self, dict_args) -> fmpz_mod_mpoly:
+        """
+        Partial evaluate this polynomial with select constants. Keys must be generator names or generator indices,
+        all values must be fmpz.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f.subs({"x1": 0})
+            2*x0 + 1
+
+        """
+        cdef:
+            fmpz_mod_mpoly res
+            slong i
+
+        args = tuple((self.ctx.variable_to_index(k), self.ctx.any_as_scalar(v)) for k, v in dict_args.items())
+        for (_, v), old in zip(args, dict_args.values()):
+            if v is NotImplemented:
+                raise TypeError(f"cannot coerce {type(old)} to fmpz")
+
+        # Partial application with args in Z. We evaluate the polynomial one variable at a time
+        res = create_fmpz_mod_mpoly(self.ctx)
+
+        fmpz_mod_mpoly_set(res.val, self.val, self.ctx.val)
+        for i, arg in args:
+            fmpz_mod_mpoly_evaluate_one_fmpz(res.val, res.val, i, (<fmpz>arg).val, self.ctx.val)
+        return res
+
+    def compose(self, *args, ctx=None) -> fmpz_mod_mpoly:
+        """
+        Compose this polynomial with other fmpz_mod_mpolys. All arguments must share the same context, it may different
+        from this polynomials context.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(1, Ordering.lex, 11, 'x')
+            >>> ctx1 = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'y')
+            >>> f = ctx.from_dict({(2,): 1})
+            >>> g = ctx1.from_dict({(1, 0): 1, (0, 1): 1})
+            >>> f
+            x^2
+            >>> g
+            y0 + y1
+            >>> f.compose(g)
+            y0^2 + 2*y0*y1 + y1^2
+
+        """
+        cdef:
+            fmpz_mod_mpoly res
+            fmpz_mod_mpoly_ctx res_ctx
+            fmpz_mod_mpoly_vec C
+            slong i, nvars = self.ctx.nvars(), nargs = len(args)
+
+        if nargs != nvars:
+            raise ValueError("number of generators does not match number of arguments")
+        elif self.ctx.nvars() == 0 and ctx is None:
+            raise ValueError("a context must be provided when composing a polynomial with no generators")
+        elif not all(typecheck(arg, fmpz_mod_mpoly) for arg in args):
+            raise TypeError("all arguments must be fmpz_mod_mpolys")
+
+        if ctx is None:
+            res_ctx = (<fmpz_mod_mpoly> args[0]).ctx
+        elif typecheck(ctx, fmpz_mod_mpoly_ctx):
+            res_ctx = <fmpz_mod_mpoly_ctx>ctx
+        else:
+            raise TypeError(f"provided context ({ctx}) is not a fmpz_mod_mpoly_ctx")
+
+        if not all((<fmpz_mod_mpoly> arg).ctx is res_ctx for arg in args):
+            raise IncompatibleContextError(
+                "all arguments must share the " + ("same" if ctx is not None else "provided") + " context"
+            )
+
+        C = fmpz_mod_mpoly_vec(args, res_ctx, double_indirect=True)
+        res = create_fmpz_mod_mpoly(res_ctx)
+        if fmpz_mod_mpoly_compose_fmpz_mod_mpoly(res.val, self.val, C.double_indirect, self.ctx.val, res_ctx.val) == 0:
+            raise ValueError("unreasonably large polynomial")  # pragma: no cover
+        return res
+
+    def context(self):
+        """
+        Return the context object for this polynomials.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2})
+            >>> ctx is p.context()
+            True
+        """
+        return self.ctx
+
+    def coefficient(self, slong i):
+        """
+        Return the coefficient at index `i`.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p.coefficient(1)
+            2
+        """
+        cdef fmpz v
+        if not 0 <= i < fmpz_mod_mpoly_length(self.val, self.ctx.val):
+            raise IndexError("term index out of range")
+        else:
+            v = fmpz.__new__(fmpz)
+            fmpz_mod_mpoly_get_term_coeff_fmpz(v.val, self.val, i, self.ctx.val)
+            return v
+
+    def monomial(self, slong i):
+        """
+        Return the exponent vector at index `i` as a tuple.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p.monomial(1)
+            (0, 1)
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        if not 0 <= i < fmpz_mod_mpoly_length(self.val, self.ctx.val):
+            raise IndexError("term index out of range")
+        res = fmpz_vec(nvars, double_indirect=True)
+        fmpz_mod_mpoly_get_term_exp_fmpz(res.double_indirect, self.val, i, self.ctx.val)
+        return res.to_tuple()
+
+    def degrees(self):
+        """
+        Return a dictionary of variable name to degree.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(1, 0, 0, 0): 1, (0, 2, 0, 0): 2, (0, 0, 3, 0): 3})
+            >>> p.degrees()
+            (1, 2, 3, 0)
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        res = fmpz_vec(nvars, double_indirect=True)
+        fmpz_mod_mpoly_degrees_fmpz(res.double_indirect, self.val, self.ctx.val)
+        return res.to_tuple()
+
+    def total_degree(self):
+        """
+        Return the total degree.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(4, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(1, 0, 0, 0): 1, (0, 2, 0, 0): 2, (0, 0, 3, 0): 3})
+            >>> p.total_degree()
+            3
+        """
+        cdef fmpz res = fmpz()
+        fmpz_mod_mpoly_total_degree_fmpz((<fmpz> res).val, self.val, self.ctx.val)
+        return res
+
+    def repr(self):
+        return f"{self.ctx}.from_dict({self.to_dict()})"
+
+    def str(self):
+        cdef bytes s = <bytes> fmpz_mod_mpoly_get_str_pretty(self.val, self.ctx.c_names, self.ctx.val)
+        res = s.decode().replace("+", " + ").replace("-", " - ")
+        if res.startswith(" - "):
+            res = "-" + res[3:]
+        return res
+
+    def gcd(self, other):
+        """
+        Return the gcd of self and other.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(1, 1): 4, (0, 0): 1})
+            >>> g = ctx.from_dict({(0, 1): 2, (1, 0): 2})
+            >>> (f * g).gcd(f)
+            4*x0*x1 + 1
+        """
+        cdef fmpz_mod_mpoly res
+        if not self.ctx.is_prime():
+            raise DomainError("gcd with non-prime modulus is not supported")
+        elif not typecheck(other, fmpz_mod_mpoly):
+            raise TypeError("argument must be a fmpz_mod_mpoly")
+        elif (<fmpz_mod_mpoly>self).ctx is not (<fmpz_mod_mpoly>other).ctx:
+            raise IncompatibleContextError(f"{(<fmpz_mod_mpoly>self).ctx} is not {(<fmpz_mod_mpoly>other).ctx}")
+        res = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_gcd(res.val, (<fmpz_mod_mpoly>self).val, (<fmpz_mod_mpoly>other).val, res.ctx.val)
+        return res
+
+    def sqrt(self):
+        """
+        Return the square root of self.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(1, 1): 4, (0, 0): 1})
+            >>> (f * f).sqrt()
+            4*x0*x1 + 1
+        """
+        cdef fmpz_mod_mpoly res
+        if not self.ctx.is_prime():
+            raise DomainError("square root with non-prime modulus is not supported")
+
+        res = create_fmpz_mod_mpoly(self.ctx)
+
+        if fmpz_mod_mpoly_sqrt(res.val, self.val, self.ctx.val):
+            return res
+        else:
+            raise ValueError("polynomial is not a perfect square")
+
+    def factor(self):
+        """
+        Factors self into irreducible factors, returning a tuple
+        (c, factors) where c is the content of the coefficients and
+        factors is a list of (poly, exp) pairs.
+
+            >>> from flint import Ordering
+            >>> Zm = fmpz_mod_mpoly
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(3, Ordering.lex, 11, 'x,y,z')
+            >>> p1 = Zm("2*x + 4", ctx)
+            >>> p2 = Zm("3*x*z + 3*x + 3*z + 3", ctx)
+            >>> (p1 * p2).factor()
+            (6, [(z + 1, 1), (x + 2, 1), (x + 1, 1)])
+            >>> (p2 * p1 * p2).factor()
+            (18, [(z + 1, 2), (x + 2, 1), (x + 1, 2)])
+        """
+        cdef:
+            fmpz_mod_mpoly_factor_t fac
+            fmpz c
+            fmpz_mod_mpoly u
+
+        if not self.ctx.is_prime():
+            raise DomainError("factorisation with non-prime modulus is not supported")
+
+        fmpz_mod_mpoly_factor_init(fac, self.ctx.val)
+        if not fmpz_mod_mpoly_factor(fac, self.val, self.ctx.val):
+            raise RuntimeError("factorisation failed")
+        res = [0] * fac.num
+
+        for i in range(fac.num):
+            u = create_fmpz_mod_mpoly(self.ctx)
+            fmpz_mod_mpoly_set((<fmpz_mod_mpoly>u).val, &fac.poly[i], self.ctx.val)
+
+            c = fmpz.__new__(fmpz)
+            fmpz_set((<fmpz>c).val, &fac.exp[i])
+
+            res[i] = (u, c)
+
+        c = fmpz.__new__(fmpz)
+        fmpz_set((<fmpz>c).val, fac.constant)
+        fmpz_mod_mpoly_factor_clear(fac, self.ctx.val)
+        return c, res
+
+    def factor_squarefree(self):
+        """
+        Factors self into irreducible factors, returning a tuple
+        (c, factors) where c is the content of the coefficients and
+        factors is a list of (poly, exp) pairs.
+
+            >>> from flint import Ordering
+            >>> Zm = fmpz_mod_mpoly
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(3, Ordering.lex, 11, 'x,y,z')
+            >>> p1 = Zm("2*x + 4", ctx)
+            >>> p2 = Zm("3*x*y + 3*x + 3*y + 3", ctx)
+            >>> (p1 * p2).factor_squarefree()
+            (6, [(y + 1, 1), (x^2 + 3*x + 2, 1)])
+            >>> (p1 * p2 * p1).factor_squarefree()
+            (12, [(y + 1, 1), (x + 1, 1), (x + 2, 2)])
+        """
+        cdef:
+            fmpz_mod_mpoly_factor_t fac
+            fmpz c
+            fmpz_mod_mpoly u
+
+        if not self.ctx.is_prime():
+            raise DomainError("factorisation with non-prime modulus is not supported")
+
+        fmpz_mod_mpoly_factor_init(fac, self.ctx.val)
+        if not fmpz_mod_mpoly_factor_squarefree(fac, self.val, self.ctx.val):
+            raise RuntimeError("factorisation failed")
+        res = [0] * fac.num
+
+        for i in range(fac.num):
+            u = create_fmpz_mod_mpoly(self.ctx)
+            fmpz_mod_mpoly_init(u.val, u.ctx.val)
+            fmpz_mod_mpoly_set((<fmpz_mod_mpoly>u).val, &fac.poly[i], self.ctx.val)
+
+            c = fmpz.__new__(fmpz)
+            fmpz_set((<fmpz>c).val, &fac.exp[i])
+
+            res[i] = (u, c)
+
+        c = fmpz.__new__(fmpz)
+        fmpz_set((<fmpz>c).val, fac.constant)
+        fmpz_mod_mpoly_factor_clear(fac, self.ctx.val)
+        return c, res
+
+    # TODO: Rethink context conversions, particularly the proposed methods in #132
+    # def coerce_to_context(self, ctx):
+    #     cdef:
+    #         fmpz_mod_mpoly res
+    #         slong *C
+    #         slong i
+
+    #     if not typecheck(ctx, fmpz_mod_mpoly_ctx):
+    #         raise ValueError("provided context is not a fmpz_mod_mpoly_ctx")
+
+    #     if self.ctx is ctx:
+    #         return self
+
+    #     C = <slong *> libc.stdlib.malloc(self.ctx.val.minfo.nvars * sizeof(slong *))
+    #     if C is NULL:
+    #         raise MemoryError("malloc returned a null pointer")
+    #     res = create_fmpz_mod_mpoly(self.ctx)
+
+    #     vars = {x: i for i, x in enumerate(ctx.py_names)}
+    #     for i, var in enumerate(self.ctx.py_names):
+    #         C[i] = <slong>vars[var]
+
+    #     fmpz_mod_mpoly_compose_fmpz_mod_mpoly_gen(res.val, self.val, C, self.ctx.val, (<fmpz_mod_mpoly_ctx>ctx).val)
+
+    #     libc.stdlib.free(C)
+    #     return res
+
+    def derivative(self, var):
+        """
+        Return the derivative of this polynomial with respect to the provided variable.
+        The argument can either be the variable as a string, or the index of the
+        variable in the context.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpz_mod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 3): 2, (2, 1): 3})
+            >>> p
+            3*x0^2*x1 + 2*x1^3
+            >>> p.derivative("x0")
+            6*x0*x1
+            >>> p.derivative(1)
+            3*x0^2 + 6*x1^2
+
+        """
+        cdef:
+            fmpz_mod_mpoly res
+            slong i = self.ctx.variable_to_index(var)
+
+        res = create_fmpz_mod_mpoly(self.ctx)
+
+        fmpz_mod_mpoly_derivative(res.val, self.val, i, self.ctx.val)
+        return res
+
+
+cdef class fmpz_mod_mpoly_vec:
+    """
+    A class representing a vector of fmpz_mod_mpolys.
+    """
+
+    def __cinit__(self, iterable_or_len, fmpz_mod_mpoly_ctx ctx, bint double_indirect = False):
+        if isinstance(iterable_or_len, int):
+            self.length = iterable_or_len
+        else:
+            self.length = len(iterable_or_len)
+
+        self.ctx = ctx
+
+        self.val = <fmpz_mod_mpoly_struct *> libc.stdlib.malloc(self.length * sizeof(fmpz_mod_mpoly_struct))
+        for i in range(self.length):
+            fmpz_mod_mpoly_init(&self.val[i], self.ctx.val)
+
+        if double_indirect:
+            self.double_indirect = <fmpz_mod_mpoly_struct **> libc.stdlib.malloc(self.length * sizeof(fmpz_mod_mpoly_struct *))
+            if self.double_indirect is NULL:
+                raise MemoryError("malloc returned a null pointer")  # pragma: no cover
+
+            for i in range(self.length):
+                self.double_indirect[i] = &self.val[i]
+        else:
+            self.double_indirect = NULL
+
+    def __init__(self, iterable_or_len, _, double_indirect: bool = False):
+        if not isinstance(iterable_or_len, int):
+            for i, x in enumerate(iterable_or_len):
+                self[i] = x
+
+    def __dealloc__(self):
+        libc.stdlib.free(self.double_indirect)
+        for i in range(self.length):
+            fmpz_mod_mpoly_clear(&self.val[i], self.ctx.val)
+        libc.stdlib.free(self.val)
+
+    def __getitem__(self, x):
+        if not isinstance(x, int):
+            raise TypeError("index is not integer")
+        elif not 0 <= x < self.length:
+            raise IndexError("index out of range")
+
+        cdef fmpz_mod_mpoly z = create_fmpz_mod_mpoly(self.ctx)
+        fmpz_mod_mpoly_set(z.val, &self.val[x], self.ctx.val)
+        return z
+
+    def __setitem__(self, x, y):
+        if not isinstance(x, int):
+            raise TypeError("index is not integer")
+        elif not 0 <= x < self.length:
+            raise IndexError("index out of range")
+        elif not typecheck(y, fmpz_mod_mpoly):
+            raise TypeError("argument is not fmpz_mod_mpoly")
+        elif (<fmpz_mod_mpoly>y).ctx is not self.ctx:
+            raise IncompatibleContextError(f"{(<fmpz_mod_mpoly>y).ctx} is not {self.ctx}")
+
+        fmpz_mod_mpoly_set(&self.val[x], (<fmpz_mod_mpoly>y).val, self.ctx.val)
+
+    def __len__(self):
+        return self.val.length
+
+    def __str__(self):
+        s = [None] * self.length
+        for i in range(self.length):
+            x = create_fmpz_mod_mpoly(self.ctx)
+            fmpz_mod_mpoly_set(x.val, &self.val[x], self.ctx.val)
+            s[i] = str(x)
+        return f"[{', '.join(s)}]"
+
+    def __repr__(self):
+        return f"fmpz_mod_mpoly_vec({self}, ctx={self.ctx})"
+
+    def to_tuple(self):
+        return tuple(self[i] for i in range(self.val.length))
diff --git a/src/flint/types/fmpz_mpoly.pyx b/src/flint/types/fmpz_mpoly.pyx
index d8e2d869..b964050d 100644
--- a/src/flint/types/fmpz_mpoly.pyx
+++ b/src/flint/types/fmpz_mpoly.pyx
@@ -24,13 +24,13 @@ from flint.flintlib.fmpz_mpoly cimport (
     fmpz_mpoly_divides,
     fmpz_mpoly_divrem,
     fmpz_mpoly_equal,
+    fmpz_mpoly_equal_fmpz,
     fmpz_mpoly_evaluate_all_fmpz,
     fmpz_mpoly_evaluate_one_fmpz,
     fmpz_mpoly_gcd,
     fmpz_mpoly_gen,
     fmpz_mpoly_get_coeff_fmpz_fmpz,
     fmpz_mpoly_get_str_pretty,
-    fmpz_mpoly_get_term,
     fmpz_mpoly_get_term_coeff_fmpz,
     fmpz_mpoly_get_term_exp_fmpz,
     fmpz_mpoly_integral,
@@ -234,13 +234,18 @@ cdef class fmpz_mpoly(flint_mpoly):
             return NotImplemented
         elif other is None:
             return op == Py_NE
-        elif typecheck(self, fmpz_mpoly) and typecheck(other, fmpz_mpoly):
+        elif typecheck(other, fmpz_mpoly):
             if (<fmpz_mpoly>self).ctx is (<fmpz_mpoly>other).ctx:
-                return (op == Py_NE) ^ bool(
-                    fmpz_mpoly_equal((<fmpz_mpoly>self).val, (<fmpz_mpoly>other).val, (<fmpz_mpoly>self).ctx.val)
-                )
+                return (op == Py_NE) ^ <bint>fmpz_mpoly_equal(self.val, (<fmpz_mpoly>other).val, self.ctx.val)
             else:
                 return op == Py_NE
+        elif typecheck(other, fmpz):
+            return (op == Py_NE) ^ <bint>fmpz_mpoly_equal_fmpz(self.val, (<fmpz>other).val, self.ctx.val)
+        elif isinstance(other, int):
+            other = any_as_fmpz(other)
+            if other is NotImplemented:
+                return NotImplemented
+            return (op == Py_NE) ^ <bint>fmpz_mpoly_equal_fmpz(self.val, (<fmpz>other).val, self.ctx.val)
         else:
             return NotImplemented
 
@@ -515,17 +520,10 @@ cdef class fmpz_mpoly(flint_mpoly):
             fmpz vres
             slong nvars = self.ctx.nvars(), nargs = len(args)
 
-        if nargs < nvars:
-            raise ValueError("not enough arguments provided")
-        elif nargs > nvars:
-            raise ValueError("more arguments provided than variables")
-
-        args_fmpz = tuple(any_as_fmpz(v) for v in args)
-        for arg in args_fmpz:
-            if arg is NotImplemented:
-                raise TypeError(f"cannot coerce argument ('{arg}') to fmpz")
+        if nargs != nvars:
+            raise ValueError("number of generators does not match number of arguments")
 
-        V = fmpz_vec(args_fmpz, double_indirect=True)
+        V = fmpz_vec(args, double_indirect=True)
         vres = fmpz.__new__(fmpz)
         if fmpz_mpoly_evaluate_all_fmpz(vres.val, self.val, V.double_indirect, self.ctx.val) == 0:
             raise ValueError("unreasonably large polynomial")  # pragma: no cover
@@ -898,7 +896,7 @@ cdef class fmpz_mpoly(flint_mpoly):
             >>> Zm = fmpz_mpoly
             >>> ctx = fmpz_mpoly_ctx.get_context(3, Ordering.lex, 'x,y,z')
             >>> p1 = Zm("2*x + 4", ctx)
-            >>> p2 = Zm("3*x*z +  + 3*x + 3*z + 3", ctx)
+            >>> p2 = Zm("3*x*z + 3*x + 3*z + 3", ctx)
             >>> (p1 * p2).factor()
             (6, [(z + 1, 1), (x + 2, 1), (x + 1, 1)])
             >>> (p2 * p1 * p2).factor()
@@ -910,7 +908,8 @@ cdef class fmpz_mpoly(flint_mpoly):
             fmpz_mpoly u
 
         fmpz_mpoly_factor_init(fac, self.ctx.val)
-        fmpz_mpoly_factor(fac, self.val, self.ctx.val)
+        if not fmpz_mpoly_factor(fac, self.val, self.ctx.val):
+            raise RuntimeError("factorisation failed")
         res = [0] * fac.num
 
         for i in range(fac.num):
diff --git a/src/flint/types/fmpz_vec.pyx b/src/flint/types/fmpz_vec.pyx
index d5d7e27e..365a1d38 100644
--- a/src/flint/types/fmpz_vec.pyx
+++ b/src/flint/types/fmpz_vec.pyx
@@ -3,6 +3,9 @@ from flint.flintlib.flint cimport slong
 from flint.flintlib.fmpz_vec cimport _fmpz_vec_init, _fmpz_vec_clear
 
 from flint.types.fmpz cimport fmpz, any_as_fmpz
+from flint.types.fmpz_mod cimport fmpz_mod
+
+from flint.utils.typecheck cimport typecheck
 
 cimport libc.stdlib
 
diff --git a/src/flint/types/meson.build b/src/flint/types/meson.build
index faf7f31e..6f32aff5 100644
--- a/src/flint/types/meson.build
+++ b/src/flint/types/meson.build
@@ -19,6 +19,7 @@ exts = [
 
   'nmod',
   'nmod_poly',
+  'nmod_mpoly',
   'nmod_mat',
   'nmod_series',
 
@@ -44,6 +45,7 @@ exts = [
   'dirichlet',
 
   'fmpz_mpoly',
+  'fmpz_mod_mpoly',
   'fmpq_mpoly',
 ]
 
diff --git a/src/flint/types/nmod_mpoly.pxd b/src/flint/types/nmod_mpoly.pxd
new file mode 100644
index 00000000..a7249f32
--- /dev/null
+++ b/src/flint/types/nmod_mpoly.pxd
@@ -0,0 +1,37 @@
+from flint.flint_base.flint_base cimport flint_mpoly, flint_mpoly_context
+
+from flint.flintlib.nmod_mpoly cimport (
+    nmod_mpoly_ctx_t,
+    nmod_mpoly_t,
+    nmod_mpoly_init,
+    nmod_mpoly_struct
+)
+from flint.flintlib.flint cimport slong, ulong
+
+cdef inline init_nmod_mpoly(nmod_mpoly var, nmod_mpoly_ctx ctx):
+    var.ctx = ctx
+    nmod_mpoly_init(var.val, ctx.val)
+    var._init = True
+
+cdef inline nmod_mpoly create_nmod_mpoly(nmod_mpoly_ctx ctx):
+    cdef nmod_mpoly var
+    var = nmod_mpoly.__new__(nmod_mpoly)
+    var.ctx = ctx
+    nmod_mpoly_init(var.val, ctx.val)
+    var._init = True
+    return var
+
+cdef class nmod_mpoly_ctx(flint_mpoly_context):
+    cdef nmod_mpoly_ctx_t val
+    cdef readonly object __prime_modulus
+
+cdef class nmod_mpoly(flint_mpoly):
+    cdef nmod_mpoly_t val
+    cdef nmod_mpoly_ctx ctx
+    cdef bint _init
+
+cdef class nmod_mpoly_vec:
+    cdef nmod_mpoly_struct *val
+    cdef slong length
+    cdef nmod_mpoly_ctx ctx
+    cdef nmod_mpoly_struct **double_indirect
diff --git a/src/flint/types/nmod_mpoly.pyx b/src/flint/types/nmod_mpoly.pyx
new file mode 100644
index 00000000..22fcacc4
--- /dev/null
+++ b/src/flint/types/nmod_mpoly.pyx
@@ -0,0 +1,1003 @@
+from flint.flint_base.flint_base cimport (
+    flint_mpoly,
+    flint_mpoly_context,
+    ordering_py_to_c,
+    ordering_c_to_py,
+)
+from flint.flint_base.flint_base import Ordering
+
+from flint.utils.typecheck cimport typecheck
+from flint.utils.flint_exceptions import DomainError, IncompatibleContextError
+
+from flint.types.fmpz cimport any_as_fmpz, fmpz
+from flint.types.fmpz_vec cimport fmpz_vec
+from flint.types.fmpz_mod cimport fmpz_mod
+
+from flint.types.nmod cimport nmod
+
+from flint.flintlib.flint cimport SIZEOF_ULONG
+from flint.flintlib.fmpz cimport fmpz_set
+from flint.flintlib.nmod cimport fmpz_get_nmod
+from flint.flintlib.nmod_mpoly cimport (
+    nmod_mpoly_add,
+    nmod_mpoly_add_ui,
+    nmod_mpoly_clear,
+    nmod_mpoly_compose_nmod_mpoly,
+    nmod_mpoly_ctx_modulus,
+    nmod_mpoly_ctx_init,
+    nmod_mpoly_degrees_fmpz,
+    nmod_mpoly_derivative,
+    nmod_mpoly_div,
+    nmod_mpoly_divides,
+    nmod_mpoly_divrem,
+    nmod_mpoly_equal,
+    nmod_mpoly_equal_ui,
+    nmod_mpoly_evaluate_all_ui,
+    nmod_mpoly_evaluate_one_ui,
+    nmod_mpoly_gcd,
+    nmod_mpoly_gen,
+    nmod_mpoly_get_coeff_ui_fmpz,
+    nmod_mpoly_get_str_pretty,
+    nmod_mpoly_get_term,
+    nmod_mpoly_get_term_coeff_ui,
+    nmod_mpoly_get_term_exp_fmpz,
+    nmod_mpoly_is_one,
+    nmod_mpoly_is_zero,
+    nmod_mpoly_length,
+    nmod_mpoly_mul,
+    nmod_mpoly_neg,
+    nmod_mpoly_pow_fmpz,
+    nmod_mpoly_push_term_ui_ffmpz,
+    nmod_mpoly_scalar_mul_ui,
+    nmod_mpoly_set,
+    nmod_mpoly_set_coeff_ui_fmpz,
+    nmod_mpoly_set_ui,
+    nmod_mpoly_set_str_pretty,
+    nmod_mpoly_sort_terms,
+    nmod_mpoly_sub,
+    nmod_mpoly_sub_ui,
+    nmod_mpoly_total_degree_fmpz,
+    nmod_mpoly_sqrt,
+)
+from flint.flintlib.nmod_mpoly_factor cimport (
+    nmod_mpoly_factor,
+    nmod_mpoly_factor_clear,
+    nmod_mpoly_factor_init,
+    nmod_mpoly_factor_squarefree,
+    nmod_mpoly_factor_t,
+)
+from flint.flintlib.ulong_extras cimport n_is_prime
+
+from cpython.object cimport Py_EQ, Py_NE
+cimport libc.stdlib
+
+cdef dict _nmod_mpoly_ctx_cache = {}
+
+
+cdef class nmod_mpoly_ctx(flint_mpoly_context):
+    """
+    A class for storing the polynomial context
+
+    :param nvars: The number of variables in the ring
+    :param ordering:  The term order for the ring
+    :param names:  A tuple containing the names of the variables of the ring.
+
+    Do not construct one of these directly, use `nmod_mpoly_ctx.get_context`.
+    """
+
+    _ctx_cache = _nmod_mpoly_ctx_cache
+
+    def __init__(self, slong nvars, ordering, names, modulus: int):
+        if modulus <= 0:
+            raise ValueError("modulus must be positive")
+
+        nmod_mpoly_ctx_init(self.val, nvars, ordering_py_to_c(ordering), modulus)
+        self.__prime_modulus = None
+        super().__init__(nvars, names)
+
+    @classmethod
+    def create_context_key(
+            cls,
+            slong nvars=1,
+            ordering=Ordering.lex,
+            modulus = None,
+            names: Optional[str] = "x",
+            nametup: Optional[tuple] = None,
+    ):
+        """
+        Create a key for the context cache via the number of variables, the ordering, the modulus, and either a
+        variable name string, or a tuple of variable names.
+        """
+        # A type hint of `ordering: Ordering` results in the error "TypeError: an integer is required" if a Ordering
+        # object is not provided. This is pretty obtuse so we check its type ourselves
+        if not isinstance(ordering, Ordering):
+            raise TypeError(f"`ordering` ('{ordering}') is not an instance of flint.Ordering")
+
+        if nametup is not None:
+            key = nvars, ordering, nametup, modulus
+        elif nametup is None and names is not None:
+            key = nvars, ordering, cls.create_variable_names(nvars, names), modulus
+        else:
+            raise ValueError("must provide either `names` or `nametup`")
+        return key
+
+    def any_as_scalar(self, other):
+        if isinstance(other, int):
+            try:
+                return <ulong>other
+            except OverflowError:
+                return <ulong>(other % self.modulus())
+        elif typecheck(other, nmod):
+            if (<nmod>other).modulus() != self.modulus():
+                raise DomainError(
+                    f"modulus does not match, got {(<nmod>other).modulus()}, expected {self.modulus()}"
+                )
+            return (<nmod>other).val
+        elif typecheck(other, fmpz):
+            return fmpz_get_nmod((<fmpz>other).val, self.val.mod)
+        elif typecheck(other, fmpz_mod):
+            if (<fmpz_mod>other).ctx.modulus() != self.modulus():
+                raise DomainError(
+                    f"modulus does not match, got {(<fmpz_mod>other).ctx.modulus()}, expected {self.modulus()}"
+                )
+            return fmpz_get_nmod((<fmpz_mod>other).val, self.val.mod)
+        else:
+            return NotImplemented
+
+    def scalar_as_mpoly(self, other: ulong):
+        # non-ulong scalars should first be converted via cls.any_as_scalar
+        return self.constant(<ulong>other)
+
+    def nvars(self):
+        """
+        Return the number of variables in the context
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.lex, 11, 'x')
+            >>> ctx.nvars()
+            4
+        """
+        return self.val.minfo.nvars
+
+    def ordering(self):
+        """
+        Return the term order of the context object.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.deglex, 11, 'w')
+            >>> ctx.ordering()
+            <Ordering.deglex: 1>
+        """
+        return ordering_c_to_py(self.val.minfo.ord)
+
+    def modulus(self):
+        """
+        Return the modulus of the context object.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.deglex, 2, 'w')
+            >>> ctx.modulus()
+            2
+
+        """
+        return nmod_mpoly_ctx_modulus(self.val)
+
+    def is_prime(self):
+        """
+        Return whether the modulus is prime
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.degrevlex, 2**127, 'z')
+            >>> ctx.is_prime()
+            False
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.degrevlex, 2**127 - 1, 'z')
+            >>> ctx.is_prime()
+            True
+        """
+        if self.__prime_modulus is None:
+            self.__prime_modulus = n_is_prime(self.modulus())
+        return self.__prime_modulus
+
+    def gen(self, slong i):
+        """
+        Return the `i`th generator of the polynomial ring
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(3, Ordering.degrevlex, 11, 'z')
+            >>> ctx.gen(1)
+            z1
+        """
+        cdef nmod_mpoly res
+        if not 0 <= i < self.val.minfo.nvars:
+            raise IndexError("generator index out of range")
+        res = create_nmod_mpoly(self)
+        nmod_mpoly_gen(res.val, i, res.ctx.val)
+        return res
+
+    def constant(self, z):
+        """
+        Create a constant polynomial in this context
+        """
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self)
+        nmod_mpoly_set_ui(res.val, z, res.ctx.val)
+        return res
+
+    def from_dict(self, d):
+        """
+        Create a nmod_mpoly from a dictionary in this context.
+
+        The dictionary's keys are tuples of ints (or anything that implicitly converts
+        to fmpz) representing exponents, and corresponding values of fmpz.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x,y')
+            >>> ctx.from_dict({(1,0):2, (1,1):3, (0,1):1})
+            3*x*y + 2*x + y
+        """
+        cdef:
+            fmpz_vec exp_vec
+            slong i, nvars = self.nvars()
+            nmod_mpoly res
+
+        if not isinstance(d, dict):
+            raise ValueError("expected a dictionary")
+
+        res = create_nmod_mpoly(self)
+
+        for i, (exps, coeff) in enumerate(d.items()):
+            if len(exps) != nvars:
+                raise ValueError(f"expected {nvars} exponents, got {len(exps)}")
+            elif not coeff:
+                continue
+
+            exp_vec = fmpz_vec(exps)
+            coeff_scalar = self.any_as_scalar(coeff)
+            if coeff_scalar is NotImplemented:
+                raise TypeError(f"cannot coerce {repr(coeff)} to nmod_mpoly coefficient")
+
+            nmod_mpoly_push_term_ui_ffmpz(res.val, coeff_scalar, exp_vec.val, self.val)
+
+        nmod_mpoly_sort_terms(res.val, self.val)
+        return res
+
+
+cdef class nmod_mpoly(flint_mpoly):
+    """
+    The *nmod_mpoly* type represents sparse multivariate polynomials over
+    the integers.
+    """
+
+    def __cinit__(self):
+        self._init = False
+
+    def __dealloc__(self):
+        if self._init:
+            nmod_mpoly_clear(self.val, self.ctx.val)
+            self._init = False
+
+    def __init__(self, val=0, ctx=None):
+        if typecheck(val, nmod_mpoly):
+            if ctx is None or ctx == (<nmod_mpoly>val).ctx:
+                init_nmod_mpoly(self, (<nmod_mpoly>val).ctx)
+                nmod_mpoly_set(self.val, (<nmod_mpoly>val).val, self.ctx.val)
+            else:
+                raise IncompatibleContextError(f"{ctx} is not {(<nmod_mpoly>val).ctx}")
+        elif isinstance(val, dict):
+            if ctx is None:
+                raise ValueError("a context is required to create a nmod_mpoly from a dict")
+            x = ctx.from_dict(val)
+            # XXX: this copy is silly, have a ctx function that assigns an nmod_mpoly_t
+            init_nmod_mpoly(self, ctx)
+            nmod_mpoly_set(self.val, (<nmod_mpoly>x).val, self.ctx.val)
+        elif isinstance(val, str):
+            if ctx is None:
+                raise ValueError("cannot parse a polynomial without context")
+            val = bytes(val, 'utf-8')
+            init_nmod_mpoly(self, ctx)
+            if nmod_mpoly_set_str_pretty(self.val, val, self.ctx.c_names, self.ctx.val) == -1:
+                raise ValueError("unable to parse nmod_mpoly from string")
+            nmod_mpoly_sort_terms(self.val, self.ctx.val)
+        elif isinstance(val, int):
+            if ctx is None:
+                raise ValueError("need context to convert int to nmod_mpoly")
+            init_nmod_mpoly(self, ctx)
+            nmod_mpoly_set_ui(self.val, val, self.ctx.val)
+        else:
+            raise TypeError(f"cannot construct a nmod_mpoly from a {type(val)}")
+
+    def _division_check(self, other):
+        super()._division_check(other)
+        if not self.ctx.is_prime():
+            raise DomainError("division with non-prime modulus is not supported")
+
+    def __bool__(self):
+        return not nmod_mpoly_is_zero(self.val, self.ctx.val)
+
+    def is_zero(self):
+        return <bint>nmod_mpoly_is_zero(self.val, self.ctx.val)
+
+    def is_one(self):
+        return <bint>nmod_mpoly_is_one(self.val, self.ctx.val)
+
+    def __richcmp__(self, other, int op):
+        if not (op == Py_EQ or op == Py_NE):
+            return NotImplemented
+        elif typecheck(self, nmod_mpoly) and typecheck(other, nmod_mpoly):
+            if (<nmod_mpoly>self).ctx is (<nmod_mpoly>other).ctx:
+                return (op == Py_NE) ^ <bint>nmod_mpoly_equal(self.val, (<nmod_mpoly>other).val, self.ctx.val)
+            else:
+                return op == Py_NE
+        elif typecheck(other, nmod):
+            if other.modulus() != self.ctx.modulus():
+                return op == Py_NE
+            return (op == Py_NE) ^ <bint>nmod_mpoly_equal_ui(self.val, int(other), self.ctx.val)
+        elif typecheck(other, fmpz):
+            return (op == Py_NE) ^ <bint>nmod_mpoly_equal_ui(
+                self.val,
+                fmpz_get_nmod((<fmpz>other).val, self.ctx.val.mod),
+                self.ctx.val
+            )
+        elif isinstance(other, int):
+            return (op == Py_NE) ^ <bint>nmod_mpoly_equal_ui(self.val, other, self.ctx.val)
+        else:
+            return NotImplemented
+
+    def __len__(self):
+        return nmod_mpoly_length(self.val, self.ctx.val)
+
+    def __getitem__(self, x):
+        """
+        Return the coefficient of the term with the exponent vector `x`.
+        Always returns a value, missing keys will return `0`.
+        Negative exponents are made positive.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p[1, 1]
+            3
+
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        if not isinstance(x, tuple):
+            raise TypeError("exponent vector index is not a tuple")
+        elif len(x) != nvars:
+            raise ValueError("exponent vector provided does not match number of variables")
+
+        exp_vec = fmpz_vec(x, double_indirect=True)
+        return nmod_mpoly_get_coeff_ui_fmpz(self.val, exp_vec.double_indirect, self.ctx.val)
+
+    def __setitem__(self, x, y):
+        """
+        Set the coefficient of the term with the exponent vector `x` to `y`.
+        Will always set a value, missing keys will create a new term.
+        Negative exponents are made positive.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p[1, 1] = 20
+            >>> p
+            20*x0*x1 + 2*x1
+
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+            ulong coeff
+
+        if not isinstance(x, tuple):
+            raise TypeError("exponent vector index is not a tuple")
+        elif len(x) != nvars:
+            raise ValueError("exponent vector provided does not match number of variables")
+
+        exp_vec = fmpz_vec(x, double_indirect=True)
+
+        coeff = self.ctx.any_as_scalar(y)
+        if coeff is NotImplemented:
+            raise TypeError("provided coefficient not coercible to ulong")
+        nmod_mpoly_set_coeff_ui_fmpz(self.val, coeff, exp_vec.double_indirect, self.ctx.val)
+
+    def __neg__(self):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_neg(res.val, (<nmod_mpoly>self).val, res.ctx.val)
+        return res
+
+    def _add_scalar_(self, other: ulong):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_add_ui(res.val, self.val, other, self.ctx.val)
+        return res
+
+    def _add_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_add(res.val, self.val, other.val, res.ctx.val)
+        return res
+
+    def _sub_scalar_(self, other: ulong):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_sub_ui(res.val, self.val, other, self.ctx.val)
+        return res
+
+    def _sub_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_sub(res.val, self.val, other.val, res.ctx.val)
+        return res
+
+    def _mul_scalar_(self, other: ulong):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_scalar_mul_ui(res.val, self.val, other, self.ctx.val)
+        return res
+
+    def _mul_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_mul(res.val, self.val, other.val, res.ctx.val)
+        return res
+
+    def _pow_(self, other: fmpz):
+        cdef nmod_mpoly res
+        res = create_nmod_mpoly(self.ctx)
+        if nmod_mpoly_pow_fmpz(res.val, self.val, (<fmpz>other).val, res.ctx.val) == 0:
+            raise ValueError("unreasonably large polynomial")  # pragma: no cover
+        return res
+
+    def _divmod_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly quotient, remainder
+        quotient = create_nmod_mpoly(self.ctx)
+        remainder = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_divrem(quotient.val, remainder.val, self.val, other.val, self.ctx.val)
+        return (quotient, remainder)
+
+    def _floordiv_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly quotient
+        quotient = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_div(quotient.val, self.val, other.val, self.ctx.val)
+        return quotient
+
+    def _truediv_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly quotient
+        quotient = create_nmod_mpoly(self.ctx)
+        if nmod_mpoly_divides(quotient.val, self.val, other.val, self.ctx.val):
+            return quotient
+        else:
+            raise DomainError("nmod_mpoly division is not exact")
+
+    def _mod_mpoly_(self, other: nmod_mpoly):
+        cdef nmod_mpoly quotient, remainder
+        quotient = create_nmod_mpoly(self.ctx)
+        remainder = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_divrem(quotient.val, remainder.val, self.val, other.val, self.ctx.val)
+        return remainder
+
+    def __call__(self, *args) -> ulong:
+        cdef:
+            slong nvars = self.ctx.nvars(), nargs = len(args)
+
+        if nargs != nvars:
+            raise ValueError("number of generators does not match number of arguments")
+
+        args = [self.ctx.any_as_scalar(x) for x in args]
+        cdef:
+            # Using sizeof(ulong) here breaks on 64 windows machines because of the `ctypedef unsigned long ulong` in
+            # flintlib/flint.pxd. Cython will inline this definition and then allocate the wrong amount of memory.
+            ulong *vals = <ulong *>libc.stdlib.malloc(nargs * SIZEOF_ULONG)
+            ulong res
+        if vals is NULL:
+            raise MemoryError("malloc returned a null pointer")  # pragma: no cover
+
+        try:
+            for i in range(nargs):
+                vals[i] = args[i]
+            res = nmod_mpoly_evaluate_all_ui(self.val, vals, self.ctx.val)
+        finally:
+            libc.stdlib.free(vals)
+
+        return res
+
+    def _iadd_scalar_(self, other: ulong):
+        nmod_mpoly_add_ui(self.val, self.val, other, self.ctx.val)
+
+    def _iadd_mpoly_(self, other: nmod_mpoly):
+        nmod_mpoly_add(self.val, self.val, other.val, self.ctx.val)
+
+    def _isub_scalar_(self, other: ulong):
+        nmod_mpoly_sub_ui(self.val, self.val, other, self.ctx.val)
+
+    def _isub_mpoly_(self, other: nmod_mpoly):
+        nmod_mpoly_sub(self.val, self.val, other.val, self.ctx.val)
+
+    def _imul_scalar_(self, other: ulong):
+        nmod_mpoly_scalar_mul_ui(self.val, self.val, other, self.ctx.val)
+
+    def _imul_mpoly_(self, other: nmod_mpoly):
+        nmod_mpoly_mul(self.val, self.val, other.val, self.ctx.val)
+
+    def monoms(self):
+        """
+        Return the exponent vectors of each term as a tuple of fmpz.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f.monoms()
+            [(1, 1), (1, 0), (0, 1), (0, 0)]
+
+        """
+        cdef:
+            slong i, nvars = self.ctx.nvars()
+            fmpz_vec vec = fmpz_vec(nvars, double_indirect=True)
+
+        res = []
+        for i in range(len(self)):
+            nmod_mpoly_get_term_exp_fmpz(vec.double_indirect, self.val, i, self.ctx.val)
+            res.append(vec.to_tuple())
+
+        return res
+
+    def coeffs(self):
+        """
+        Return the coefficients of each term as a fmpz
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f.coeffs()
+            [4, 2, 3, 1]
+
+        """
+        return [nmod_mpoly_get_term_coeff_ui(self.val, i, self.ctx.val) for i in range(len(self))]
+
+    # def terms(self):
+    #     """
+    #     Return the terms of this polynomial as a list of nmod_mpolys.
+
+    #         >>> from flint import Ordering
+    #         >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+    #         >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+    #         >>> f.terms()
+    #         [4*x0*x1, 2*x0, 3*x1, 1]
+
+    #     """
+    #     cdef:
+    #         nmod_mpoly term
+    #         slong i
+
+    #     res = []
+    #     for i in range(len(self)):
+    #         term = create_nmod_mpoly(self.ctx)
+    #         nmod_mpoly_get_term(term.val, self.val, i, self.ctx.val)
+    #         res.append(term)
+
+    #     return res
+
+    def subs(self, dict_args) -> nmod_mpoly:
+        """
+        Partial evaluate this polynomial with select constants. Keys must be generator names or generator indices,
+        all values must be fmpz.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(0, 0): 1, (1, 0): 2, (0, 1): 3, (1, 1): 4})
+            >>> f.subs({"x1": 0})
+            2*x0 + 1
+
+        """
+        cdef:
+            nmod_mpoly res
+            slong i
+
+        args = tuple((self.ctx.variable_to_index(k), self.ctx.any_as_scalar(v)) for k, v in dict_args.items())
+        for (_, v), old in zip(args, dict_args.values()):
+            if v is NotImplemented:
+                raise TypeError(f"cannot coerce {type(old)} to ulong")
+
+        # Partial application with args in Z. We evaluate the polynomial one variable at a time
+        res = create_nmod_mpoly(self.ctx)
+
+        nmod_mpoly_set(res.val, self.val, self.ctx.val)
+        for i, arg in args:
+            nmod_mpoly_evaluate_one_ui(res.val, res.val, i, arg, self.ctx.val)
+        return res
+
+    def compose(self, *args, ctx=None) -> nmod_mpoly:
+        """
+        Compose this polynomial with other nmod_mpolys. All arguments must share the same context, it may different
+        from this polynomials context.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(1, Ordering.lex, 11, 'x')
+            >>> ctx1 = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'y')
+            >>> f = ctx.from_dict({(2,): 1})
+            >>> g = ctx1.from_dict({(1, 0): 1, (0, 1): 1})
+            >>> f
+            x^2
+            >>> g
+            y0 + y1
+            >>> f.compose(g)
+            y0^2 + 2*y0*y1 + y1^2
+
+        """
+        cdef:
+            nmod_mpoly res
+            nmod_mpoly_ctx res_ctx
+            nmod_mpoly_vec C
+            slong nvars = self.ctx.nvars(), nargs = len(args)
+
+        if nargs != nvars:
+            raise ValueError("number of generators does not match number of arguments")
+        elif self.ctx.nvars() == 0 and ctx is None:
+            raise ValueError("a context must be provided when composing a polynomial with no generators")
+        elif not all(typecheck(arg, nmod_mpoly) for arg in args):
+            raise TypeError("all arguments must be nmod_mpolys")
+
+        if ctx is None:
+            res_ctx = (<nmod_mpoly> args[0]).ctx
+        elif typecheck(ctx, nmod_mpoly_ctx):
+            res_ctx = <nmod_mpoly_ctx>ctx
+        else:
+            raise TypeError(f"provided context ({ctx}) is not a nmod_mpoly_ctx")
+
+        if not all((<nmod_mpoly> arg).ctx is res_ctx for arg in args):
+            raise IncompatibleContextError(
+                "all arguments must share the " + ("same" if ctx is not None else "provided") + " context"
+            )
+
+        C = nmod_mpoly_vec(args, res_ctx, double_indirect=True)
+        res = create_nmod_mpoly(res_ctx)
+        if nmod_mpoly_compose_nmod_mpoly(res.val, self.val, C.double_indirect, self.ctx.val, res_ctx.val) == 0:
+            raise ValueError("unreasonably large polynomial")  # pragma: no cover
+        return res
+
+    def context(self):
+        """
+        Return the context object for this polynomials.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2})
+            >>> ctx is p.context()
+            True
+        """
+        return self.ctx
+
+    def coefficient(self, slong i):
+        """
+        Return the coefficient at index `i`.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p.coefficient(1)
+            2
+        """
+        if not 0 <= i < nmod_mpoly_length(self.val, self.ctx.val):
+            raise IndexError("term index out of range")
+        else:
+            return nmod_mpoly_get_term_coeff_ui(self.val, i, self.ctx.val)
+
+    def monomial(self, slong i):
+        """
+        Return the exponent vector at index `i` as a tuple.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 1): 2, (1, 1): 3})
+            >>> p.monomial(1)
+            (0, 1)
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        if not 0 <= i < nmod_mpoly_length(self.val, self.ctx.val):
+            raise IndexError("term index out of range")
+        res = fmpz_vec(nvars, double_indirect=True)
+        nmod_mpoly_get_term_exp_fmpz(res.double_indirect, self.val, i, self.ctx.val)
+        return res.to_tuple()
+
+    def degrees(self):
+        """
+        Return a dictionary of variable name to degree.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(1, 0, 0, 0): 1, (0, 2, 0, 0): 2, (0, 0, 3, 0): 3})
+            >>> p.degrees()
+            (1, 2, 3, 0)
+        """
+        cdef:
+            slong nvars = self.ctx.nvars()
+
+        res = fmpz_vec(nvars, double_indirect=True)
+        nmod_mpoly_degrees_fmpz(res.double_indirect, self.val, self.ctx.val)
+        return res.to_tuple()
+
+    def total_degree(self):
+        """
+        Return the total degree.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(4, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(1, 0, 0, 0): 1, (0, 2, 0, 0): 2, (0, 0, 3, 0): 3})
+            >>> p.total_degree()
+            3
+        """
+        cdef fmpz res = fmpz()
+        nmod_mpoly_total_degree_fmpz((<fmpz> res).val, self.val, self.ctx.val)
+        return res
+
+    def repr(self):
+        return f"{self.ctx}.from_dict({self.to_dict()})"
+
+    def str(self):
+        cdef bytes s = <bytes> nmod_mpoly_get_str_pretty(self.val, self.ctx.c_names, self.ctx.val)
+        res = s.decode().replace("+", " + ").replace("-", " - ")
+        if res.startswith(" - "):
+            res = "-" + res[3:]
+        return res
+
+    def gcd(self, other):
+        """
+        Return the gcd of self and other.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(1, 1): 4, (0, 0): 1})
+            >>> g = ctx.from_dict({(0, 1): 2, (1, 0): 2})
+            >>> (f * g).gcd(f)
+            4*x0*x1 + 1
+        """
+        cdef nmod_mpoly res
+        if not self.ctx.is_prime():
+            raise DomainError("gcd with non-prime modulus is not supported")
+        elif not typecheck(other, nmod_mpoly):
+            raise TypeError("argument must be a nmod_mpoly")
+        elif (<nmod_mpoly>self).ctx is not (<nmod_mpoly>other).ctx:
+            raise IncompatibleContextError(f"{(<nmod_mpoly>self).ctx} is not {(<nmod_mpoly>other).ctx}")
+        res = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_gcd(res.val, (<nmod_mpoly>self).val, (<nmod_mpoly>other).val, res.ctx.val)
+        return res
+
+    def sqrt(self):
+        """
+        Return the square root of self.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> f = ctx.from_dict({(1, 1): 4, (0, 0): 1})
+            >>> (f * f).sqrt()
+            4*x0*x1 + 1
+        """
+        cdef nmod_mpoly res
+        if not self.ctx.is_prime():
+            raise DomainError("square root with non-prime modulus is not supported")
+
+        res = create_nmod_mpoly(self.ctx)
+
+        if nmod_mpoly_sqrt(res.val, self.val, self.ctx.val):
+            return res
+        else:
+            raise ValueError("polynomial is not a perfect square")
+
+    def factor(self):
+        """
+        Factors self into irreducible factors, returning a tuple
+        (c, factors) where c is the content of the coefficients and
+        factors is a list of (poly, exp) pairs.
+
+            >>> from flint import Ordering
+            >>> Zm = nmod_mpoly
+            >>> ctx = nmod_mpoly_ctx.get_context(3, Ordering.lex, 11, 'x,y,z')
+            >>> p1 = Zm("2*x + 4", ctx)
+            >>> p2 = Zm("3*x*z + 3*x + 3*z + 3", ctx)
+            >>> (p1 * p2).factor()
+            (6, [(z + 1, 1), (x + 2, 1), (x + 1, 1)])
+            >>> (p2 * p1 * p2).factor()
+            (18, [(z + 1, 2), (x + 2, 1), (x + 1, 2)])
+        """
+        cdef:
+            nmod_mpoly_factor_t fac
+            fmpz c
+            nmod_mpoly u
+
+        if not self.ctx.is_prime():
+            raise DomainError("factorisation with non-prime modulus is not supported")
+
+        nmod_mpoly_factor_init(fac, self.ctx.val)
+        if not nmod_mpoly_factor(fac, self.val, self.ctx.val):
+            raise RuntimeError("factorisation failed")
+        res = [0] * fac.num
+
+        for i in range(fac.num):
+            u = create_nmod_mpoly(self.ctx)
+            nmod_mpoly_set((<nmod_mpoly>u).val, &fac.poly[i], self.ctx.val)
+
+            c = fmpz.__new__(fmpz)
+            fmpz_set((<fmpz>c).val, &fac.exp[i])
+
+            res[i] = (u, c)
+
+        constant = fac.constant
+        nmod_mpoly_factor_clear(fac, self.ctx.val)
+        return constant, res
+
+    def factor_squarefree(self):
+        """
+        Factors self into irreducible factors, returning a tuple
+        (c, factors) where c is the content of the coefficients and
+        factors is a list of (poly, exp) pairs.
+
+            >>> from flint import Ordering
+            >>> Zm = nmod_mpoly
+            >>> ctx = nmod_mpoly_ctx.get_context(3, Ordering.lex, 11, 'x,y,z')
+            >>> p1 = Zm("2*x + 4", ctx)
+            >>> p2 = Zm("3*x*y + 3*x + 3*y + 3", ctx)
+            >>> (p1 * p2).factor_squarefree()
+            (6, [(y + 1, 1), (x^2 + 3*x + 2, 1)])
+            >>> (p1 * p2 * p1).factor_squarefree()
+            (12, [(y + 1, 1), (x + 1, 1), (x + 2, 2)])
+        """
+        cdef:
+            nmod_mpoly_factor_t fac
+            fmpz c
+            nmod_mpoly u
+
+        if not self.ctx.is_prime():
+            raise DomainError("factorisation with non-prime modulus is not supported")
+
+        nmod_mpoly_factor_init(fac, self.ctx.val)
+        if not nmod_mpoly_factor_squarefree(fac, self.val, self.ctx.val):
+            raise RuntimeError("factorisation failed")
+        res = [0] * fac.num
+
+        for i in range(fac.num):
+            u = create_nmod_mpoly(self.ctx)
+            nmod_mpoly_init(u.val, u.ctx.val)
+            nmod_mpoly_set((<nmod_mpoly>u).val, &fac.poly[i], self.ctx.val)
+
+            c = fmpz.__new__(fmpz)
+            fmpz_set((<fmpz>c).val, &fac.exp[i])
+
+            res[i] = (u, c)
+
+        constant = fac.constant
+        nmod_mpoly_factor_clear(fac, self.ctx.val)
+        return constant, res
+
+    # TODO: Rethink context conversions, particularly the proposed methods in #132
+    # def coerce_to_context(self, ctx):
+    #     cdef:
+    #         nmod_mpoly res
+    #         slong *C
+    #         slong i
+
+    #     if not typecheck(ctx, nmod_mpoly_ctx):
+    #         raise ValueError("provided context is not a nmod_mpoly_ctx")
+
+    #     if self.ctx is ctx:
+    #         return self
+
+    #     C = <slong *> libc.stdlib.malloc(self.ctx.val.minfo.nvars * sizeof(slong *))
+    #     if C is NULL:
+    #         raise MemoryError("malloc returned a null pointer")
+    #     res = create_nmod_mpoly(self.ctx)
+
+    #     vars = {x: i for i, x in enumerate(ctx.py_names)}
+    #     for i, var in enumerate(self.ctx.py_names):
+    #         C[i] = <slong>vars[var]
+
+    #     nmod_mpoly_compose_nmod_mpoly_gen(res.val, self.val, C, self.ctx.val, (<nmod_mpoly_ctx>ctx).val)
+
+    #     libc.stdlib.free(C)
+    #     return res
+
+    def derivative(self, var):
+        """
+        Return the derivative of this polynomial with respect to the provided variable.
+        The argument can either be the variable as a string, or the index of the
+        variable in the context.
+
+            >>> from flint import Ordering
+            >>> ctx = nmod_mpoly_ctx.get_context(2, Ordering.lex, 11, 'x')
+            >>> p = ctx.from_dict({(0, 3): 2, (2, 1): 3})
+            >>> p
+            3*x0^2*x1 + 2*x1^3
+            >>> p.derivative("x0")
+            6*x0*x1
+            >>> p.derivative(1)
+            3*x0^2 + 6*x1^2
+
+        """
+        cdef:
+            nmod_mpoly res
+            slong i = self.ctx.variable_to_index(var)
+
+        res = create_nmod_mpoly(self.ctx)
+
+        nmod_mpoly_derivative(res.val, self.val, i, self.ctx.val)
+        return res
+
+
+cdef class nmod_mpoly_vec:
+    """
+    A class representing a vector of nmod_mpolys.
+    """
+
+    def __cinit__(self, iterable_or_len, nmod_mpoly_ctx ctx, bint double_indirect = False):
+        if isinstance(iterable_or_len, int):
+            self.length = iterable_or_len
+        else:
+            self.length = len(iterable_or_len)
+
+        self.ctx = ctx
+
+        self.val = <nmod_mpoly_struct *> libc.stdlib.malloc(self.length * sizeof(nmod_mpoly_struct))
+        for i in range(self.length):
+            nmod_mpoly_init(&self.val[i], self.ctx.val)
+
+        if double_indirect:
+            self.double_indirect = <nmod_mpoly_struct **> libc.stdlib.malloc(self.length * sizeof(nmod_mpoly_struct *))
+            if self.double_indirect is NULL:
+                raise MemoryError("malloc returned a null pointer")  # pragma: no cover
+
+            for i in range(self.length):
+                self.double_indirect[i] = &self.val[i]
+        else:
+            self.double_indirect = NULL
+
+    def __init__(self, iterable_or_len, _, double_indirect: bool = False):
+        if not isinstance(iterable_or_len, int):
+            for i, x in enumerate(iterable_or_len):
+                self[i] = x
+
+    def __dealloc__(self):
+        libc.stdlib.free(self.double_indirect)
+        for i in range(self.length):
+            nmod_mpoly_clear(&self.val[i], self.ctx.val)
+        libc.stdlib.free(self.val)
+
+    def __getitem__(self, x):
+        if not isinstance(x, int):
+            raise TypeError("index is not integer")
+        elif not 0 <= x < self.length:
+            raise IndexError("index out of range")
+
+        cdef nmod_mpoly z = create_nmod_mpoly(self.ctx)
+        nmod_mpoly_set(z.val, &self.val[x], self.ctx.val)
+        return z
+
+    def __setitem__(self, x, y):
+        if not isinstance(x, int):
+            raise TypeError("index is not integer")
+        elif not 0 <= x < self.length:
+            raise IndexError("index out of range")
+        elif not typecheck(y, nmod_mpoly):
+            raise TypeError("argument is not nmod_mpoly")
+        elif (<nmod_mpoly>y).ctx is not self.ctx:
+            raise IncompatibleContextError(f"{(<nmod_mpoly>y).ctx} is not {self.ctx}")
+
+        nmod_mpoly_set(&self.val[x], (<nmod_mpoly>y).val, self.ctx.val)
+
+    def __len__(self):
+        return self.val.length
+
+    def __str__(self):
+        s = [None] * self.length
+        for i in range(self.length):
+            x = create_nmod_mpoly(self.ctx)
+            nmod_mpoly_set(x.val, &self.val[x], self.ctx.val)
+            s[i] = str(x)
+        return f"[{', '.join(s)}]"
+
+    def __repr__(self):
+        return f"nmod_mpoly_vec({self}, ctx={self.ctx})"
+
+    def to_tuple(self):
+        return tuple(self[i] for i in range(self.val.length))