diff --git a/doc/source/fmpz_mod.rst b/doc/source/fmpz_mod.rst
new file mode 100644
index 00000000..67e3206d
--- /dev/null
+++ b/doc/source/fmpz_mod.rst
@@ -0,0 +1,13 @@
+**fmpz_mod** -- integers mod n
+===============================================================================
+
+.. autoclass :: flint.fmpz_mod_ctx
+  :members:
+  :inherited-members:
+  :undoc-members:
+
+.. autoclass :: flint.fmpz_mod
+  :members:
+  :inherited-members:
+  :undoc-members:
+
diff --git a/doc/source/general.rst b/doc/source/general.rst
index 6462df8a..52156f0f 100644
--- a/doc/source/general.rst
+++ b/doc/source/general.rst
@@ -161,8 +161,6 @@ determined from the available data.
 The following convenience functions are provided for numerical evaluation
 with adaptive working precision.
 
-.. autofunction :: flint.good
-
 .. autofunction :: flint.showgood
 
 Power series
diff --git a/doc/source/index.rst b/doc/source/index.rst
index bc7c7ae1..193e5701 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -43,6 +43,7 @@ Scalar types
 
    fmpz.rst
    fmpq.rst
+   fmpz_mod.rst
    nmod.rst
    arb.rst
    acb.rst
diff --git a/doc/source/nmod.rst b/doc/source/nmod.rst
index 4f8bda48..febc625a 100644
--- a/doc/source/nmod.rst
+++ b/doc/source/nmod.rst
@@ -1,4 +1,4 @@
-**nmod** -- integers mod n
+**nmod** -- integers mod wordsize n
 ===============================================================================
 
 .. autoclass :: flint.nmod
diff --git a/setup.py b/setup.py
index d7df0a0b..1972be10 100644
--- a/setup.py
+++ b/setup.py
@@ -94,6 +94,7 @@
     ("flint.types.acb_mat", ["src/flint/types/acb_mat.pyx"]),
     ("flint.types.acb_series", ["src/flint/types/acb_series.pyx"]),
     ("flint.types.fmpz_mpoly", ["src/flint/types/fmpz_mpoly.pyx"]),
+    ("flint.types.fmpz_mod", ["src/flint/types/fmpz_mod.pyx"]),
     ("flint.types.dirichlet", ["src/flint/types/dirichlet.pyx"]),
     ("flint.flint_base.flint_base", ["src/flint/flint_base/flint_base.pyx"]),
     ("flint.flint_base.flint_context", ["src/flint/flint_base/flint_context.pyx"]),
diff --git a/src/flint/__init__.py b/src/flint/__init__.py
index 74b306a0..ef3f78e9 100644
--- a/src/flint/__init__.py
+++ b/src/flint/__init__.py
@@ -21,6 +21,8 @@
 from .types.acb_mat import *
 from .types.acb_series import *
 from .types.fmpz_mpoly import *
+from .types.fmpz_mod import *
+from .types.dirichlet import *
 from .functions.showgood import showgood
 
 __version__ = '0.4.4'
diff --git a/src/flint/flintlib/fmpz.pxd b/src/flint/flintlib/fmpz.pxd
index 4779537d..4decaaa9 100644
--- a/src/flint/flintlib/fmpz.pxd
+++ b/src/flint/flintlib/fmpz.pxd
@@ -1,11 +1,14 @@
-from flint.flintlib.flint cimport fmpz_struct, ulong, mp_limb_t
+from flint.flintlib.flint cimport fmpz_struct, ulong, mp_limb_t, mp_ptr
 from flint.flintlib.flint cimport mp_size_t, mp_bitcnt_t, slong, flint_rand_t, flint_bitcnt_t
-# from flint.flintlib.nmod cimport nmod_t
-# from flint.flintlib.fmpz_factor cimport fmpz_factor_t
 
 cdef extern from "flint/fmpz.h":
     ctypedef fmpz_struct fmpz_t[1]
 
+    ctypedef struct fmpz_preinvn_struct:
+        mp_ptr dinv
+        slong n
+        flint_bitcnt_t norm
+    ctypedef fmpz_preinvn_struct fmpz_preinvn_t[1]
 
 # from here on is parsed
     # fmpz_struct PTR_TO_COEFF(__mpz_struct * ptr)
diff --git a/src/flint/flintlib/fmpz_mod.pxd b/src/flint/flintlib/fmpz_mod.pxd
new file mode 100644
index 00000000..a96a0e3a
--- /dev/null
+++ b/src/flint/flintlib/fmpz_mod.pxd
@@ -0,0 +1,45 @@
+from flint.flintlib.flint cimport ulong, slong
+from flint.flintlib.fmpz cimport fmpz_t, fmpz_preinvn_struct
+from flint.flintlib.nmod cimport nmod_t
+
+# unimported types  {'fmpz_mod_discrete_log_pohlig_hellman_t'}
+
+cdef extern from "flint/fmpz_mod.h":
+    ctypedef struct fmpz_mod_ctx_struct:
+        fmpz_t n
+        nmod_t mod
+        ulong n_limbs[3]
+        ulong ninv_limbs[3]
+        fmpz_preinvn_struct * ninv_huge
+    ctypedef fmpz_mod_ctx_struct fmpz_mod_ctx_t[1]
+
+    # Parsed from here
+    void fmpz_mod_ctx_init(fmpz_mod_ctx_t ctx, const fmpz_t n)
+    void fmpz_mod_ctx_clear(fmpz_mod_ctx_t ctx)
+    void fmpz_mod_ctx_set_modulus(fmpz_mod_ctx_t ctx, const fmpz_t n)
+    void fmpz_mod_set_fmpz(fmpz_t a, const fmpz_t b, const fmpz_mod_ctx_t ctx)
+    int fmpz_mod_is_canonical(const fmpz_t a, const fmpz_mod_ctx_t ctx)
+    int fmpz_mod_is_one(const fmpz_t a, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_add(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_add_fmpz(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_add_ui(fmpz_t a, const fmpz_t b, ulong c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_add_si(fmpz_t a, const fmpz_t b, slong c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_sub(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_sub_fmpz(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_sub_ui(fmpz_t a, const fmpz_t b, ulong c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_sub_si(fmpz_t a, const fmpz_t b, slong c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_fmpz_sub(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_ui_sub(fmpz_t a, ulong b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_si_sub(fmpz_t a, slong b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_neg(fmpz_t a, const fmpz_t b, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_mul(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_inv(fmpz_t a, const fmpz_t b, const fmpz_mod_ctx_t ctx)
+    int fmpz_mod_divides(fmpz_t a, const fmpz_t b, const fmpz_t c, const fmpz_mod_ctx_t ctx)
+    void fmpz_mod_pow_ui(fmpz_t a, const fmpz_t b, ulong e, const fmpz_mod_ctx_t ctx)
+    int fmpz_mod_pow_fmpz(fmpz_t a, const fmpz_t b, const fmpz_t e, const fmpz_mod_ctx_t ctx)
+    # void fmpz_mod_discrete_log_pohlig_hellman_init(fmpz_mod_discrete_log_pohlig_hellman_t L)
+    # void fmpz_mod_discrete_log_pohlig_hellman_clear(fmpz_mod_discrete_log_pohlig_hellman_t L)
+    # double fmpz_mod_discrete_log_pohlig_hellman_precompute_prime(fmpz_mod_discrete_log_pohlig_hellman_t L, const fmpz_t p)
+    # const fmpz_struct * fmpz_mod_discrete_log_pohlig_hellman_primitive_root(const fmpz_mod_discrete_log_pohlig_hellman_t L)
+    # void fmpz_mod_discrete_log_pohlig_hellman_run(fmpz_t x, const fmpz_mod_discrete_log_pohlig_hellman_t L, const fmpz_t y)
+    int fmpz_next_smooth_prime(fmpz_t a, const fmpz_t b)
diff --git a/src/flint/test/test.py b/src/flint/test/test.py
index 9c9dc8fe..8121d1ed 100644
--- a/src/flint/test/test.py
+++ b/src/flint/test/test.py
@@ -1584,6 +1584,192 @@ def test_pickling():
         obj2 = pickle.loads(s)
         assert obj == obj2
 
+def test_fmpz_mod():
+    from flint import fmpz_mod_ctx, fmpz, fmpz_mod
+    
+    p_sml = 163
+    p_med = 2**127 - 1
+    p_big = 2**255 - 19
+
+    F_sml = fmpz_mod_ctx(p_sml)
+    F_med = fmpz_mod_ctx(p_med)
+    F_big = fmpz_mod_ctx(p_big)
+
+    # Context tests
+    assert raises(lambda: fmpz_mod_ctx("AAA"), TypeError)
+    assert raises(lambda: fmpz_mod_ctx(-1), ValueError)
+    assert F_sml.modulus() == p_sml
+    assert F_med.modulus() == p_med
+    assert F_big.modulus() == p_big
+
+    F_big_copy = fmpz_mod_ctx(p_big)
+    assert F_big_copy == F_big
+    assert F_big != F_sml
+    assert hash(F_big_copy) == hash(F_big)
+    assert hash(F_big) != hash(F_sml)
+    assert F_big_copy != F_sml
+    assert F_big_copy != "A"
+
+    assert repr(F_sml) == "fmpz_mod_ctx(163)"
+    assert str(F_sml) == "Context for fmpz_mod with modulus: 163"
+
+    # Type tests
+    assert raises(lambda: fmpz_mod(1, "AAA"), TypeError)
+
+    # Test for small, medium and large char.
+    for F_test in [F_sml, F_med, F_big]:
+        test_mod = int(F_test.modulus())
+        test_x = (-123) % test_mod # canonical value
+        test_y = ((-456) % test_mod)**2 # non-canoncial value
+
+        F_test_copy = fmpz_mod_ctx(test_mod)
+        F_other = fmpz_mod_ctx(11)
+
+        assert raises(lambda: F_test(test_x) > 0, TypeError)
+        assert raises(lambda: F_test(test_x) >= 0, TypeError)
+        assert raises(lambda: F_test(test_x) < 0, TypeError)
+        assert raises(lambda: F_test(test_x) <= 0, TypeError)
+
+        assert (test_x == F_test(test_x)) is True, f"{test_x}, {F_test(test_x)}"
+        assert (124 != F_test(test_x)) is True
+        assert (F_test(test_x) == test_x) is True
+        assert (F_test(test_x) == test_x + test_mod) is True
+        assert (F_test(test_x) == 1) is False
+        assert (F_test(test_x) != 1) is True
+        assert (F_test(test_x) == F_test(test_x)) is True
+        assert (F_test(test_x) == F_test(test_x + test_mod)) is True
+        assert (F_test(test_x) == F_test(1)) is False
+        assert (F_test(test_x) != F_test(1)) is True
+
+        assert (hash(F_test(test_x)) == hash(test_x)) is True
+        assert (hash(F_test(F_test(test_x))) == hash(test_x)) is True
+        assert (hash(F_test(test_x)) == hash(1)) is False
+        assert (hash(F_test(test_x)) != hash(1)) is True
+        assert (hash(F_test(test_x)) == hash(F_test(test_x))) is True
+        assert (hash(F_test(test_x)) == hash(F_test(test_x + test_mod))) is True
+        assert (hash(F_test(test_x)) == hash(F_test(1))) is False
+        assert (hash(F_test(test_x)) != hash(F_test(1))) is True
+
+        # Is one, zero
+        assert (F_test(0) == 0) is True
+        assert F_test(0).is_zero() is True
+        assert not F_test(0)
+        assert not F_test(test_mod)
+        assert F_test(1).is_one() is True
+        assert F_test(test_mod + 1).is_one() is True
+        assert F_test(1).is_one() is True
+        assert F_test(2).is_one() is False
+
+        # int, str, repr
+        assert str(F_test(11)) == "11"
+        assert str(F_test(-1)) == str(test_mod - 1)
+        assert repr(F_test(11)) == f"fmpz_mod(11, {test_mod})"
+        assert repr(F_test(-1)) == f"fmpz_mod({test_mod - 1}, {test_mod})"
+
+        assert +F_test(5) == F_test(5)
+
+        # Arithmetic tests
+
+        # Negation
+        assert -F_test(test_x) == F_test(-test_x) == (-test_x % test_mod)
+        assert -F_test(1) == F_test(-1) == F_test(test_mod - 1)
+
+        # Addition
+        assert F_test(test_x) + F_test(test_y) == F_test(test_x + test_y)
+        assert F_test(test_x) + F_test_copy(test_y) == F_test(test_x + test_y)
+        assert F_test(test_x) + F_test(test_y) == F_test_copy(test_x + test_y)
+        assert raises(lambda: F_test(test_x) + "AAA", TypeError)
+
+        assert F_test(test_x) + F_test(test_y) == F_test(test_y) + F_test(test_x)
+        assert F_test(test_x) + test_y == F_test(test_x + test_y)
+        assert test_y + F_test(test_x) == F_test(test_x + test_y)
+        assert F_test(test_x) + fmpz(test_y) == F_test(test_y) + F_test(test_x)
+        assert raises(lambda: F_test(test_x) + F_other(test_y), ValueError)
+
+        # Subtraction
+
+        assert F_test(test_x) - F_test(test_y) == F_test(test_x - test_y)
+        assert F_test(test_x) - test_y == F_test(test_x - test_y)
+        assert F_test(test_x) - test_y == F_test(test_x) - F_test(test_y)
+        assert F_test(test_y) - test_x == F_test(test_y) - F_test(test_x)
+        assert test_x - F_test(test_y) == F_test(test_x) - F_test(test_y)
+        assert test_y - F_test(test_x) == F_test(test_y) - F_test(test_x)
+        assert F_test(test_x) - fmpz(test_y) == F_test(test_x) - F_test(test_y)
+        assert raises(lambda: F_test(test_x) - F_other(test_y), ValueError)
+        assert raises(lambda: F_test(test_x) - "AAA", TypeError)
+
+        # Multiplication
+
+        assert F_test(test_x) * F_test(test_y) == (test_x * test_y) % test_mod
+        assert F_test(test_x) * test_y == (test_x * test_y) % test_mod
+        assert test_y * F_test(test_x) == (test_x * test_y) % test_mod
+
+        assert F_test(1) * F_test(test_x) == F_test(1 * test_x)
+        assert F_test(2) * F_test(test_x) == F_test(2 * test_x)
+        assert F_test(3) * F_test(test_x) == F_test(3 * test_x)
+        assert 1 * F_test(test_x) == F_test(1 * test_x)
+        assert 2 * F_test(test_x) == F_test(2 * test_x)
+        assert 3 * F_test(test_x) == F_test(3 * test_x)
+        assert F_test(test_x) * 1 == F_test(1 * test_x)
+        assert F_test(test_x) * 2 == F_test(2 * test_x)
+        assert F_test(test_x) * 3 == F_test(3 * test_x)
+        assert fmpz(1) * F_test(test_x) == F_test(1 * test_x)
+        assert fmpz(2) * F_test(test_x) == F_test(2 * test_x)
+        assert fmpz(3) * F_test(test_x) == F_test(3 * test_x)
+        assert raises(lambda: F_test(test_x) * "AAA", TypeError)
+        assert raises(lambda: F_test(test_x) * F_other(test_x), ValueError)
+
+        # Exponentiation 
+
+        assert F_test(0)**0 == pow(0, 0, test_mod)
+        assert F_test(0)**1 == pow(0, 1, test_mod)
+        assert F_test(0)**2 == pow(0, 2, test_mod)
+        assert raises(lambda: F_test(0)**(-1), ZeroDivisionError)
+        assert raises(lambda: F_test(0)**("AA"), NotImplementedError)
+
+        assert F_test(test_x)**fmpz(0) == pow(test_x, 0, test_mod)
+        assert F_test(test_x)**fmpz(1) == pow(test_x, 1, test_mod)
+        assert F_test(test_x)**fmpz(2) == pow(test_x, 2, test_mod)
+        assert F_test(test_x)**fmpz(3) == pow(test_x, 3, test_mod)
+
+        assert F_test(test_x)**0 == pow(test_x, 0, test_mod)
+        assert F_test(test_x)**1 == pow(test_x, 1, test_mod)
+        assert F_test(test_x)**2 == pow(test_x, 2, test_mod)
+        assert F_test(test_x)**3 == pow(test_x, 3, test_mod)
+        assert F_test(test_x)**100 == pow(test_x, 100, test_mod)
+
+        assert F_test(test_x)**(-1) == pow(test_x, -1, test_mod)
+        assert F_test(test_x)**(-2) == pow(test_x, -2, test_mod)
+        assert F_test(test_x)**(-3) == pow(test_x, -3, test_mod)
+        assert F_test(test_x)**(-4) == pow(test_x, -4, test_mod)
+
+        # Inversion
+
+        assert raises(lambda: ~F_test(0), ZeroDivisionError)
+        assert ~F_test(test_x) == pow(test_x, -1, test_mod)
+        assert ~F_test(1) == pow(1, -1, test_mod)
+        assert ~F_test(2) == pow(2, -1, test_mod), f"Broken!! {~F_test(2)}, {pow(2, -1, test_mod)}"
+
+        assert F_test(1).inverse(check=False) == pow(1, -1, test_mod)
+        assert F_test(2).inverse(check=False) == pow(2, -1, test_mod)
+        assert F_test(test_x).inverse(check=False) == pow(test_x, -1, test_mod)
+
+        # Division
+        assert raises(lambda: F_test(1) / F_test(0), ZeroDivisionError)
+        assert F_test(test_x) / F_test(test_y) == (test_x * pow(test_y, -1, test_mod)) % test_mod
+        assert F_test(test_x) / fmpz(test_y) == (test_x * pow(test_y, -1, test_mod)) % test_mod
+        assert F_test(test_x) / test_y == (test_x * pow(test_y, -1, test_mod)) % test_mod
+        assert raises(lambda: F_test(test_x) / "AAA", TypeError)
+        assert raises(lambda: "AAA" / F_test(test_x), TypeError)
+        assert raises(lambda: F_other(test_x) / F_test(test_x), ValueError)
+        assert raises(lambda: F_test(test_x) // F_test(test_x), TypeError)
+        assert 1 / F_test(2) ==  pow(2, -1, test_mod)
+        assert 1 / F_test(test_x) ==  pow(test_x, -1, test_mod)
+        assert 1 / F_test(test_y) ==  pow(test_y, -1, test_mod)
+
+        assert fmpz(test_y) / F_test(test_x) == (test_y * pow(test_x, -1, test_mod)) % test_mod
+        assert test_y / F_test(test_x) == (test_y * pow(test_x, -1, test_mod)) % test_mod
+            
 
 all_tests = [
     test_pyflint,
@@ -1603,4 +1789,5 @@ def test_pickling():
     test_nmod_poly,
     test_nmod_mat,
     test_arb,
+    test_fmpz_mod,
 ]
diff --git a/src/flint/types/fmpz.pxd b/src/flint/types/fmpz.pxd
index 4317c89f..074aa411 100644
--- a/src/flint/types/fmpz.pxd
+++ b/src/flint/types/fmpz.pxd
@@ -10,6 +10,7 @@ from flint.flintlib.fmpz cimport fmpz_t, fmpz_set_str, fmpz_set_si
 from cpython.version cimport PY_MAJOR_VERSION
 
 cdef int fmpz_set_any_ref(fmpz_t x, obj)
+cdef fmpz_get_intlong(fmpz_t x)
 
 cdef inline int fmpz_set_pylong(fmpz_t x, obj):
     cdef int overflow
diff --git a/src/flint/types/fmpz_mod.pxd b/src/flint/types/fmpz_mod.pxd
new file mode 100644
index 00000000..264d35d7
--- /dev/null
+++ b/src/flint/types/fmpz_mod.pxd
@@ -0,0 +1,14 @@
+from flint.flint_base.flint_base cimport flint_scalar
+from flint.flintlib.fmpz cimport fmpz_t
+from flint.flintlib.fmpz_mod cimport fmpz_mod_ctx_t
+
+
+cdef class fmpz_mod_ctx:
+    cdef fmpz_mod_ctx_t val
+
+cdef class fmpz_mod(flint_scalar):
+    cdef fmpz_mod_ctx ctx
+    cdef fmpz_t val
+
+    cdef any_as_fmpz_mod(self, obj)
+
diff --git a/src/flint/types/fmpz_mod.pyx b/src/flint/types/fmpz_mod.pyx
new file mode 100644
index 00000000..362a9161
--- /dev/null
+++ b/src/flint/types/fmpz_mod.pyx
@@ -0,0 +1,364 @@
+from flint.flintlib.fmpz cimport (
+    fmpz_t,
+    fmpz_one,
+    fmpz_set,
+    fmpz_init,
+    fmpz_clear,
+    fmpz_equal
+)
+from flint.flintlib.fmpz_mod cimport *
+
+from flint.utils.typecheck cimport typecheck
+from flint.flint_base.flint_base cimport flint_scalar
+from flint.types.fmpz cimport (
+    fmpz,
+    any_as_fmpz,
+    fmpz_get_intlong
+)
+
+
+cdef class fmpz_mod_ctx:
+    r"""
+    Context object for creating :class:`~.fmpz_mod` initalised 
+    with a modulus :math:`N`.
+
+        >>> fmpz_mod_ctx(2**127 - 1)
+        fmpz_mod_ctx(170141183460469231731687303715884105727)
+
+    """
+    def __cinit__(self):
+        cdef fmpz one = fmpz.__new__(fmpz)
+        fmpz_one(one.val)
+        fmpz_mod_ctx_init(self.val, one.val)
+
+    def __dealloc__(self):
+        fmpz_mod_ctx_clear(self.val)
+
+    def __init__(self, mod):
+        # Ensure modulus is fmpz type
+        if not typecheck(mod, fmpz):
+            mod = any_as_fmpz(mod)
+            if mod is NotImplemented:
+                raise TypeError("Context modulus must be able to be case to an `fmpz` type")
+
+        # Ensure modulus is positive
+        if mod < 1:
+            raise ValueError("Modulus is expected to be positive")
+
+        # Init the context
+        fmpz_mod_ctx_init(self.val, (<fmpz>mod).val)
+    
+    def modulus(self):
+        """
+        Return the modulus from the context as an fmpz
+        type
+
+            >>> mod_ctx = fmpz_mod_ctx(2**127 - 1)
+            >>> mod_ctx.modulus()
+            170141183460469231731687303715884105727
+
+        """
+        n = fmpz()
+        fmpz_set(n.val, (<fmpz_t>self.val.n))
+        return n
+
+    def __eq__(self, other):
+        if typecheck(other, fmpz_mod_ctx):
+            return fmpz_equal(self.val.n, (<fmpz_mod_ctx>other).val.n)
+        return False
+
+    def __hash__(self):
+        return hash(self.modulus())
+
+    def __str__(self):
+        return f"Context for fmpz_mod with modulus: {self.modulus()}"
+
+    def __repr__(self):
+        return f"fmpz_mod_ctx({self.modulus()})"
+
+    def __call__(self, val):
+        return fmpz_mod(val, self)
+
+cdef class fmpz_mod(flint_scalar):
+    """
+    The *fmpz_mod* type represents integer modulo an 
+    arbitrary-size modulus. For wordsize modulus, see
+    :class:`~.nmod`.
+
+    An *fmpz_mod* element is constructed from an :class:`~.fmpz_mod_ctx`
+    either by passing it as an argument to the type, or
+    by directly calling the context
+
+        >>> fmpz_mod(-1, fmpz_mod_ctx(2**127 - 1))
+        fmpz_mod(170141183460469231731687303715884105726, 170141183460469231731687303715884105727)
+        >>> ZmodN = fmpz_mod_ctx(2**127 - 1)
+        >>> ZmodN(-2)
+        fmpz_mod(170141183460469231731687303715884105725, 170141183460469231731687303715884105727)
+
+    """
+
+    def __cinit__(self):
+        fmpz_init(self.val)
+
+    def __dealloc__(self):
+        fmpz_clear(self.val)
+
+    def __init__(self, val, ctx):
+        if not typecheck(ctx, fmpz_mod_ctx):
+            raise TypeError
+        self.ctx = ctx
+
+        # When the input is also an fmpz_mod we just need
+        # moduli to match
+        if typecheck(val, fmpz_mod):
+            if self.ctx != (<fmpz_mod>val).ctx:
+                raise ValueError("moduli must match")
+            # fmpz_mod_set_fmpz(self.val, (<fmpz_mod>val).val, self.ctx.val)
+            fmpz_set(self.val, (<fmpz_mod>val).val)
+            return
+
+        # For all other cases, the easiest is to first convert to
+        # fmpz type and set this way
+        if not typecheck(val, fmpz):
+            val = any_as_fmpz(val)
+            if val is NotImplemented:
+                raise NotImplementedError
+        fmpz_mod_set_fmpz(self.val, (<fmpz>val).val, self.ctx.val)
+
+    cdef any_as_fmpz_mod(self, obj):
+        try:
+            return self.ctx(obj)
+        except NotImplementedError:
+            return NotImplemented
+
+    def is_zero(self):
+        """
+        Return whether an element is equal to zero
+
+            >>> mod_ctx = fmpz_mod_ctx(163)
+            >>> mod_ctx(0).is_zero()
+            True
+            >>> mod_ctx(1).is_zero()
+            False
+        """
+        return self == 0
+    
+    def is_one(self):
+        """
+        Return whether an element is equal to one
+
+            >>> mod_ctx = fmpz_mod_ctx(163)
+            >>> mod_ctx(0).is_one()
+            False
+            >>> mod_ctx(1).is_zero()
+            True
+        """
+
+        cdef bint res
+        res = fmpz_mod_is_one(self.val, self.ctx.val)
+        return res == 1
+
+    def __richcmp__(self, other, int op):
+        cdef bint res
+        if op != 2 and op != 3:
+            raise TypeError("fmpz_mod cannot be ordered")
+
+        if not typecheck(other, fmpz_mod):
+            other = self.any_as_fmpz_mod(other)
+
+        if typecheck(self, fmpz_mod) and typecheck(other, fmpz_mod):
+            res = fmpz_equal(self.val, (<fmpz_mod>other).val) and \
+                  (self.ctx == (<fmpz_mod>other).ctx)
+            if op == 2:
+                return res
+            else:
+                return not res
+        else:
+            return NotImplemented
+
+    def __bool__(self):
+        return not self.is_zero()
+
+    def __repr__(self):
+        return "fmpz_mod({}, {})".format(
+            fmpz_get_intlong(self.val),
+            self.ctx.modulus()
+        )
+
+    def __hash__(self):
+        return  hash((int(self)))
+
+    def __int__(self):
+        return fmpz_get_intlong(self.val)
+
+    def __str__(self):
+        return str(int(self))
+
+    # ---------------- #
+    #    Arithmetic    #
+    # ---------------- #
+
+    def __pos__(self):
+        return self
+
+    def __neg__(self):
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        res.ctx = self.ctx
+        fmpz_mod_neg(res.val, self.val, self.ctx.val)
+        return res
+
+    def __add__(self, other):
+        other = self.any_as_fmpz_mod(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        res.ctx = self.ctx
+        fmpz_mod_add(
+            res.val, self.val, (<fmpz_mod>other).val, self.ctx.val
+        )
+        return res
+
+    def __radd__(self, other):
+        return self.__add__(other)
+
+    def __sub__(self, other):
+        other = self.any_as_fmpz_mod(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        res.ctx = self.ctx
+
+        fmpz_mod_sub(
+            res.val, self.val, (<fmpz_mod>other).val, self.ctx.val
+        )
+        return res
+
+    def __rsub__(self, other):
+        return self.__sub__(other).__neg__()
+
+    def __mul__(self, other):
+        other = self.any_as_fmpz_mod(other)
+        if other is NotImplemented:
+            return NotImplemented
+
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        res.ctx = self.ctx
+
+        fmpz_mod_mul(
+            res.val, self.val, (<fmpz_mod>other).val, self.ctx.val
+        )
+        return res
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
+    @staticmethod
+    def _div_(left, right):
+        cdef bint check
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        
+        # Division when left and right are fmpz_mod
+        if typecheck(left, fmpz_mod) and typecheck(right, fmpz_mod):
+            res.ctx = (<fmpz_mod>left).ctx
+            if not (<fmpz_mod>left).ctx == (<fmpz_mod>right).ctx:
+                raise ValueError("moduli must match")
+            check = fmpz_mod_divides(
+                res.val, (<fmpz_mod>left).val, (<fmpz_mod>right).val, res.ctx.val
+            ) 
+        
+        # Case when only left is fmpz_mod
+        elif typecheck(left, fmpz_mod):
+            res.ctx = (<fmpz_mod>left).ctx
+            right = any_as_fmpz(right)
+            if right is NotImplemented:
+                return NotImplemented   
+            check = fmpz_mod_divides(
+                res.val, (<fmpz_mod>left).val, (<fmpz>right).val, res.ctx.val
+            ) 
+
+        # Case when right is an fmpz_mod
+        else:
+            res.ctx = (<fmpz_mod>right).ctx
+            left = any_as_fmpz(left)
+            if left is NotImplemented:
+                return NotImplemented   
+            check = fmpz_mod_divides(
+                res.val, (<fmpz>left).val, (<fmpz_mod>right).val, res.ctx.val
+            ) 
+        
+        if check == 0:
+            raise ZeroDivisionError(f"{right} is not invertible modulo {res.ctx.modulus()}")
+
+        return res
+
+    def __truediv__(s, t):
+        return fmpz_mod._div_(s, t)
+
+    def __rtruediv__(s, t):
+        return fmpz_mod._div_(t, s)
+
+    def __floordiv__(self, other):
+        return NotImplemented
+
+    def inverse(self, check=True):
+        r"""
+        Computes :math:`a^{-1} \pmod N`
+
+        When check=False, the solutions is assumed to exist and Flint will abort on
+        failure. 
+
+            >>> mod_ctx = fmpz_mod_ctx(163)
+            >>> mod_ctx(2).inverse()
+            fmpz_mod(82, 163)
+            >>> mod_ctx(2).inverse(check=False)
+            fmpz_mod(82, 163)
+        """
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        res.ctx = self.ctx
+
+        if check is False:
+            fmpz_mod_inv(res.val, self.val, self.ctx.val)
+            return res
+
+        cdef bint r
+        cdef fmpz one = fmpz.__new__(fmpz)
+        fmpz_one(one.val)
+
+        r = fmpz_mod_divides(
+            res.val, one.val, self.val, self.ctx.val
+        )
+        if r == 0:
+            raise ZeroDivisionError(f"{self} is not invertible modulo {self.ctx.modulus()}")
+
+        return res
+
+    def __invert__(self):
+        return self.inverse()
+
+    def __pow__(self, e):
+        cdef bint check
+        cdef fmpz_mod res
+        res = fmpz_mod.__new__(fmpz_mod)
+        res.ctx = self.ctx
+
+        # Attempt to convert exponent to fmpz
+        e = any_as_fmpz(e)
+        if e is NotImplemented:
+            raise NotImplementedError
+
+        check = fmpz_mod_pow_fmpz(
+            res.val, self.val, (<fmpz>e).val, self.ctx.val
+        )
+
+        if check == 0:
+            raise ZeroDivisionError(f"{self} is not invertible modulo {self.ctx.modulus()}")
+
+        return res
diff --git a/src/flint/types/nmod.pyx b/src/flint/types/nmod.pyx
index c0a01bb8..53bee07b 100644
--- a/src/flint/types/nmod.pyx
+++ b/src/flint/types/nmod.pyx
@@ -47,9 +47,6 @@ cdef class nmod(flint_scalar):
 
     """
 
-    # cdef mp_limb_t val
-    # cdef nmod_t mod
-
     def __init__(self, val, mod):
         cdef mp_limb_t m
         m = mod