diff --git a/doc/source/fmpz_mod.rst 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
@@ -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
@@ -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
@@ -1,4 +1,4 @@
-**nmod** -- integers mod n
+**nmod** -- integers mod wordsize n
 ===============================================================================
 
 .. autoclass :: flint.nmod

diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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)
+