diff --git a/README.md b/README.md
index 0b2b8e40..5a5368a3 100644
--- a/README.md
+++ b/README.md
@@ -164,6 +164,7 @@ Contributors (0.7.0):
 - Giacomo Pope (GP)
 - Joris Roos (JR)
 - Edgar Costa (EC)
+- Frédéric Chapoton (FC)
 - Oscar Benjamin (OB)
 
 Highlights (0.7.0):
@@ -175,18 +176,34 @@ Highlights (0.7.0):
 - [gh-132](https://github.com/flintlib/python-flint/pull/132),
   [gh-164](https://github.com/flintlib/python-flint/pull/164),
   [gh-190](https://github.com/flintlib/python-flint/pull/190),
+  [gh-191](https://github.com/flintlib/python-flint/pull/191):
   [gh-192](https://github.com/flintlib/python-flint/pull/192):
+  [gh-216](https://github.com/flintlib/python-flint/pull/216):
+  [gh-225](https://github.com/flintlib/python-flint/pull/225):
   Add `fmpz_mpoly`, `fmpq_mpoly`, `nmod_poly` and `fmpz_mod_poly`
   types for multivariate polynomials with integer, rational or
-  integers mod n coefficients. (JM)
+  integers mod `n` coefficients. (JM)
 - [gh-142](https://github.com/flintlib/python-flint/pull/142)
   Add `acb_theta` module for the numerical evaluation of [theta
   functions](https://flintlib.org/doc/acb_theta.html) (only
   available for `Flint >= 3.1`). (EC)
+- [gh-218](https://github.com/flintlib/python-flint/pull/218)
+  An experimental interface for FLINT's generic rings has been
+  added. This provides access to many of FLINT's types that
+  are not yet wrapped by python-flint such as Gaussian integer,
+  number fields, qqbar, calcium, as well as both univariate and
+  multivariate polynomials and series over these rings (no
+  matrices yet though). (OB)
 - [gh-129](https://github.com/flintlib/python-flint/pull/129)
+  [gh-208](https://github.com/flintlib/python-flint/pull/208)
   Use meson/meson-python instead of setuptools as the build system
   for parallel builds and better detection of build and dependency
   requirements. (OB)
+- [gh-201](https://github.com/flintlib/python-flint/pull/201)
+  [gh-202](https://github.com/flintlib/python-flint/pull/202)
+  The documentation has been updated and is now at
+  [readthedocs](https://python-flint.readthedocs.io/en/latest/).
+  (OB)
 
 Compatibility break (0.7.0):
 
@@ -198,11 +215,23 @@ Compatibility break (0.7.0):
 
 Other changes (0.7.0):
 
+- [gh-215](https://github.com/flintlib/python-flint/pull/215)
+  [gh-219](https://github.com/flintlib/python-flint/pull/219)
+  The FLINT binding declarations are now fully generated
+  automatically from the FLINT docs. (OB)
+- [gh-203](https://github.com/flintlib/python-flint/pull/203)
+  [gh-204](https://github.com/flintlib/python-flint/pull/204)
+  [gh-205](https://github.com/flintlib/python-flint/pull/205)
+  [gh-206](https://github.com/flintlib/python-flint/pull/206)
+  [gh-207](https://github.com/flintlib/python-flint/pull/207)
+  [gh-211](https://github.com/flintlib/python-flint/pull/211)
+  [gh-212](https://github.com/flintlib/python-flint/pull/212)
+  Various linting fixes and codebase improvements (FC and GP).
 - [gh-189](https://github.com/flintlib/python-flint/pull/189)
   All scalar and poly types now have `sqrt`. All poly types now
   have `factor_squarefree` and `leading_coefficient` methods.
   Exception types raised in a number of places were changed to
-  `DomainError` for better consistency.
+  `DomainError` for better consistency. (OB)
 - [gh-196](https://github.com/flintlib/python-flint/pull/196)
   Supported Python versions are 3.10-3.13 (3.9 dropped). CI
   Testing added for 3.13 free-threaded CPython.
diff --git a/bin/all_rst_to_pxd.sh b/bin/all_rst_to_pxd.sh
index b6278216..86fb1123 100755
--- a/bin/all_rst_to_pxd.sh
+++ b/bin/all_rst_to_pxd.sh
@@ -15,15 +15,15 @@ modules=(
     "mpoly"
     # "thread_pool"
     # "machine_vectors"
-    # "gr"
-    # "gr_domains"
-    # "gr_generic"
-    # "gr_implementing"
-    # "gr_mat"
-    # "gr_mpoly"
-    # "gr_poly"
-    # "gr_special"
-    # "gr_vec"
+    "gr"
+    "gr_domains"
+    "gr_generic"
+    "gr_implementing"
+    "gr_mat"
+    "gr_mpoly"
+    "gr_poly"
+    "gr_special"
+    "gr_vec"
     "ulong_extras"
     "fmpz"
     # "fmpz_extras"
diff --git a/bin/rst_to_pxd.py b/bin/rst_to_pxd.py
index 5797b9b2..57f325fa 100755
--- a/bin/rst_to_pxd.py
+++ b/bin/rst_to_pxd.py
@@ -52,6 +52,9 @@
 c_types = set(["void", "char", "short", "long", "int", "float", "double"])
 type_modifers = re.compile(r"\*|(\bconst\b)|(\bunsigned\b)|(\bsigned\b)")
 import_dict = {}
+# gr_domains.rst lists functions that are in gr.h
+doc_to_header = {'flint/gr_domains': 'flint/gr'}
+
 
 def get_cython_struct_types(file):
     """
@@ -144,6 +147,9 @@ def get_all_types(function_list):
 
 
 def has_types(line, types):
+    #ts = set(t for t in get_parameter_types(line) if t in types)
+    #if ts:
+    #    print(ts, file=sys.stderr)
     return any(t in types for t in get_parameter_types(line))
 
 
@@ -171,7 +177,9 @@ def generate_pxd_file(h_name, opts):
     fill_import_dict(opts.flint_lib_dir)
     l=[]
     docdir = opts.arb_doc_dir
-    name = h_name 
+    name = h_name
+    h_name = doc_to_header.get(h_name, h_name)
+
     if name[:6] == "flint/":
         docdir = opts.flint_doc_dir
         name = name[6:]
diff --git a/doc/source/_gr.rst b/doc/source/_gr.rst
new file mode 100644
index 00000000..673f0918
--- /dev/null
+++ b/doc/source/_gr.rst
@@ -0,0 +1,151 @@
+**_gr** -- generic rings (unstable interface)
+===============================================================================
+
+.. note::
+   This module provides a preliminary experimental interface for FLINT's
+   generic rings. This is largely untested and the interface should be
+   considered unstable. In future the generic rings code will be integrated
+   with the rest of python-flint and this module may be removed. For now it
+   provides access to many FLINT types that are not wrapped yet in the rest of
+   python-flint.
+
+The generic rings module provides access to the generic rings interface in
+FLINT. Usage of this interface consists of creating a context object to
+represent a particular domain and then using that context object to create
+elements of that domain. For example to create a context for polynomials in two
+variables over the Gaussian integers :math:`\mathbb{Z}[i][x,y]` we would do::
+
+    >>> from flint.types._gr import gr_fmpzi_ctx, gr_gr_mpoly_ctx
+    >>> ctx = gr_gr_mpoly_ctx.new(gr_fmpzi_ctx, ["x", "y"])
+    >>> ctx.gens()
+    [x, y]
+    >>> ctx.gens_recursive()
+    [I, x, y]
+    >>> I, x, y = ctx.gens_recursive()
+    >>> p = (x + y + I)**2
+    >>> p
+    x^2 + 2*x*y + (2*I)*x + y^2 + (2*I)*y - 1
+
+Some domains such as ``gr_fmpzi_ctx`` are global and do not need to be created.
+Others such as ``gr_gr_mpoly_ctx`` are created using :meth:`gr_ctx.new`.
+
+.. autoclass :: flint.types._gr.gr_ctx
+  :members:
+  :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_scalar_ctx
+  :members:
+  :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_poly_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_mpoly_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr._gr_fmpz_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr._gr_fmpq_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr._gr_fmpzi_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr._gr_fexpr_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_nmod_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_fmpz_mod_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_fq_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_fq_nmod_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_fq_zech_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_nf_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_nf_fmpz_poly_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_real_qqbar_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_complex_qqbar_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_real_ca_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_complex_ca_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_real_algebraic_ca_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_complex_algebraic_ca_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_complex_extended_ca_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_real_float_arf_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_complex_float_acf_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_real_arb_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_complex_acb_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_gr_poly_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_gr_mpoly_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr_series_ctx
+   :members:
+   :undoc-members:
+
+.. autoclass :: flint.types._gr.gr
+   :members:
+   :inherited-members:
+   :undoc-members:
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 89345a87..193f047c 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -100,7 +100,7 @@
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'default'
+html_theme = 'furo'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
diff --git a/doc/source/index.rst b/doc/source/index.rst
index 887c68a7..d8bb93e8 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -96,3 +96,10 @@ Power series types
    arb_series.rst
    acb_series.rst
 
+Experimental generic rings interface
+....................................
+
+.. toctree::
+   :maxdepth: 2
+
+   _gr.rst
diff --git a/requirements-dev.txt b/requirements-dev.txt
index ecb5291f..1ee3ec44 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -8,3 +8,4 @@ coverage
 pytest-cov
 sphinx
 sphinx-rtd-theme
+furo
diff --git a/src/flint/flint_base/flint_base.pxd b/src/flint/flint_base/flint_base.pxd
index daf22a96..2f9ede7e 100644
--- a/src/flint/flint_base/flint_base.pxd
+++ b/src/flint/flint_base/flint_base.pxd
@@ -1,5 +1,8 @@
 from flint.flintlib.types.mpoly cimport ordering_t
 
+cdef class flint_ctx:
+    pass
+
 cdef class flint_elem:
     pass
 
diff --git a/src/flint/flintlib/functions/gr.pxd b/src/flint/flintlib/functions/gr.pxd
new file mode 100644
index 00000000..fc7feefc
--- /dev/null
+++ b/src/flint/flintlib/functions/gr.pxd
@@ -0,0 +1,187 @@
+from flint.flintlib.types.flint cimport flint_rand_t, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_ptr, gr_srcptr, gr_stream_t, gr_vec_t, truth_t
+
+# unknown type fexpr_t
+# unknown type fexpr_vec_t
+
+# .. macro:: GR_SUCCESS
+# .. macro:: GR_DOMAIN
+# .. macro:: GR_UNABLE
+# .. macro:: GR_TEST_FAIL
+# .. macro:: GR_MUST_SUCCEED(expr)
+# .. macro:: GR_TMP_INIT_VEC(vec, len, ctx)
+# .. macro:: GR_TMP_INIT(x1, ctx)
+# .. macro:: GR_TMP_CLEAR(x1, ctx)
+
+cdef extern from "flint/gr.h":
+    slong gr_ctx_sizeof_elem(gr_ctx_t ctx)
+    int gr_ctx_clear(gr_ctx_t ctx)
+    int gr_ctx_write(gr_stream_t out, gr_ctx_t ctx)
+    int gr_ctx_print(gr_ctx_t ctx)
+    int gr_ctx_println(gr_ctx_t ctx)
+    int gr_ctx_get_str(char ** s, gr_ctx_t ctx)
+    int gr_ctx_set_gen_name(gr_ctx_t ctx, const char * s)
+    int gr_ctx_set_gen_names(gr_ctx_t ctx, const char ** s)
+    void gr_init(gr_ptr res, gr_ctx_t ctx)
+    void gr_clear(gr_ptr res, gr_ctx_t ctx)
+    void gr_swap(gr_ptr x, gr_ptr y, gr_ctx_t ctx)
+    void gr_set_shallow(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    gr_ptr gr_heap_init(gr_ctx_t ctx)
+    void gr_heap_clear(gr_ptr x, gr_ctx_t ctx)
+    gr_ptr gr_heap_init_vec(slong len, gr_ctx_t ctx)
+    void gr_heap_clear_vec(gr_ptr x, slong len, gr_ctx_t ctx)
+    int gr_randtest(gr_ptr res, flint_rand_t state, gr_ctx_t ctx)
+    int gr_randtest_not_zero(gr_ptr res, flint_rand_t state, gr_ctx_t ctx)
+    int gr_randtest_small(gr_ptr res, flint_rand_t state, gr_ctx_t ctx)
+    int gr_write(gr_stream_t out, gr_srcptr x, gr_ctx_t ctx)
+    int gr_print(gr_srcptr x, gr_ctx_t ctx)
+    int gr_println(gr_srcptr x, gr_ctx_t ctx)
+    int gr_get_str(char ** s, gr_srcptr x, gr_ctx_t ctx)
+    int gr_set_str(gr_ptr res, const char * x, gr_ctx_t ctx)
+    int gr_write_n(gr_stream_t out, gr_srcptr x, slong n, gr_ctx_t ctx)
+    int gr_get_str_n(char ** s, gr_srcptr x, slong n, gr_ctx_t ctx)
+    int gr_set(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_set_other(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int gr_set_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_set_si(gr_ptr res, slong x, gr_ctx_t ctx)
+    int gr_set_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_set_fmpq(gr_ptr res, const fmpq_t x, gr_ctx_t ctx)
+    int gr_set_d(gr_ptr res, double x, gr_ctx_t ctx)
+    int gr_get_si(slong * res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_get_ui(ulong * res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_get_fmpz(fmpz_t res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_get_fmpq(fmpq_t res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_get_d(double * res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_set_fmpz_2exp_fmpz(gr_ptr res, const fmpz_t a, const fmpz_t b, gr_ctx_t ctx)
+    int gr_get_fmpz_2exp_fmpz(fmpz_t res1, fmpz_t res2, gr_srcptr x, gr_ctx_t ctx)
+    int gr_set_fmpz_10exp_fmpz(gr_ptr res, const fmpz_t a, const fmpz_t b, gr_ctx_t ctx)
+    # int gr_get_fexpr(fexpr_t res, gr_srcptr x, gr_ctx_t ctx)
+    # int gr_get_fexpr_serialize(fexpr_t res, gr_srcptr x, gr_ctx_t ctx)
+    # int gr_set_fexpr(gr_ptr res, fexpr_vec_t inputs, gr_vec_t outputs, const fexpr_t x, gr_ctx_t ctx)
+    int gr_zero(gr_ptr res, gr_ctx_t ctx)
+    int gr_one(gr_ptr res, gr_ctx_t ctx)
+    int gr_neg_one(gr_ptr res, gr_ctx_t ctx)
+    int gr_gen(gr_ptr res, gr_ctx_t ctx)
+    int gr_gens(gr_vec_t res, gr_ctx_t ctx)
+    int gr_gens_recursive(gr_vec_t res, gr_ctx_t ctx)
+    truth_t gr_is_zero(gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_is_one(gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_is_neg_one(gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_equal(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_is_integer(gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_is_rational(gr_srcptr x, gr_ctx_t ctx)
+    int gr_neg(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_add(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_add_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_add_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_add_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_add_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_add_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_other_add(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_sub(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_sub_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_sub_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_sub_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_sub_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_sub_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_other_sub(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_mul(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_mul_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_mul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_mul_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_mul_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_mul_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_other_mul(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_addmul(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_addmul_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_addmul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_addmul_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_addmul_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_addmul_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_submul(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_submul_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_submul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_submul_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_submul_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_submul_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_mul_two(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sqr(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_mul_2exp_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_mul_2exp_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    truth_t gr_is_invertible(gr_srcptr x, gr_ctx_t ctx)
+    int gr_inv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_div_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_div_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_div_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_div_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_div_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_other_div(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_div_nonunique(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_divides(gr_srcptr d, gr_srcptr x, gr_ctx_t ctx)
+    int gr_divexact(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_divexact_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_divexact_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_divexact_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_divexact_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_other_divexact(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_euclidean_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_euclidean_rem(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_euclidean_divrem(gr_ptr res1, gr_ptr res2, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_pow(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_pow_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_pow_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_pow_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_pow_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_pow_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_other_pow(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_is_square(gr_srcptr x, gr_ctx_t ctx)
+    int gr_sqrt(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_rsqrt(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_gcd(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_lcm(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_factor(gr_ptr c, gr_vec_t factors, gr_vec_t exponents, gr_srcptr x, int flags, gr_ctx_t ctx)
+    int gr_numerator(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_denominator(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_floor(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_ceil(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_trunc(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_nint(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_abs(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_i(gr_ptr res, gr_ctx_t ctx)
+    int gr_conj(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_re(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_im(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sgn(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_csgn(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_arg(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_pos_inf(gr_ptr res, gr_ctx_t ctx)
+    int gr_neg_inf(gr_ptr res, gr_ctx_t ctx)
+    int gr_uinf(gr_ptr res, gr_ctx_t ctx)
+    int gr_undefined(gr_ptr res, gr_ctx_t ctx)
+    int gr_unknown(gr_ptr res, gr_ctx_t ctx)
+    int gr_cmp(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_cmp_other(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_cmpabs(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_cmpabs_other(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    truth_t gr_le(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_lt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_ge(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_gt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_abs_le(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_abs_lt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_abs_ge(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    truth_t gr_abs_gt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_min(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_max(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_set_interval_mid_rad(gr_ptr res, gr_srcptr m, gr_srcptr r, gr_ctx_t ctx)
+    int gr_ctx_fq_prime(fmpz_t p, gr_ctx_t ctx)
+    int gr_ctx_fq_degree(slong * deg, gr_ctx_t ctx)
+    int gr_ctx_fq_order(fmpz_t q, gr_ctx_t ctx)
+    int gr_fq_frobenius(gr_ptr res, gr_srcptr x, slong e, gr_ctx_t ctx)
+    int gr_fq_multiplicative_order(fmpz_t res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fq_norm(fmpz_t res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fq_trace(fmpz_t res, gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_fq_is_primitive(gr_srcptr x, gr_ctx_t ctx)
+    int gr_fq_pth_root(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
diff --git a/src/flint/flintlib/functions/gr_domains.pxd b/src/flint/flintlib/functions/gr_domains.pxd
new file mode 100644
index 00000000..3988ec74
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_domains.pxd
@@ -0,0 +1,73 @@
+from flint.flintlib.types.flint cimport flint_rand_t, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_poly_t
+from flint.flintlib.types.fmpz cimport fmpz_poly_t
+from flint.flintlib.types.gr cimport gr_ctx_t, truth_t
+from flint.flintlib.types.mpoly cimport ordering_t
+
+
+
+cdef extern from "flint/gr.h":
+    int gr_ctx_cmp_coercion(gr_ctx_t ctx1, gr_ctx_t ctx2)
+    truth_t gr_ctx_is_finite(gr_ctx_t ctx)
+    truth_t gr_ctx_is_multiplicative_group(gr_ctx_t ctx)
+    truth_t gr_ctx_is_ring(gr_ctx_t ctx)
+    truth_t gr_ctx_is_commutative_ring(gr_ctx_t ctx)
+    truth_t gr_ctx_is_integral_domain(gr_ctx_t ctx)
+    truth_t gr_ctx_is_unique_factorization_domain(gr_ctx_t ctx)
+    truth_t gr_ctx_is_field(gr_ctx_t ctx)
+    truth_t gr_ctx_is_algebraically_closed(gr_ctx_t ctx)
+    truth_t gr_ctx_is_finite_characteristic(gr_ctx_t ctx)
+    truth_t gr_ctx_is_ordered_ring(gr_ctx_t ctx)
+    truth_t gr_ctx_is_zero_ring(gr_ctx_t ctx)
+    truth_t gr_ctx_is_exact(gr_ctx_t ctx)
+    truth_t gr_ctx_is_canonical(gr_ctx_t ctx)
+    truth_t gr_ctx_has_real_prec(gr_ctx_t ctx)
+    int gr_ctx_set_real_prec(gr_ctx_t ctx, slong prec)
+    int gr_ctx_get_real_prec(slong * prec, gr_ctx_t ctx)
+    void gr_ctx_init_perm(gr_ctx_t ctx, ulong n)
+    void gr_ctx_init_psl2z(gr_ctx_t ctx)
+    int gr_ctx_init_dirichlet_group(gr_ctx_t ctx, ulong q)
+    void gr_ctx_init_random(gr_ctx_t ctx, flint_rand_t state)
+    void gr_ctx_init_fmpz(gr_ctx_t ctx)
+    void gr_ctx_init_fmpq(gr_ctx_t ctx)
+    void gr_ctx_init_fmpzi(gr_ctx_t ctx)
+    int gr_ctx_set_is_field(gr_ctx_t ctx, truth_t is_field)
+    void gr_ctx_init_nmod(gr_ctx_t ctx, ulong n)
+    void gr_ctx_init_nmod8(gr_ctx_t ctx, unsigned char n)
+    void gr_ctx_init_nmod32(gr_ctx_t ctx, unsigned int n)
+    void gr_ctx_init_fmpz_mod(gr_ctx_t ctx, const fmpz_t n)
+    void gr_ctx_init_fq(gr_ctx_t ctx, const fmpz_t p, slong d, const char * var)
+    void gr_ctx_init_fq_nmod(gr_ctx_t ctx, ulong p, slong d, const char * var)
+    void gr_ctx_init_fq_zech(gr_ctx_t ctx, ulong p, slong d, const char * var)
+    void gr_ctx_init_nf(gr_ctx_t ctx, const fmpq_poly_t poly)
+    void gr_ctx_init_nf_fmpz_poly(gr_ctx_t ctx, const fmpz_poly_t poly)
+    void gr_ctx_init_real_qqbar(gr_ctx_t ctx)
+    void gr_ctx_init_complex_qqbar(gr_ctx_t ctx)
+    void _gr_ctx_qqbar_set_limits(gr_ctx_t ctx, slong deg_limit, slong bits_limit)
+    void gr_ctx_init_real_arb(gr_ctx_t ctx, slong prec)
+    void gr_ctx_init_complex_acb(gr_ctx_t ctx, slong prec)
+    void gr_ctx_arb_set_prec(gr_ctx_t ctx, slong prec)
+    slong gr_ctx_arb_get_prec(gr_ctx_t ctx)
+    void gr_ctx_init_real_ca(gr_ctx_t ctx)
+    void gr_ctx_init_complex_ca(gr_ctx_t ctx)
+    void gr_ctx_init_real_algebraic_ca(gr_ctx_t ctx)
+    void gr_ctx_init_complex_algebraic_ca(gr_ctx_t ctx)
+    void gr_ctx_ca_set_option(gr_ctx_t ctx, slong option, slong value)
+    slong gr_ctx_ca_get_option(gr_ctx_t ctx, slong option)
+    void gr_ctx_init_complex_extended_ca(gr_ctx_t ctx)
+    void gr_ctx_init_real_float_arf(gr_ctx_t ctx, slong prec)
+    void gr_ctx_init_complex_float_acf(gr_ctx_t ctx, slong prec)
+    void gr_ctx_init_vector_gr_vec(gr_ctx_t ctx, gr_ctx_t base_type)
+    void gr_ctx_init_vector_space_gr_vec(gr_ctx_t ctx, gr_ctx_t base_type, slong n)
+    void gr_ctx_init_matrix_domain(gr_ctx_t ctx, gr_ctx_t base_ring)
+    void gr_ctx_init_matrix_space(gr_ctx_t ctx, gr_ctx_t base_ring, slong n, slong m)
+    void gr_ctx_init_matrix_ring(gr_ctx_t ctx, gr_ctx_t base_ring, slong n)
+    void gr_ctx_init_fmpz_poly(gr_ctx_t ctx)
+    void gr_ctx_init_fmpq_poly(gr_ctx_t ctx)
+    void gr_ctx_init_gr_poly(gr_ctx_t ctx, gr_ctx_t base_ring)
+    void gr_ctx_init_fmpz_mpoly(gr_ctx_t ctx, slong nvars, const ordering_t ord)
+    void gr_ctx_init_gr_mpoly(gr_ctx_t ctx, gr_ctx_t base_ring, slong nvars, const ordering_t ord)
+    void gr_ctx_init_series_mod_gr_poly(gr_ctx_t ctx, gr_ctx_t base_ring, slong n)
+    void gr_ctx_init_gr_series(gr_ctx_t ctx, gr_ctx_t base_ring, slong prec)
+    void gr_ctx_init_fmpz_mpoly_q(gr_ctx_t ctx, slong nvars, const ordering_t ord)
+    void gr_ctx_init_fexpr(gr_ctx_t ctx)
diff --git a/src/flint/flintlib/functions/gr_generic.pxd b/src/flint/flintlib/functions/gr_generic.pxd
new file mode 100644
index 00000000..9f3bfa58
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_generic.pxd
@@ -0,0 +1,214 @@
+from flint.flintlib.types.flint cimport flint_rand_t, fmpz_struct, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_t
+from flint.flintlib.types.fmpz cimport fmpz_mpoly_ctx_t, fmpz_mpoly_t, fmpz_poly_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_ptr, gr_srcptr, gr_stream_t, truth_t
+
+
+# .. macro :: GR_PARSE_BALANCE_ADDITIONS
+# .. macro :: GR_PARSE_RING_EXPONENTS
+
+cdef extern from "flint/gr_generic.h":
+    void gr_generic_init(void)
+    void gr_generic_clear(void)
+    void gr_generic_swap(void)
+    void gr_generic_randtest(void)
+    void gr_generic_write(void)
+    void gr_generic_zero(void)
+    void gr_generic_one(void)
+    void gr_generic_equal(void)
+    void gr_generic_set(void)
+    void gr_generic_set_si(void)
+    void gr_generic_set_ui(void)
+    void gr_generic_set_fmpz(void)
+    void gr_generic_neg(void)
+    void gr_generic_add(void)
+    void gr_generic_sub(void)
+    void gr_generic_mul(void)
+    int gr_generic_ctx_clear(gr_ctx_t ctx)
+    void gr_generic_set_shallow(gr_ptr res, gr_srcptr x, const gr_ctx_t ctx)
+    int gr_generic_write_n(gr_stream_t out, gr_srcptr x, slong n, gr_ctx_t ctx)
+    int gr_generic_randtest_not_zero(gr_ptr x, flint_rand_t state, gr_ctx_t ctx)
+    int gr_generic_randtest_small(gr_ptr x, flint_rand_t state, gr_ctx_t ctx)
+    truth_t gr_generic_is_zero(gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_generic_is_one(gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_generic_is_neg_one(gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_neg_one(gr_ptr res, gr_ctx_t ctx)
+    int gr_generic_set_other(gr_ptr res, gr_srcptr x, gr_ctx_t xctx, gr_ctx_t ctx)
+    int gr_generic_set_fmpq(gr_ptr res, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_set_str_expr(gr_ptr res, const char * s, int flags, gr_ctx_t ctx)
+    int gr_generic_set_str(gr_ptr res, const char * s, gr_ctx_t ctx)
+    int gr_generic_set_str_balance_additions(gr_ptr res, const char * s, gr_ctx_t ctx)
+    int gr_generic_set_str_ring_exponents(gr_ptr res, const char * s, gr_ctx_t ctx)
+    int gr_generic_add_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_add_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_generic_add_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_add_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_add_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_other_add(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_sub_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_generic_sub_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_sub_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_sub_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_sub_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_other_sub(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_mul_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_mul_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_generic_mul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_mul_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_mul_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_other_mul(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_addmul(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_addmul_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_generic_addmul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_addmul_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_addmul_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_addmul_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_submul(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_submul_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_generic_submul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_submul_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_submul_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_submul_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_mul_two(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_sqr(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_mul_2exp_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_mul_2exp_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_set_fmpz_2exp_fmpz(gr_ptr res, const fmpz_t x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_get_fmpz_2exp_fmpz(fmpz_t res1, fmpz_t res2, gr_ptr x, gr_ctx_t ctx)
+    int gr_generic_inv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    truth_t gr_generic_is_invertible(gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_div_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
+    int gr_generic_div_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_generic_div_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
+    int gr_generic_div_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_div_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_other_div(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_divexact(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_pow_fmpz_sliding(gr_ptr f, gr_srcptr g, const fmpz_t pow, gr_ctx_t ctx)
+    int gr_generic_pow_ui_sliding(gr_ptr f, gr_srcptr g, ulong pow, gr_ctx_t ctx)
+    int gr_generic_pow_fmpz_binexp(gr_ptr res, gr_srcptr x, const fmpz_t exp, gr_ctx_t ctx)
+    int gr_generic_pow_ui_binexp(gr_ptr res, gr_srcptr x, ulong e, gr_ctx_t ctx)
+    int gr_generic_pow_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t e, gr_ctx_t ctx)
+    int gr_generic_pow_si(gr_ptr res, gr_srcptr x, slong e, gr_ctx_t ctx)
+    int gr_generic_pow_ui(gr_ptr res, gr_srcptr x, ulong e, gr_ctx_t ctx)
+    int gr_generic_pow_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
+    int gr_generic_pow_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_other_pow(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)
+    int _gr_fmpz_poly_evaluate_horner(gr_ptr res, const fmpz_struct * f, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fmpz_poly_evaluate_horner(gr_ptr res, const fmpz_poly_t f, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_fmpz_poly_evaluate_rectangular(gr_ptr res, const fmpz_struct * f, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fmpz_poly_evaluate_rectangular(gr_ptr res, const fmpz_poly_t f, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_fmpz_poly_evaluate(gr_ptr res, const fmpz_struct * f, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fmpz_poly_evaluate(gr_ptr res, const fmpz_poly_t f, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fmpz_mpoly_evaluate_iter(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
+    int gr_fmpz_mpoly_evaluate_horner(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
+    int gr_fmpz_mpoly_evaluate(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
+    truth_t gr_generic_is_square(gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_sqrt(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_rsqrt(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_numerator(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_denominator(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_generic_cmp(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_cmpabs(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_generic_cmp_other(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_cmpabs_other(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
+    int gr_generic_bernoulli_ui(gr_ptr res, ulong n, gr_ctx_t ctx)
+    int gr_generic_bernoulli_fmpz(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_generic_bernoulli_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_generic_eulernum_ui(gr_ptr res, ulong n, gr_ctx_t ctx)
+    int gr_generic_eulernum_fmpz(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_generic_eulernum_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_generic_stirling_s1u_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_generic_stirling_s1_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_generic_stirling_s2_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_generic_stirling_s1u_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    int gr_generic_stirling_s1_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    int gr_generic_stirling_s2_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    void gr_generic_vec_init(gr_ptr vec, slong len, gr_ctx_t ctx)
+    void gr_generic_vec_clear(gr_ptr vec, slong len, gr_ctx_t ctx)
+    void gr_generic_vec_swap(gr_ptr vec1, gr_ptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_zero(gr_ptr vec, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_set(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_neg(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_normalise(slong * res, gr_srcptr vec, slong len, gr_ctx_t ctx)
+    slong gr_generic_vec_normalise_weak(gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar_2exp_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_scalar_addmul(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_scalar_submul(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_scalar_addmul_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_scalar_submul_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    truth_t gr_generic_vec_equal(gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_is_zero(gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_dot(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_dot_rev(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_dot_ui(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, const ulong * vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_dot_si(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, const slong * vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_dot_fmpz(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, const fmpz_struct * vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_set_powers(gr_ptr res, gr_srcptr x, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_reciprocals(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_add(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_sub(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_mul(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_div(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_divexact(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_pow(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_add_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_sub_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_div_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_pow_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_generic_vec_add_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_sub_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_div_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_pow_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int gr_generic_vec_add_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int gr_generic_vec_sub_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int gr_generic_vec_div_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int gr_generic_vec_pow_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int gr_generic_vec_add_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int gr_generic_vec_sub_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int gr_generic_vec_div_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int gr_generic_vec_pow_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int gr_generic_vec_add_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int gr_generic_vec_sub_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int gr_generic_vec_div_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int gr_generic_vec_pow_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int gr_generic_scalar_add_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_sub_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_mul_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_div_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_divexact_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_pow_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_add_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_sub_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_mul_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_div_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_pow_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int gr_generic_other_add_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int gr_generic_other_sub_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int gr_generic_other_mul_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int gr_generic_other_div_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int gr_generic_other_divexact_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int gr_generic_other_pow_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int gr_generic_vec_add_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int gr_generic_vec_sub_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int gr_generic_vec_mul_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int gr_generic_vec_div_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int gr_generic_vec_divexact_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int gr_generic_vec_pow_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int gr_generic_scalar_other_add_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_other_sub_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_other_mul_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_other_div_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_other_divexact_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int gr_generic_scalar_other_pow_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
diff --git a/src/flint/flintlib/functions/gr_implementing.pxd b/src/flint/flintlib/functions/gr_implementing.pxd
new file mode 100644
index 00000000..c97fd3aa
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_implementing.pxd
@@ -0,0 +1,14 @@
+from flint.flintlib.types.flint cimport slong
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_funcptr, truth_t
+
+# unknown type gr_method_tab_input
+
+
+cdef extern from "flint/gr_implementing.h":
+    # void gr_method_tab_init(gr_funcptr * methods, gr_method_tab_input * tab)
+    int gr_not_implemented(void)
+    int gr_not_in_domain(void)
+    truth_t gr_generic_ctx_predicate(gr_ctx_t ctx)
+    truth_t gr_generic_ctx_predicate_true(gr_ctx_t ctx)
+    truth_t gr_generic_ctx_predicate_false(gr_ctx_t ctx)
+    void gr_test_ring(gr_ctx_t R, slong iters, int test_flags)
diff --git a/src/flint/flintlib/functions/gr_mat.pxd b/src/flint/flintlib/functions/gr_mat.pxd
new file mode 100644
index 00000000..539b1483
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_mat.pxd
@@ -0,0 +1,184 @@
+from flint.flintlib.types.flint cimport flint_rand_t, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_mat_t, fmpq_t
+from flint.flintlib.types.fmpz cimport fmpz_mat_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_mat_t, gr_poly_t, gr_ptr, gr_srcptr, gr_stream_t, gr_vec_t, truth_t
+
+# unknown type gr_method_binary_op
+# unknown type gr_method_binary_predicate
+# unknown type gr_method_mat_binary_op
+# unknown type gr_method_mat_binary_op_with_flag
+# unknown type gr_method_mat_lu_op
+# unknown type gr_method_mat_unary_op_get_scalar
+# unknown type gr_method_unary_op
+# unknown type gr_method_unary_predicate
+
+# .. macro:: GR_MAT_ENTRY(mat, i, j, sz)
+# .. macro:: gr_mat_nrows(mat, ctx)
+# .. macro:: gr_mat_ncols(mat, ctx)
+
+cdef extern from "flint/gr_mat.h":
+    gr_ptr gr_mat_entry_ptr(gr_mat_t mat, slong i, slong j, gr_ctx_t ctx)
+    void gr_mat_init(gr_mat_t mat, slong rows, slong cols, gr_ctx_t ctx)
+    int gr_mat_init_set(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    void gr_mat_clear(gr_mat_t mat, gr_ctx_t ctx)
+    void gr_mat_swap(gr_mat_t mat1, gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_swap_entrywise(gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    void gr_mat_window_init(gr_mat_t window, const gr_mat_t mat, slong r1, slong c1, slong r2, slong c2, gr_ctx_t ctx)
+    void gr_mat_window_clear(gr_mat_t window, gr_ctx_t ctx)
+    int gr_mat_write(gr_stream_t out, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_print(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_equal(const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    truth_t gr_mat_is_zero(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_one(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_neg_one(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_scalar(const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_zero(gr_mat_t res, gr_ctx_t ctx)
+    int gr_mat_one(gr_mat_t res, gr_ctx_t ctx)
+    int gr_mat_set(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_set_fmpz_mat(gr_mat_t res, const fmpz_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_set_fmpq_mat(gr_mat_t res, const fmpq_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_set_gr_mat_other(gr_mat_t res, const gr_mat_t mat, gr_ctx_t mat_ctx, gr_ctx_t ctx)
+    int gr_mat_set_scalar(gr_mat_t res, gr_srcptr c, gr_ctx_t ctx)
+    int gr_mat_set_ui(gr_mat_t res, ulong c, gr_ctx_t ctx)
+    int gr_mat_set_si(gr_mat_t res, slong c, gr_ctx_t ctx)
+    int gr_mat_set_fmpz(gr_mat_t res, const fmpz_t c, gr_ctx_t ctx)
+    int gr_mat_set_fmpq(gr_mat_t res, const fmpq_t c, gr_ctx_t ctx)
+    int gr_mat_concat_horizontal(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_concat_vertical(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_transpose(gr_mat_t B, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_swap_rows(gr_mat_t mat, slong * perm, slong r, slong s, gr_ctx_t ctx)
+    int gr_mat_swap_cols(gr_mat_t mat, slong * perm, slong r, slong s, gr_ctx_t ctx)
+    int gr_mat_invert_rows(gr_mat_t mat, slong * perm, gr_ctx_t ctx)
+    int gr_mat_invert_cols(gr_mat_t mat, slong * perm, gr_ctx_t ctx)
+    truth_t gr_mat_is_empty(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_square(const gr_mat_t mat, gr_ctx_t ctx)
+    # int gr_mat_entrywise_unary_op(gr_mat_t res, gr_method_unary_op f, const gr_mat_t mat, gr_ctx_t ctx)
+    # int gr_mat_entrywise_binary_op(gr_mat_t res, gr_method_binary_op f, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    # int gr_mat_entrywise_binary_op_scalar(gr_mat_t res, gr_method_binary_op f, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    # truth_t gr_mat_entrywise_unary_predicate_all(gr_method_unary_predicate f, const gr_mat_t mat, gr_ctx_t ctx)
+    # truth_t gr_mat_entrywise_unary_predicate_any(gr_method_unary_predicate f, const gr_mat_t mat, gr_ctx_t ctx)
+    # truth_t gr_mat_entrywise_binary_predicate_all(gr_method_binary_predicate f, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_norm_max(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_norm_1(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_norm_inf(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_norm_frobenius(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_neg(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_add(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_sub(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_mul_classical(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_mul_strassen(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_mul_generic(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_mul(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+    int gr_mat_sqr(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_add_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    int gr_mat_sub_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    int gr_mat_mul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    int gr_mat_addmul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    int gr_mat_submul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    int gr_mat_div_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_mat_gr_poly_evaluate(gr_mat_t res, gr_srcptr poly, slong len, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_gr_poly_evaluate(gr_mat_t res, const gr_poly_t poly, const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_upper_triangular(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_lower_triangular(const gr_mat_t mat, gr_ctx_t ctx)
+    truth_t gr_mat_is_diagonal(const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_mul_diag(gr_mat_t res, const gr_mat_t A, const gr_vec_t D, gr_ctx_t ctx)
+    int gr_mat_diag_mul(gr_mat_t res, const gr_vec_t D, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_find_nonzero_pivot_large_abs(slong * pivot_row, gr_mat_t mat, slong start_row, slong end_row, slong column, gr_ctx_t ctx)
+    int gr_mat_find_nonzero_pivot_generic(slong * pivot_row, gr_mat_t mat, slong start_row, slong end_row, slong column, gr_ctx_t ctx)
+    int gr_mat_find_nonzero_pivot(slong * pivot_row, gr_mat_t mat, slong start_row, slong end_row, slong column, gr_ctx_t ctx)
+    int gr_mat_lu_classical(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
+    int gr_mat_lu_recursive(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
+    int gr_mat_lu_generic(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
+    int gr_mat_lu(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
+    int gr_mat_fflu(slong * rank, slong * P, gr_mat_t LU, gr_ptr den, const gr_mat_t A, int rank_check, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_tril_classical(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_tril_recursive(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_tril_generic(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_tril(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_triu_classical(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_triu_recursive(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_triu_generic(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_triu(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_fflu(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_lu(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_fflu_precomp(gr_mat_t X, const slong * perm, const gr_mat_t LU, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_lu_precomp(gr_mat_t X, const slong * perm, const gr_mat_t LU, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_den_fflu(gr_mat_t X, gr_ptr den, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_nonsingular_solve_den(gr_mat_t X, gr_ptr den, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_solve_field(gr_mat_t X, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+    int gr_mat_det_fflu(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_det_berkowitz(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_det_lu(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_det_cofactor(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_det_generic_field(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_det_generic_integral_domain(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_det_generic(gr_ptr res, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_det(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_trace(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_rank_fflu(slong * rank, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_rank_lu(slong * rank, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_rank(slong * rank, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_rref_lu(slong * rank, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_rref_fflu(slong * rank, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_rref(slong * rank, gr_mat_t R, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_rref_den_fflu(slong * rank, gr_mat_t R, gr_ptr den, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_rref_den(slong * rank, gr_mat_t R, gr_ptr den, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_nullspace(gr_mat_t X, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_inv(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_adjugate_charpoly(gr_mat_t adj, gr_ptr det, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_adjugate_cofactor(gr_mat_t adj, gr_ptr det, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_adjugate(gr_mat_t adj, gr_ptr det, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_berkowitz(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_berkowitz(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_danilevsky_inplace(gr_ptr res, gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_danilevsky(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_danilevsky(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_gauss(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_gauss(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_householder(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_householder(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_faddeev(gr_ptr res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_faddeev(gr_poly_t res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_faddeev_bsgs(gr_ptr res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_faddeev_bsgs(gr_poly_t res, gr_mat_t adj, const gr_mat_t mat, gr_ctx_t ctx)
+    int _gr_mat_charpoly_from_hessenberg(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_charpoly_from_hessenberg(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_minpoly_field(gr_poly_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_apply_row_similarity(gr_mat_t M, slong r, gr_ptr d, gr_ctx_t ctx)
+    int gr_mat_eigenvalues(gr_vec_t lambda_, gr_vec_t mult, const gr_mat_t mat, int flags, gr_ctx_t ctx)
+    int gr_mat_eigenvalues_other(gr_vec_t lambda_, gr_vec_t mult, const gr_mat_t mat, gr_ctx_t mat_ctx, int flags, gr_ctx_t ctx)
+    int gr_mat_diagonalization_precomp(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_t A, const gr_vec_t eigenvalues, const gr_vec_t mult, gr_ctx_t ctx)
+    int gr_mat_diagonalization_generic(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_t A, int flags, gr_ctx_t ctx)
+    int gr_mat_diagonalization(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_t A, int flags, gr_ctx_t ctx)
+    int gr_mat_set_jordan_blocks(gr_mat_t mat, const gr_vec_t lambda_, slong num_blocks, slong * block_lambda, slong * block_size, gr_ctx_t ctx)
+    int gr_mat_jordan_blocks(gr_vec_t lambda_, slong * num_blocks, slong * block_lambda, slong * block_size, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_jordan_transformation(gr_mat_t mat, const gr_vec_t lambda_, slong num_blocks, slong * block_lambda, slong * block_size, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_jordan_form(gr_mat_t J, gr_mat_t P, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_exp_jordan(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_exp(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_log_jordan(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
+    int gr_mat_log(gr_mat_t res, const gr_mat_t A, gr_ctx_t ctx)
+    truth_t gr_mat_is_hessenberg(const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_hessenberg_gauss(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_hessenberg_householder(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_hessenberg(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
+    int gr_mat_randtest(gr_mat_t res, flint_rand_t state, gr_ctx_t ctx)
+    int gr_mat_randops(gr_mat_t mat, flint_rand_t state, slong count, gr_ctx_t ctx)
+    int gr_mat_randpermdiag(int * parity, gr_mat_t mat, flint_rand_t state, gr_ptr diag, slong n, gr_ctx_t ctx)
+    int gr_mat_randrank(gr_mat_t mat, flint_rand_t state, slong rank, gr_ctx_t ctx)
+    int gr_mat_ones(gr_mat_t res, gr_ctx_t ctx)
+    int gr_mat_pascal(gr_mat_t res, int triangular, gr_ctx_t ctx)
+    int gr_mat_stirling(gr_mat_t res, int kind, gr_ctx_t ctx)
+    int gr_mat_hilbert(gr_mat_t res, gr_ctx_t ctx)
+    int gr_mat_hadamard(gr_mat_t res, gr_ctx_t ctx)
+    int gr_mat_reduce_row(slong * column, gr_mat_t A, slong * P, slong * L, slong m, gr_ctx_t ctx)
+    # void gr_mat_test_mul(gr_method_mat_binary_op mul_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void gr_mat_test_lu(gr_method_mat_lu_op lu_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void gr_mat_test_det(gr_method_mat_unary_op_get_scalar det_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void gr_mat_test_nonsingular_solve_tril(gr_method_mat_binary_op_with_flag solve_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void gr_mat_test_nonsingular_solve_triu(gr_method_mat_binary_op_with_flag solve_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void gr_mat_test_approx_mul_max_norm(gr_method_mat_binary_op mul_impl, gr_srcptr rel_tol, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void gr_mat_test_approx_mul_pos_entrywise_accurate(gr_method_mat_binary_op mul_impl, gr_srcptr rel_tol, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
diff --git a/src/flint/flintlib/functions/gr_mpoly.pxd b/src/flint/flintlib/functions/gr_mpoly.pxd
new file mode 100644
index 00000000..03724074
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_mpoly.pxd
@@ -0,0 +1,59 @@
+from flint.flintlib.types.flint cimport flint_bitcnt_t, flint_rand_t, fmpz_struct, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_mpoly_t, gr_ptr, gr_srcptr, gr_stream_t, truth_t
+from flint.flintlib.types.mpoly cimport mpoly_ctx_t
+
+
+
+cdef extern from "flint/gr_mpoly.h":
+    void gr_mpoly_init(gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_init3(gr_mpoly_t A, slong alloc, flint_bitcnt_t bits, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_init2(gr_mpoly_t A, slong alloc, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_clear(gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_swap(gr_mpoly_t A, gr_mpoly_t B, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set(gr_mpoly_t A, const gr_mpoly_t B, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_zero(gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    truth_t gr_mpoly_is_zero(const gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_gen(gr_mpoly_t A, slong var, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    truth_t gr_mpoly_is_gen(const gr_mpoly_t A, slong var, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    truth_t gr_mpoly_equal(const gr_mpoly_t A, const gr_mpoly_t B, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_randtest_bits(gr_mpoly_t A, flint_rand_t state, slong length, flint_bitcnt_t exp_bits, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_write_pretty(gr_stream_t out, const gr_mpoly_t A, const char ** x, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_print_pretty(const gr_mpoly_t A, const char ** x, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_get_coeff_scalar_fmpz(gr_ptr c, const gr_mpoly_t A, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_get_coeff_scalar_ui(gr_ptr c, const gr_mpoly_t A, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_scalar_fmpz(gr_mpoly_t A, gr_srcptr c, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_ui_fmpz(gr_mpoly_t A, ulong c, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_si_fmpz(gr_mpoly_t A, slong c, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_fmpz_fmpz(gr_mpoly_t A, const fmpz_t c, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_fmpq_fmpz(gr_mpoly_t A, const fmpq_t c, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_scalar_ui(gr_mpoly_t poly, gr_srcptr c, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_ui_ui(gr_mpoly_t A, ulong c, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_si_ui(gr_mpoly_t A, slong c, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_fmpz_ui(gr_mpoly_t A, const fmpz_t c, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_set_coeff_fmpq_ui(gr_mpoly_t A, const fmpq_t c, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_neg(gr_mpoly_t A, const gr_mpoly_t B, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_add(gr_mpoly_t A, const gr_mpoly_t B, const gr_mpoly_t C, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_sub(gr_mpoly_t A, const gr_mpoly_t B, const gr_mpoly_t C, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul(gr_mpoly_t A, const gr_mpoly_t B, const gr_mpoly_t C, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_johnson(gr_mpoly_t A, const gr_mpoly_t B, const gr_mpoly_t C, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_monomial(gr_mpoly_t A, const gr_mpoly_t B, const gr_mpoly_t C, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_scalar(gr_mpoly_t A, const gr_mpoly_t B, gr_srcptr c, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_si(gr_mpoly_t A, const gr_mpoly_t B, slong c, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_ui(gr_mpoly_t A, const gr_mpoly_t B, ulong c, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_fmpz(gr_mpoly_t A, const gr_mpoly_t B, const fmpz_t c, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_mul_fmpq(gr_mpoly_t A, const gr_mpoly_t B, const fmpq_t c, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void _gr_mpoly_fit_length(gr_ptr * coeffs, slong * coeffs_alloc, ulong ** exps, slong * exps_alloc, slong N, slong length, gr_ctx_t cctx)
+    void gr_mpoly_fit_length(gr_mpoly_t A, slong len, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_fit_bits(gr_mpoly_t A, flint_bitcnt_t bits, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_fit_length_fit_bits(gr_mpoly_t A, slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_fit_length_reset_bits(gr_mpoly_t A, slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void _gr_mpoly_set_length(gr_mpoly_t A, slong newlen, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void _gr_mpoly_push_exp_ui(gr_mpoly_t A, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_push_term_scalar_ui(gr_mpoly_t A, gr_srcptr c, const ulong * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void _gr_mpoly_push_exp_fmpz(gr_mpoly_t A, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_push_term_scalar_fmpz(gr_mpoly_t A, gr_srcptr c, const fmpz_struct * exp, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_sort_terms(gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    int gr_mpoly_combine_like_terms(gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    truth_t gr_mpoly_is_canonical(const gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
+    void gr_mpoly_assert_canonical(const gr_mpoly_t A, const mpoly_ctx_t mctx, gr_ctx_t cctx)
diff --git a/src/flint/flintlib/functions/gr_poly.pxd b/src/flint/flintlib/functions/gr_poly.pxd
new file mode 100644
index 00000000..e0c594ad
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_poly.pxd
@@ -0,0 +1,290 @@
+from flint.flintlib.types.flint cimport flint_rand_t, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_poly_t, fmpq_t
+from flint.flintlib.types.fmpz cimport fmpz_poly_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_poly_t, gr_ptr, gr_srcptr, gr_stream_t, gr_vec_t, truth_t
+
+# unknown type gr_method_poly_binary_binary_op
+# unknown type gr_method_poly_binary_op
+# unknown type gr_method_poly_binary_trunc_op
+# unknown type gr_method_poly_gcd_op
+# unknown type gr_method_poly_unary_trunc_op
+# unknown type gr_method_poly_xgcd_op
+
+
+cdef extern from "flint/gr_poly.h":
+    void gr_poly_init(gr_poly_t poly, gr_ctx_t ctx)
+    void gr_poly_init2(gr_poly_t poly, slong len, gr_ctx_t ctx)
+    void gr_poly_clear(gr_poly_t poly, gr_ctx_t ctx)
+    gr_ptr gr_poly_entry_ptr(gr_poly_t poly, slong i, gr_ctx_t ctx)
+    gr_srcptr gr_poly_entry_srcptr(const gr_poly_t poly, slong i, gr_ctx_t ctx)
+    slong gr_poly_length(const gr_poly_t poly, gr_ctx_t ctx)
+    void gr_poly_swap(gr_poly_t poly1, gr_poly_t poly2, gr_ctx_t ctx)
+    void gr_poly_fit_length(gr_poly_t poly, slong len, gr_ctx_t ctx)
+    void _gr_poly_set_length(gr_poly_t poly, slong len, gr_ctx_t ctx)
+    void _gr_poly_normalise(gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_set(gr_poly_t res, const gr_poly_t src, gr_ctx_t ctx)
+    int gr_poly_get_fmpz_poly(gr_poly_t res, const fmpz_poly_t src, gr_ctx_t ctx)
+    int gr_poly_set_fmpq_poly(gr_poly_t res, const fmpq_poly_t src, gr_ctx_t ctx)
+    int gr_poly_set_gr_poly_other(gr_poly_t res, const gr_poly_t x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int _gr_poly_reverse(gr_ptr res, gr_srcptr poly, slong len, slong n, gr_ctx_t ctx)
+    int gr_poly_reverse(gr_poly_t res, const gr_poly_t poly, slong n, gr_ctx_t ctx)
+    int gr_poly_truncate(gr_poly_t res, const gr_poly_t poly, slong newlen, gr_ctx_t ctx)
+    int gr_poly_zero(gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_one(gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_neg_one(gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_gen(gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_write(gr_stream_t out, const gr_poly_t poly, const char * x, gr_ctx_t ctx)
+    int gr_poly_print(const gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_randtest(gr_poly_t poly, flint_rand_t state, slong len, gr_ctx_t ctx)
+    truth_t _gr_poly_equal(gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    truth_t gr_poly_equal(const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    truth_t gr_poly_is_zero(const gr_poly_t poly, gr_ctx_t ctx)
+    truth_t gr_poly_is_one(const gr_poly_t poly, gr_ctx_t ctx)
+    truth_t gr_poly_is_gen(const gr_poly_t poly, gr_ctx_t ctx)
+    truth_t gr_poly_is_scalar(const gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_set_scalar(gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int gr_poly_set_si(gr_poly_t poly, slong c, gr_ctx_t ctx)
+    int gr_poly_set_ui(gr_poly_t poly, ulong c, gr_ctx_t ctx)
+    int gr_poly_set_fmpz(gr_poly_t poly, const fmpz_t c, gr_ctx_t ctx)
+    int gr_poly_set_fmpq(gr_poly_t poly, const fmpq_t c, gr_ctx_t ctx)
+    int gr_poly_set_coeff_scalar(gr_poly_t poly, slong n, gr_srcptr c, gr_ctx_t ctx)
+    int gr_poly_set_coeff_si(gr_poly_t poly, slong n, slong c, gr_ctx_t ctx)
+    int gr_poly_set_coeff_ui(gr_poly_t poly, slong n, ulong c, gr_ctx_t ctx)
+    int gr_poly_set_coeff_fmpz(gr_poly_t poly, slong n, const fmpz_t c, gr_ctx_t ctx)
+    int gr_poly_set_coeff_fmpq(gr_poly_t poly, slong n, const fmpq_t c, gr_ctx_t ctx)
+    int gr_poly_get_coeff_scalar(gr_ptr res, const gr_poly_t poly, slong n, gr_ctx_t ctx)
+    int gr_poly_neg(gr_poly_t res, const gr_poly_t src, gr_ctx_t ctx)
+    int _gr_poly_add(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_add(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_sub(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_sub(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_mul(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_mul(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_mullow_generic(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong len, gr_ctx_t ctx)
+    int _gr_poly_mullow(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong len, gr_ctx_t ctx)
+    int gr_poly_mullow(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong len, gr_ctx_t ctx)
+    int gr_poly_mul_scalar(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_poly_mul_karatsuba(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_mul_karatsuba(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_pow_series_ui_binexp(gr_ptr res, gr_srcptr f, slong flen, ulong exp, slong len, gr_ctx_t ctx)
+    int gr_poly_pow_series_ui_binexp(gr_poly_t res, const gr_poly_t poly, ulong exp, slong len, gr_ctx_t ctx)
+    int _gr_poly_pow_series_ui(gr_ptr res, gr_srcptr f, slong flen, ulong exp, slong len, gr_ctx_t ctx)
+    int gr_poly_pow_series_ui(gr_poly_t res, const gr_poly_t poly, ulong exp, slong len, gr_ctx_t ctx)
+    int _gr_poly_pow_ui_binexp(gr_ptr res, gr_srcptr f, slong flen, ulong exp, gr_ctx_t ctx)
+    int gr_poly_pow_ui_binexp(gr_poly_t res, const gr_poly_t poly, ulong exp, gr_ctx_t ctx)
+    int _gr_poly_pow_ui(gr_ptr res, gr_srcptr f, slong flen, ulong exp, gr_ctx_t ctx)
+    int gr_poly_pow_ui(gr_poly_t res, const gr_poly_t poly, ulong exp, gr_ctx_t ctx)
+    int gr_poly_pow_fmpz(gr_poly_t res, const gr_poly_t poly, const fmpz_t exp, gr_ctx_t ctx)
+    int _gr_poly_pow_series_fmpq_recurrence(gr_ptr h, gr_srcptr f, slong flen, const fmpq_t exp, slong len, int precomp, gr_ctx_t ctx)
+    int gr_poly_pow_series_fmpq_recurrence(gr_poly_t res, const gr_poly_t poly, const fmpq_t exp, slong len, gr_ctx_t ctx)
+    int _gr_poly_shift_left(gr_ptr res, gr_srcptr poly, slong len, slong n, gr_ctx_t ctx)
+    int gr_poly_shift_left(gr_poly_t res, const gr_poly_t poly, slong n, gr_ctx_t ctx)
+    int _gr_poly_shift_right(gr_ptr res, gr_srcptr poly, slong len, slong n, gr_ctx_t ctx)
+    int gr_poly_shift_right(gr_poly_t res, const gr_poly_t poly, slong n, gr_ctx_t ctx)
+    int gr_poly_div_scalar(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_poly_divrem_divconquer_preinv1(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_srcptr invB, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_divrem_divconquer_noinv(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_divrem_divconquer(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_divrem_divconquer(gr_poly_t Q, gr_poly_t R, const gr_poly_t A, const gr_poly_t B, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_divrem_basecase_preinv1(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_srcptr invB, gr_ctx_t ctx)
+    int _gr_poly_divrem_basecase_noinv(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int _gr_poly_divrem_basecase(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_divrem_basecase(gr_poly_t Q, gr_poly_t R, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_divrem_newton(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_divrem_newton(gr_poly_t Q, gr_poly_t R, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_divrem(gr_ptr Q, gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_divrem(gr_poly_t Q, gr_poly_t R, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_div_divconquer_preinv1(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_srcptr invB, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_div_divconquer_noinv(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_div_divconquer(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_div_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_div_basecase_preinv1(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_srcptr invB, gr_ctx_t ctx)
+    int _gr_poly_div_basecase_noinv(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int _gr_poly_div_basecase(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_div_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_div_newton(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_div_newton(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_div(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_div(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_rem(gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_rem(gr_poly_t R, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_inv_series_newton(gr_ptr res, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_inv_series_newton(gr_poly_t res, const gr_poly_t A, slong len, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_inv_series_basecase_preinv1(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr Ainv, slong len, gr_ctx_t ctx)
+    int _gr_poly_inv_series_basecase(gr_ptr res, gr_srcptr A, slong Alen, slong len, gr_ctx_t ctx)
+    int gr_poly_inv_series_basecase(gr_poly_t res, const gr_poly_t A, slong len, gr_ctx_t ctx)
+    int _gr_poly_inv_series(gr_ptr res, gr_srcptr A, slong Alen, slong len, gr_ctx_t ctx)
+    int gr_poly_inv_series(gr_poly_t res, const gr_poly_t A, slong len, gr_ctx_t ctx)
+    int _gr_poly_div_series_newton(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_div_series_newton(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_div_series_divconquer(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_div_series_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_div_series_invmul(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, gr_ctx_t ctx)
+    int gr_poly_div_series_invmul(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
+    int _gr_poly_div_series_basecase_preinv1(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_srcptr Binv, slong len, gr_ctx_t ctx)
+    int _gr_poly_div_series_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+    int _gr_poly_div_series_basecase(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+    int gr_poly_div_series_basecase(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
+    int _gr_poly_div_series(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+    int gr_poly_div_series(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
+    int _gr_poly_divexact_basecase_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx)
+    int gr_poly_divexact_basecase_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_divexact_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx)
+    int gr_poly_divexact_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_divexact_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx)
+    int _gr_poly_divexact_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx)
+    int gr_poly_divexact_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_divexact_series_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+    int _gr_poly_divexact_series_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx)
+    int gr_poly_divexact_series_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx)
+    int _gr_poly_sqrt_series_newton(gr_ptr res, gr_srcptr f, slong flen, slong len, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_sqrt_series_newton(gr_poly_t res, const gr_poly_t f, slong len, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_sqrt_series_basecase(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_sqrt_series_basecase(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_sqrt_series_miller(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_sqrt_series_miller(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_sqrt_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_sqrt_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_rsqrt_series_newton(gr_ptr res, gr_srcptr f, slong flen, slong len, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_rsqrt_series_newton(gr_poly_t res, const gr_poly_t f, slong len, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_rsqrt_series_basecase(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_rsqrt_series_basecase(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_rsqrt_series_miller(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_rsqrt_series_miller(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_rsqrt_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_rsqrt_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_evaluate_rectangular(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_poly_evaluate_rectangular(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_poly_evaluate_modular(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_poly_evaluate_modular(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_poly_evaluate_horner(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_poly_evaluate_horner(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_poly_evaluate(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
+    int gr_poly_evaluate(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_poly_evaluate_other_horner(gr_ptr res, gr_srcptr f, slong len, const gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int gr_poly_evaluate_other_horner(gr_ptr res, const gr_poly_t f, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int _gr_poly_evaluate_other_rectangular(gr_ptr res, gr_srcptr f, slong len, const gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int gr_poly_evaluate_other_rectangular(gr_ptr res, const gr_poly_t f, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int _gr_poly_evaluate_other(gr_ptr res, gr_srcptr f, slong len, const gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    int gr_poly_evaluate_other(gr_ptr res, const gr_poly_t f, gr_srcptr x, gr_ctx_t x_ctx, gr_ctx_t ctx)
+    gr_ptr * _gr_poly_tree_alloc(slong len, gr_ctx_t ctx)
+    void _gr_poly_tree_free(gr_ptr * tree, slong len, gr_ctx_t ctx)
+    int _gr_poly_tree_build(gr_ptr * tree, gr_srcptr roots, slong len, gr_ctx_t ctx)
+    int _gr_poly_evaluate_vec_fast_precomp(gr_ptr vs, gr_srcptr poly, slong plen, gr_ptr * tree, slong len, gr_ctx_t ctx)
+    int _gr_poly_evaluate_vec_fast(gr_ptr ys, gr_srcptr poly, slong plen, gr_srcptr xs, slong n, gr_ctx_t ctx)
+    int gr_poly_evaluate_vec_fast(gr_vec_t ys, const gr_poly_t poly, const gr_vec_t xs, gr_ctx_t ctx)
+    int _gr_poly_evaluate_vec_iter(gr_ptr ys, gr_srcptr poly, slong plen, gr_srcptr xs, slong n, gr_ctx_t ctx)
+    int gr_poly_evaluate_vec_iter(gr_vec_t ys, const gr_poly_t poly, const gr_vec_t xs, gr_ctx_t ctx)
+    int _gr_poly_taylor_shift_horner(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_poly_taylor_shift_horner(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_poly_taylor_shift_divconquer(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_poly_taylor_shift_divconquer(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_poly_taylor_shift_convolution(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_poly_taylor_shift_convolution(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_poly_taylor_shift(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int gr_poly_taylor_shift(gr_poly_t res, const gr_poly_t poly, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_poly_compose_horner(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_compose_horner(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_compose_divconquer(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_compose_divconquer(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_compose(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_compose(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, gr_ctx_t ctx)
+    int _gr_poly_compose_series_horner(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
+    int gr_poly_compose_series_horner(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx)
+    int _gr_poly_compose_series_brent_kung(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
+    int gr_poly_compose_series_brent_kung(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx)
+    int _gr_poly_compose_series_divconquer(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
+    int gr_poly_compose_series_divconquer(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx)
+    int _gr_poly_compose_series(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, slong n, gr_ctx_t ctx)
+    int gr_poly_compose_series(gr_poly_t res, const gr_poly_t poly1, const gr_poly_t poly2, slong n, gr_ctx_t ctx)
+    int _gr_poly_revert_series_lagrange(gr_ptr res, gr_srcptr f, slong flen, slong n, gr_ctx_t ctx)
+    int gr_poly_revert_series_lagrange(gr_poly_t res, const gr_poly_t f, slong n, gr_ctx_t ctx)
+    int _gr_poly_revert_series_lagrange_fast(gr_ptr res, gr_srcptr f, slong flen, slong n, gr_ctx_t ctx)
+    int gr_poly_revert_series_lagrange_fast(gr_poly_t res, const gr_poly_t f, slong n, gr_ctx_t ctx)
+    int _gr_poly_revert_series_newton(gr_ptr res, gr_srcptr f, slong flen, slong n, gr_ctx_t ctx)
+    int gr_poly_revert_series_newton(gr_poly_t res, const gr_poly_t f, slong n, gr_ctx_t ctx)
+    int _gr_poly_revert_series(gr_ptr res, gr_srcptr f, slong flen, slong n, gr_ctx_t ctx)
+    int gr_poly_revert_series(gr_poly_t res, const gr_poly_t f, slong n, gr_ctx_t ctx)
+    int _gr_poly_derivative(gr_ptr res, gr_srcptr poly, slong len, gr_ctx_t ctx)
+    int gr_poly_derivative(gr_poly_t res, const gr_poly_t poly, gr_ctx_t ctx)
+    int _gr_poly_nth_derivative(gr_ptr res, gr_srcptr poly, ulong n, slong len, gr_ctx_t ctx)
+    int gr_poly_nth_derivative(gr_poly_t res, const gr_poly_t poly, ulong n, gr_ctx_t ctx)
+    int _gr_poly_integral(gr_ptr res, gr_srcptr poly, slong len, gr_ctx_t ctx)
+    int gr_poly_integral(gr_poly_t res, const gr_poly_t poly, gr_ctx_t ctx)
+    int _gr_poly_make_monic(gr_ptr res, gr_srcptr poly, slong len, gr_ctx_t ctx)
+    int gr_poly_make_monic(gr_poly_t res, const gr_poly_t src, gr_ctx_t ctx)
+    truth_t _gr_poly_is_monic(gr_srcptr poly, slong len, gr_ctx_t ctx)
+    truth_t gr_poly_is_monic(const gr_poly_t res, gr_ctx_t ctx)
+    int _gr_poly_hgcd(gr_ptr r, slong * sgn, gr_ptr * M, slong * lenM, gr_ptr A, slong * lenA, gr_ptr B, slong * lenB, gr_srcptr a, slong lena, gr_srcptr b, slong lenb, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_gcd_hgcd(gr_ptr G, slong * _lenG, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong inner_cutoff, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_gcd_hgcd(gr_poly_t G, const gr_poly_t A, const gr_poly_t B, slong inner_cutoff, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_gcd_euclidean(gr_ptr G, slong * lenG, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_gcd_euclidean(gr_poly_t G, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_gcd_generic(gr_ptr G, slong * lenG, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int _gr_poly_gcd(gr_ptr G, slong * lenG, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_gcd(gr_poly_t G, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_xgcd_euclidean(slong * lenG, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_xgcd_euclidean(gr_poly_t G, gr_poly_t S, gr_poly_t T, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_xgcd_hgcd(slong * Glen, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong hgcd_cutoff, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_xgcd_hgcd(gr_poly_t G, gr_poly_t S, gr_poly_t T, const gr_poly_t A, const gr_poly_t B, slong hgcd_cutoff, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_xgcd_generic(slong * lenG, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int _gr_poly_xgcd(slong * lenG, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx)
+    int gr_poly_xgcd(gr_poly_t G, gr_poly_t S, gr_poly_t T, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx)
+    int _gr_poly_resultant_euclidean(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_resultant_euclidean(gr_ptr res, const gr_poly_t f, const gr_poly_t g, gr_ctx_t ctx)
+    int _gr_poly_resultant_hgcd(gr_ptr res, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong inner_cutoff, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_resultant_hgcd(gr_ptr res, const gr_poly_t f, const gr_poly_t g, slong inner_cutoff, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_resultant_sylvester(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_resultant_sylvester(gr_ptr res, const gr_poly_t f, const gr_poly_t g, gr_ctx_t ctx)
+    int _gr_poly_resultant_small(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_resultant_small(gr_ptr res, const gr_poly_t f, const gr_poly_t g, gr_ctx_t ctx)
+    int _gr_poly_resultant(gr_ptr res, gr_srcptr poly1, slong len1, gr_srcptr poly2, slong len2, gr_ctx_t ctx)
+    int gr_poly_resultant(gr_ptr res, const gr_poly_t f, const gr_poly_t g, gr_ctx_t ctx)
+    int gr_poly_factor_squarefree(gr_ptr c, gr_vec_t fac, gr_vec_t exp, const gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_squarefree_part(gr_poly_t res, const gr_poly_t poly, gr_ctx_t ctx)
+    int gr_poly_roots(gr_vec_t roots, gr_vec_t mult, const gr_poly_t poly, int flags, gr_ctx_t ctx)
+    int gr_poly_roots_other(gr_vec_t roots, gr_vec_t mult, const gr_poly_t poly, gr_ctx_t poly_ctx, int flags, gr_ctx_t ctx)
+    int _gr_poly_refine_roots_aberth(gr_ptr w, gr_srcptr f, gr_srcptr f_prime, slong deg, gr_srcptr z, int progressive, gr_ctx_t ctx)
+    int _gr_poly_refine_roots_wdk(gr_ptr w, gr_srcptr f, slong deg, gr_srcptr z, int progressive, gr_ctx_t ctx)
+    int _gr_poly_asin_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_asin_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_asinh_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_asinh_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_acos_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_acos_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_acosh_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_acosh_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_atan_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_atan_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_atanh_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_atanh_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_log_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_log_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_log1p_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_log1p_series(gr_poly_t res, const gr_poly_t f, slong len, gr_ctx_t ctx)
+    int _gr_poly_exp_series_basecase(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+    int gr_poly_exp_series_basecase(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+    int _gr_poly_exp_series_basecase_mul(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+    int gr_poly_exp_series_basecase_mul(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+    int _gr_poly_exp_series_newton(gr_ptr f, gr_ptr g, gr_srcptr h, slong hlen, slong n, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_exp_series_newton(gr_poly_t f, const gr_poly_t h, slong n, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_exp_series_generic(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+    int _gr_poly_exp_series(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx)
+    int gr_poly_exp_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+    int _gr_poly_sin_cos_series_basecase(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, int times_pi, gr_ctx_t ctx)
+    int gr_poly_sin_cos_series_basecase(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, int times_pi, gr_ctx_t ctx)
+    int _gr_poly_sin_cos_series_tangent(gr_ptr s, gr_ptr c, gr_srcptr h, slong hlen, slong n, int times_pi, gr_ctx_t ctx)
+    int gr_poly_sin_cos_series_tangent(gr_poly_t s, gr_poly_t c, const gr_poly_t h, slong n, int times_pi, gr_ctx_t ctx)
+    int _gr_poly_tan_series_basecase(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+    int gr_poly_tan_series_basecase(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+    int _gr_poly_tan_series_newton(gr_ptr f, gr_srcptr h, slong hlen, slong n, slong cutoff, gr_ctx_t ctx)
+    int gr_poly_tan_series_newton(gr_poly_t f, const gr_poly_t h, slong n, slong cutoff, gr_ctx_t ctx)
+    int _gr_poly_tan_series(gr_ptr f, gr_srcptr h, slong hlen, slong n, gr_ctx_t ctx)
+    int gr_poly_tan_series(gr_poly_t f, const gr_poly_t h, slong n, gr_ctx_t ctx)
+    # void _gr_poly_test_mullow(gr_method_poly_binary_trunc_op mullow_impl, gr_method_poly_binary_trunc_op mullow_ref, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void _gr_poly_test_divrem(gr_method_poly_binary_binary_op divrem_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void _gr_poly_test_div(gr_method_poly_binary_op div_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void _gr_poly_test_inv_series(gr_method_poly_unary_trunc_op inv_series_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void _gr_poly_test_div_series(gr_method_poly_binary_trunc_op div_series_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void _gr_poly_test_gcd(gr_method_poly_gcd_op gcd_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+    # void _gr_poly_test_xgcd(gr_method_poly_xgcd_op xgcd_impl, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
diff --git a/src/flint/flintlib/functions/gr_special.pxd b/src/flint/flintlib/functions/gr_special.pxd
new file mode 100644
index 00000000..caa12236
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_special.pxd
@@ -0,0 +1,228 @@
+from flint.flintlib.types.dirichlet cimport dirichlet_char_t, dirichlet_group_t
+from flint.flintlib.types.flint cimport fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_ptr, gr_srcptr, gr_vec_t
+
+
+
+cdef extern from "flint/gr_special.h":
+    int gr_pi(gr_ptr res, gr_ctx_t ctx)
+    int gr_euler(gr_ptr res, gr_ctx_t ctx)
+    int gr_catalan(gr_ptr res, gr_ctx_t ctx)
+    int gr_khinchin(gr_ptr res, gr_ctx_t ctx)
+    int gr_glaisher(gr_ptr res, gr_ctx_t ctx)
+    int gr_exp(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_expm1(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_exp2(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_exp10(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_exp_pi_i(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log1p(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log2(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log10(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log_pi_i(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sin(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cos(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sin_cos(gr_ptr res1, gr_ptr res2, gr_srcptr x, gr_ctx_t ctx)
+    int gr_tan(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cot(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sec(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_csc(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sin_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cos_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sin_cos_pi(gr_ptr res1, gr_ptr res2, gr_srcptr x, gr_ctx_t ctx)
+    int gr_tan_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cot_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sec_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_csc_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sinc(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sinc_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sinh(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cosh(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sinh_cosh(gr_ptr res1, gr_ptr res2, gr_srcptr x, gr_ctx_t ctx)
+    int gr_tanh(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_coth(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sech(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_csch(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_asin(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acos(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_atan(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_atan2(gr_ptr res, gr_srcptr y, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acot(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_asec(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acsc(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_asin_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acos_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_atan_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acot_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_asec_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acsc_pi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_asinh(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acosh(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_atanh(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acoth(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_asech(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_acsch(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_lambertw(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_lambertw_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t k, gr_ctx_t ctx)
+    int gr_fac(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fac_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_fac_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_fac_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_rfac(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_rfac_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_rfac_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_rfac_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_bin(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bin_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_bin_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_bin_vec(gr_ptr res, gr_srcptr x, slong len, gr_ctx_t ctx)
+    int gr_bin_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    int gr_rising(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_rising_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_falling(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_falling_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
+    int gr_gamma(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_gamma_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_gamma_fmpq(gr_ptr res, const fmpq_t x, gr_ctx_t ctx)
+    int gr_rgamma(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_lgamma(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_digamma(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_barnes_g(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log_barnes_g(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_beta(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_doublefac(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_doublefac_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_harmonic(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_harmonic_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_bernoulli_ui(gr_ptr res, ulong n, gr_ctx_t ctx)
+    int gr_bernoulli_fmpz(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_bernoulli_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_eulernum_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_eulernum_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_eulernum_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_fib_ui(gr_ptr res, ulong n, gr_ctx_t ctx)
+    int gr_fib_fmpz(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_fib_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_stirling_s1u_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_stirling_s1_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_stirling_s2_uiui(gr_ptr res, ulong x, ulong y, gr_ctx_t ctx)
+    int gr_stirling_s1u_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    int gr_stirling_s1_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    int gr_stirling_s2_ui_vec(gr_ptr res, ulong x, slong len, gr_ctx_t ctx)
+    int gr_bellnum_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_bellnum_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_bellnum_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_partitions_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_partitions_fmpz(gr_ptr res, const fmpz_t x, gr_ctx_t ctx)
+    int gr_partitions_vec(gr_ptr res, slong len, gr_ctx_t ctx)
+    int gr_erf(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_erfc(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_erfcx(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_erfi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_erfinv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_erfcinv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_fresnel_s(gr_ptr res, gr_srcptr x, int normalized, gr_ctx_t ctx)
+    int gr_fresnel_c(gr_ptr res, gr_srcptr x, int normalized, gr_ctx_t ctx)
+    int gr_fresnel(gr_ptr res1, gr_ptr res2, gr_srcptr x, int normalized, gr_ctx_t ctx)
+    int gr_gamma_upper(gr_ptr res, gr_srcptr x, gr_srcptr y, int regularized, gr_ctx_t ctx)
+    int gr_gamma_lower(gr_ptr res, gr_srcptr x, gr_srcptr y, int regularized, gr_ctx_t ctx)
+    int gr_beta_lower(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, int regularized, gr_ctx_t ctx)
+    int gr_exp_integral(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_exp_integral_ei(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sin_integral(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cos_integral(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_sinh_integral(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_cosh_integral(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_log_integral(gr_ptr res, gr_srcptr x, int offset, gr_ctx_t ctx)
+    int gr_dilog(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_chebyshev_t_fmpz(gr_ptr res, const fmpz_t n, gr_srcptr x, gr_ctx_t ctx)
+    int gr_chebyshev_t(gr_ptr res, gr_srcptr n, gr_srcptr x, gr_ctx_t ctx)
+    int gr_chebyshev_u_fmpz(gr_ptr res, const fmpz_t n, gr_srcptr x, gr_ctx_t ctx)
+    int gr_chebyshev_u(gr_ptr res, gr_srcptr n, gr_srcptr x, gr_ctx_t ctx)
+    int gr_jacobi_p(gr_ptr res, gr_srcptr n, gr_srcptr a, gr_srcptr b, gr_srcptr z, gr_ctx_t ctx)
+    int gr_gegenbauer_c(gr_ptr res, gr_srcptr n, gr_srcptr m, gr_srcptr z, gr_ctx_t ctx)
+    int gr_laguerre_l(gr_ptr res, gr_srcptr n, gr_srcptr m, gr_srcptr z, gr_ctx_t ctx)
+    int gr_hermite_h(gr_ptr res, gr_srcptr n, gr_srcptr z, gr_ctx_t ctx)
+    int gr_legendre_p(gr_ptr res, gr_srcptr n, gr_srcptr m, gr_srcptr z, int type, gr_ctx_t ctx)
+    int gr_legendre_q(gr_ptr res, gr_srcptr n, gr_srcptr m, gr_srcptr z, int type, gr_ctx_t ctx)
+    int gr_spherical_y_si(gr_ptr res, slong n, slong m, gr_srcptr theta, gr_srcptr phi, gr_ctx_t ctx)
+    int gr_legendre_p_root_ui(gr_ptr root, gr_ptr weight, ulong n, ulong k, gr_ctx_t ctx)
+    int gr_bessel_j(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bessel_y(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bessel_i(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bessel_k(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bessel_j_y(gr_ptr res1, gr_ptr res2, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bessel_i_scaled(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_bessel_k_scaled(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_airy(gr_ptr res1, gr_ptr res2, gr_ptr res3, gr_ptr res4, gr_srcptr x, gr_ctx_t ctx)
+    int gr_airy_ai(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_airy_bi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_airy_ai_prime(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_airy_bi_prime(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_airy_ai_zero(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_airy_bi_zero(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_airy_ai_prime_zero(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_airy_bi_prime_zero(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_coulomb(gr_ptr res1, gr_ptr res2, gr_ptr res3, gr_ptr res4, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_ctx_t ctx)
+    int gr_coulomb_f(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_ctx_t ctx)
+    int gr_coulomb_g(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_ctx_t ctx)
+    int gr_coulomb_hpos(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_ctx_t ctx)
+    int gr_coulomb_hneg(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_ctx_t ctx)
+    int gr_hypgeom_0f1(gr_ptr res, gr_srcptr a, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_hypgeom_1f1(gr_ptr res, gr_srcptr a, gr_srcptr b, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_hypgeom_u(gr_ptr res, gr_srcptr a, gr_srcptr b, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_hypgeom_2f1(gr_ptr res, gr_srcptr a, gr_srcptr b, gr_srcptr c, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_hypgeom_pfq(gr_ptr res, const gr_vec_t a, const gr_vec_t b, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_zeta(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_zeta_ui(gr_ptr res, ulong x, gr_ctx_t ctx)
+    int gr_hurwitz_zeta(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_polygamma(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_polylog(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_lerch_phi(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_ctx_t ctx)
+    int gr_stieltjes(gr_ptr res, const fmpz_t x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_dirichlet_eta(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_riemann_xi(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_zeta_zero(gr_ptr res, const fmpz_t n, gr_ctx_t ctx)
+    int gr_zeta_zero_vec(gr_ptr res, const fmpz_t n, slong len, gr_ctx_t ctx)
+    int gr_zeta_nzeros(gr_ptr res, gr_srcptr t, gr_ctx_t ctx)
+    int gr_dirichlet_chi_fmpz(gr_ptr res, const dirichlet_group_t G, const dirichlet_char_t chi, const fmpz_t n, gr_ctx_t ctx)
+    int gr_dirichlet_chi_vec(gr_ptr res, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, gr_ctx_t ctx)
+    int gr_dirichlet_l(gr_ptr res, const dirichlet_group_t G, const dirichlet_char_t chi, gr_srcptr s, gr_ctx_t ctx)
+    int gr_dirichlet_l_all(gr_vec_t res, const dirichlet_group_t G, gr_srcptr s, gr_ctx_t ctx)
+    int gr_dirichlet_hardy_theta(gr_ptr res, const dirichlet_group_t G, const dirichlet_char_t chi, gr_srcptr t, gr_ctx_t ctx)
+    int gr_dirichlet_hardy_z(gr_ptr res, const dirichlet_group_t G, const dirichlet_char_t chi, gr_srcptr t, gr_ctx_t ctx)
+    int gr_agm1(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
+    int gr_agm(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+    int gr_elliptic_k(gr_ptr res, gr_srcptr m, gr_ctx_t ctx)
+    int gr_elliptic_e(gr_ptr res, gr_srcptr m, gr_ctx_t ctx)
+    int gr_elliptic_pi(gr_ptr res, gr_srcptr n, gr_srcptr m, gr_ctx_t ctx)
+    int gr_elliptic_f(gr_ptr res, gr_srcptr phi, gr_srcptr m, int pi, gr_ctx_t ctx)
+    int gr_elliptic_e_inc(gr_ptr res, gr_srcptr phi, gr_srcptr m, int pi, gr_ctx_t ctx)
+    int gr_elliptic_pi_inc(gr_ptr res, gr_srcptr n, gr_srcptr phi, gr_srcptr m, int pi, gr_ctx_t ctx)
+    int gr_carlson_rc(gr_ptr res, gr_srcptr x, gr_srcptr y, int flags, gr_ctx_t ctx)
+    int gr_carlson_rf(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_carlson_rd(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_carlson_rg(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, int flags, gr_ctx_t ctx)
+    int gr_carlson_rj(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_srcptr z, gr_srcptr w, int flags, gr_ctx_t ctx)
+    int gr_jacobi_theta(gr_ptr res1, gr_ptr res2, gr_ptr res3, gr_ptr res4, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_jacobi_theta_1(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_jacobi_theta_2(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_jacobi_theta_3(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_jacobi_theta_4(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_dedekind_eta(gr_ptr res, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_dedekind_eta_q(gr_ptr res, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_modular_j(gr_ptr res, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_modular_lambda(gr_ptr res, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_modular_delta(gr_ptr res, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_hilbert_class_poly(gr_ptr res, slong D, gr_srcptr x, gr_ctx_t ctx)
+    int gr_eisenstein_e(gr_ptr res, ulong n, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_eisenstein_g(gr_ptr res, ulong n, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_eisenstein_g_vec(gr_ptr res, gr_srcptr tau, slong len, gr_ctx_t ctx)
+    int gr_elliptic_invariants(gr_ptr res1, gr_ptr res2, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_elliptic_roots(gr_ptr res1, gr_ptr res2, gr_ptr res3, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_weierstrass_p(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_weierstrass_p_prime(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_weierstrass_p_inv(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_weierstrass_zeta(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
+    int gr_weierstrass_sigma(gr_ptr res, gr_srcptr z, gr_srcptr tau, gr_ctx_t ctx)
diff --git a/src/flint/flintlib/functions/gr_vec.pxd b/src/flint/flintlib/functions/gr_vec.pxd
new file mode 100644
index 00000000..3f876dfa
--- /dev/null
+++ b/src/flint/flintlib/functions/gr_vec.pxd
@@ -0,0 +1,112 @@
+from flint.flintlib.types.flint cimport flint_rand_t, fmpz_struct, fmpz_t, slong, ulong
+from flint.flintlib.types.fmpq cimport fmpq_t
+from flint.flintlib.types.gr cimport gr_ctx_t, gr_ptr, gr_srcptr, gr_stream_t, gr_vec_t, truth_t
+
+
+# .. macro:: GR_VEC_ENTRY(vec, i, sz)
+# .. macro:: GR_ENTRY(vec, i, size)
+
+cdef extern from "flint/gr_vec.h":
+    void gr_vec_init(gr_vec_t vec, slong len, gr_ctx_t ctx)
+    void gr_vec_clear(gr_vec_t vec, gr_ctx_t ctx)
+    gr_ptr gr_vec_entry_ptr(gr_vec_t vec, slong i, gr_ctx_t ctx)
+    slong gr_vec_length(const gr_vec_t vec, gr_ctx_t ctx)
+    void gr_vec_fit_length(gr_vec_t vec, slong len, gr_ctx_t ctx)
+    void gr_vec_set_length(gr_vec_t vec, slong len, gr_ctx_t ctx)
+    int gr_vec_set(gr_vec_t res, const gr_vec_t src, gr_ctx_t ctx)
+    int gr_vec_append(gr_vec_t vec, gr_srcptr x, gr_ctx_t ctx)
+    int _gr_vec_write(gr_stream_t out, gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int gr_vec_write(gr_stream_t out, const gr_vec_t vec, gr_ctx_t ctx)
+    int gr_vec_print(const gr_vec_t vec, gr_ctx_t ctx)
+    void _gr_vec_init(gr_ptr vec, slong len, gr_ctx_t ctx)
+    void _gr_vec_clear(gr_ptr vec, slong len, gr_ctx_t ctx)
+    void _gr_vec_swap(gr_ptr vec1, gr_ptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_randtest(gr_ptr res, flint_rand_t state, slong len, gr_ctx_t ctx)
+    int _gr_vec_set(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
+    truth_t _gr_vec_equal(gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_zero(gr_ptr vec, slong len, gr_ctx_t ctx)
+    truth_t _gr_vec_is_zero(gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int _gr_vec_normalise(slong * res, gr_srcptr vec, slong len, gr_ctx_t ctx)
+    slong _gr_vec_normalise_weak(gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int _gr_vec_neg(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
+    int _gr_vec_add(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int _gr_vec_sub(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int _gr_vec_mul(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int _gr_vec_div(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int _gr_vec_divexact(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int _gr_vec_pow(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
+    int _gr_vec_add_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_sub_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_div_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_divexact_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_pow_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_scalar_add_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_sub_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_mul_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_div_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_divexact_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_pow_vec(gr_ptr vec1, gr_srcptr c, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_add_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int _gr_vec_sub_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int _gr_vec_mul_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int _gr_vec_div_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int _gr_vec_divexact_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int _gr_vec_pow_other(gr_ptr vec1, gr_srcptr vec2, gr_srcptr vec3, gr_ctx_t ctx3, slong len, gr_ctx_t ctx)
+    int _gr_other_add_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int _gr_other_sub_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int _gr_other_mul_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int _gr_other_div_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int _gr_other_divexact_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int _gr_other_pow_vec(gr_ptr vec1, gr_srcptr vec2, gr_ctx_t ctx2, gr_srcptr vec3, slong len, gr_ctx_t ctx)
+    int _gr_vec_add_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int _gr_vec_sub_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int _gr_vec_div_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int _gr_vec_divexact_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int _gr_vec_pow_scalar_other(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t cctx, gr_ctx_t ctx)
+    int _gr_scalar_other_add_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_other_sub_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_other_mul_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_other_div_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_other_divexact_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_scalar_other_pow_vec(gr_ptr vec1, gr_srcptr c, gr_ctx_t cctx, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_add_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_sub_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_div_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_divexact_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_pow_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_add_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int _gr_vec_sub_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int _gr_vec_div_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int _gr_vec_divexact_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int _gr_vec_pow_scalar_ui(gr_ptr vec1, gr_srcptr vec2, slong len, ulong c, gr_ctx_t ctx)
+    int _gr_vec_add_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int _gr_vec_sub_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int _gr_vec_div_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int _gr_vec_divexact_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int _gr_vec_pow_scalar_fmpz(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpz_t c, gr_ctx_t ctx)
+    int _gr_vec_add_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int _gr_vec_sub_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int _gr_vec_div_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int _gr_vec_divexact_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int _gr_vec_pow_scalar_fmpq(gr_ptr vec1, gr_srcptr vec2, slong len, const fmpq_t c, gr_ctx_t ctx)
+    int _gr_vec_addmul_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_submul_scalar(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
+    int _gr_vec_addmul_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_submul_scalar_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_mul_scalar_2exp_si(gr_ptr res, gr_srcptr vec, slong len, slong c, gr_ctx_t ctx)
+    int _gr_vec_sum(gr_ptr res, gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int _gr_vec_product(gr_ptr res, gr_srcptr vec, slong len, gr_ctx_t ctx)
+    int _gr_vec_dot(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_dot_si(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, const slong * vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_dot_ui(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, const ulong * vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_dot_fmpz(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, const fmpz_struct * vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_dot_rev(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
+    int _gr_vec_step(gr_ptr vec, gr_srcptr start, gr_srcptr step, slong len, gr_ctx_t ctx)
+    int _gr_vec_reciprocals(gr_ptr res, slong len, gr_ctx_t ctx)
+    int _gr_vec_set_powers(gr_ptr res, gr_srcptr x, slong len, gr_ctx_t ctx)
diff --git a/src/flint/flintlib/types/flint.pxd b/src/flint/flintlib/types/flint.pxd
index 8db56127..77488e28 100644
--- a/src/flint/flintlib/types/flint.pxd
+++ b/src/flint/flintlib/types/flint.pxd
@@ -61,14 +61,23 @@ cdef extern from "flint/flint.h":
     #define SIZEOF_ULONG sizeof(ulong)
     """
     int SIZEOF_ULONG
+
+    ctypedef struct __FLINT_FILE:
+        pass
+    ctypedef __FLINT_FILE FLINT_FILE
+
     const char * FLINT_VERSION
     const int __FLINT_RELEASE
+
     const int FLINT_BITS
+
     ctypedef void * flint_rand_t
     void flint_rand_init(flint_rand_t state)
     void flint_rand_clear(flint_rand_t state)
+
     void flint_set_num_threads(long)
     long flint_get_num_threads()
+
     void flint_cleanup()
 
     ctypedef struct nmod_t:
diff --git a/src/flint/flintlib/types/gr.pxd b/src/flint/flintlib/types/gr.pxd
new file mode 100644
index 00000000..b62cc01e
--- /dev/null
+++ b/src/flint/flintlib/types/gr.pxd
@@ -0,0 +1,76 @@
+from flint.flintlib.types.flint cimport (
+    FLINT_FILE,
+    slong,
+    ulong,
+    flint_bitcnt_t,
+)
+
+
+cdef extern from "flint/gr.h":
+
+    ctypedef int truth_t
+    cdef int T_TRUE
+    cdef int T_FALSE
+    cdef int T_UNKNOWN
+
+    # Macros
+    cdef int GR_SUCCESS
+    cdef int GR_DOMAIN
+    cdef int GR_UNABLE
+
+    ctypedef struct gr_stream_struct:
+        FLINT_FILE * fp
+        char * s
+        slong len
+        slong alloc
+
+    ctypedef gr_stream_struct gr_stream_t[1]
+
+    ctypedef int (*gr_funcptr)()
+
+    cdef const int GR_CTX_STRUCT_DATA_BYTES
+
+    cdef struct gr_ctx_struct:
+        char data[GR_CTX_STRUCT_DATA_BYTES]
+        ulong which_ring
+        slong sizeof_elem
+        gr_funcptr * methods
+        ulong size_limit
+
+    ctypedef gr_ctx_struct gr_ctx_t[1]
+
+    ctypedef void * gr_ptr
+    ctypedef const void * gr_srcptr
+    ctypedef void * gr_ctx_ptr
+
+    ctypedef struct gr_vec_struct:
+        gr_ptr entries
+        slong alloc
+        slong length
+
+    ctypedef gr_vec_struct gr_vec_t[1]
+
+    ctypedef struct gr_mat_struct:
+        gr_ptr entries
+        slong r
+        slong c
+        gr_ptr * rows
+
+    ctypedef gr_mat_struct gr_mat_t[1]
+
+    ctypedef struct gr_poly_struct:
+        gr_ptr coeffs
+        slong alloc
+        slong length
+
+    ctypedef gr_poly_struct gr_poly_t[1]
+
+    ctypedef struct gr_mpoly_struct:
+        gr_ptr coeffs
+        ulong * exps
+        slong length
+        flint_bitcnt_t bits  # number of bits per exponent
+        slong coeffs_alloc   # abs size in ulong units
+        slong exps_alloc     # abs size in ulong units
+
+    ctypedef gr_mpoly_struct gr_mpoly_t[1]
diff --git a/src/flint/types/_gr.pxd b/src/flint/types/_gr.pxd
new file mode 100644
index 00000000..f615e991
--- /dev/null
+++ b/src/flint/types/_gr.pxd
@@ -0,0 +1,1219 @@
+cimport cython
+cimport libc.stdlib
+
+from flint.flintlib.types.flint cimport (
+    slong,
+    ulong,
+    fmpz_t,
+)
+from flint.flintlib.types.fmpz cimport (
+    fmpz_poly_t,
+)
+from flint.flintlib.types.fmpq cimport (
+    fmpq_poly_t,
+)
+from flint.flintlib.types.mpoly cimport (
+    ordering_t,
+)
+from flint.flintlib.functions.ulong_extras cimport (
+    n_is_prime,
+)
+from flint.flintlib.functions.fmpz cimport (
+    fmpz_init_set,
+)
+from flint.flintlib.functions.fmpz_poly cimport (
+    fmpz_poly_init,
+    fmpz_poly_set,
+)
+from flint.flintlib.functions.fmpq_poly cimport (
+    fmpq_poly_init,
+    fmpq_poly_set,
+)
+
+
+from flint.flintlib.types.gr cimport (
+    truth_t,
+    T_TRUE,
+    T_FALSE,
+    T_UNKNOWN,
+
+    GR_SUCCESS,
+    GR_DOMAIN,
+    GR_UNABLE,
+
+    gr_ctx_t,
+    gr_ptr,
+    gr_vec_t,
+)
+from flint.flintlib.functions.gr_domains cimport (
+    gr_ctx_init_fmpz,
+    gr_ctx_init_fmpq,
+    gr_ctx_init_fmpzi,
+    gr_ctx_init_fexpr,
+
+    gr_ctx_init_nmod,
+    gr_ctx_init_fmpz_mod,
+    # gr_ctx_set_is_field,
+
+    gr_ctx_init_fq,
+    gr_ctx_init_fq_nmod,
+    gr_ctx_init_fq_zech,
+
+    gr_ctx_init_nf,
+    gr_ctx_init_nf_fmpz_poly,
+
+    gr_ctx_init_real_qqbar,
+    gr_ctx_init_complex_qqbar,
+    _gr_ctx_qqbar_set_limits,
+
+    gr_ctx_init_real_ca,
+    gr_ctx_init_complex_ca,
+    gr_ctx_init_real_algebraic_ca,
+    gr_ctx_init_complex_algebraic_ca,
+    gr_ctx_init_complex_extended_ca,
+
+    gr_ctx_init_real_float_arf,
+    gr_ctx_init_complex_float_acf,
+
+    gr_ctx_init_real_arb,
+    gr_ctx_init_complex_acb,
+
+    gr_ctx_init_gr_poly,
+    gr_ctx_init_gr_mpoly,
+    # gr_ctx_init_fmpz_mpoly_q,
+    # gr_ctx_init_series_mod_gr_poly,
+    gr_ctx_init_gr_series,
+)
+from flint.flintlib.functions.gr cimport (
+    gr_heap_init,
+    gr_set_str,
+    gr_get_str,
+    gr_set,
+    gr_get_si,
+
+    gr_zero,
+    gr_one,
+    gr_gen,
+    gr_gens,
+    gr_gens_recursive,
+    gr_ctx_set_gen_names,
+
+    gr_i,
+    gr_pos_inf,
+    gr_neg_inf,
+    gr_uinf,
+    gr_undefined,
+    gr_unknown,
+
+    gr_is_zero,
+    gr_is_one,
+    gr_is_neg_one,
+    # gr_is_integer,
+    # gr_is_rational,
+
+    gr_equal,
+
+    gr_neg,
+    gr_add,
+    gr_add_si,
+    gr_sub,
+    gr_sub_si,
+    gr_mul,
+    gr_mul_si,
+    gr_inv,
+    gr_div,
+    gr_div_si,
+    gr_pow_si,
+
+    gr_is_square,
+    gr_sqrt,
+    gr_rsqrt,
+
+    gr_gcd,
+    gr_lcm,
+    gr_factor,
+
+    gr_floor,
+    gr_ceil,
+    gr_trunc,
+    gr_nint,
+
+    gr_abs,
+    gr_conj,
+    gr_re,
+    gr_im,
+    gr_sgn,
+    gr_csgn,
+    gr_arg,
+
+    # gr_le,
+    # gr_lt,
+    # gr_ge,
+    # gr_gt,
+
+    gr_numerator,
+    gr_denominator,
+)
+from flint.flintlib.functions.gr_vec cimport (
+    gr_vec_init,
+    gr_vec_clear,
+    gr_vec_length,
+    gr_vec_entry_ptr,
+)
+
+from flint.flint_base.flint_base cimport (
+    flint_ctx,
+    flint_scalar,
+    ordering_py_to_c,
+)
+
+from flint.types.fmpz cimport fmpz
+from flint.types.fmpz_poly cimport fmpz_poly
+from flint.types.fmpq_poly cimport fmpq_poly
+
+
+# XXX: Can't import from Python modules in a .pxd file
+# from flint.utils.flint_exceptions import DomainError, UnableError, UknownError
+
+
+cdef inline truth_to_py(truth_t t):
+    if t == T_TRUE:
+        return True
+    elif t == T_FALSE:
+        return False
+    else:
+        return None
+
+
+cdef inline truth_to_bool(truth_t t):
+    if t == T_TRUE:
+        return True
+    elif t == T_FALSE:
+        return False
+    else:
+        # raise UnknownError("Unknown truth value")
+        raise AssertionError("Unknown truth value")
+
+
+@cython.no_gc
+cdef class gr_ctx(flint_ctx):
+    cdef gr_ctx_t ctx_t
+    cdef bint _init
+
+    @cython.final
+    cdef inline gr new_gr(self):
+        cdef gr py_val
+        cdef gr_ptr pval
+        pval = gr_heap_init(self.ctx_t)
+        if pval == NULL:
+            raise MemoryError("Failed to allocate memory for gr object")
+        py_val = gr.__new__(gr)
+        py_val.pval = pval
+        py_val.ctx = self
+        py_val._init = True
+        return py_val
+
+    @cython.final
+    cdef inline gr from_str(self, s: str):
+        cdef gr py_val
+        cdef bytes b
+        b = s.encode('utf-8')
+        py_val = self.new_gr()
+        err = gr_set_str(py_val.pval, b, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to parse string")
+        return py_val
+
+    @cython.final
+    cdef inline str to_str(self, val: gr):
+        cdef str py_str
+        cdef char *s
+        err = gr_get_str(&s, val.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to convert to string")
+        py_str = (<bytes>s).decode("utf-8")
+        # XXX: Is this a memory leak?
+        # flint_free(s)
+        return py_str
+
+    @cython.final
+    cdef inline _error(self, int err, str msg):
+        if err & GR_DOMAIN:
+            # return DomainError(msg)
+            return AssertionError(msg)
+        elif err & GR_UNABLE:
+            # return UnableError(msg)
+            return AssertionError(msg)
+        else:
+            return AssertionError("Bad error code")
+
+    @cython.final
+    cdef inline gr _zero(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_zero(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create zero")
+        return res
+
+    @cython.final
+    cdef inline gr _one(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_one(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create one")
+        return res
+
+    @cython.final
+    cdef inline gr _i(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_i(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create sqrt(-1)")
+        return res
+
+    @cython.final
+    cdef inline gr _pos_inf(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_pos_inf(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create +inf")
+        return res
+
+    @cython.final
+    cdef inline gr _neg_inf(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_neg_inf(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create -inf")
+        return res
+
+    @cython.final
+    cdef inline gr _uinf(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_uinf(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create uinf")
+        return res
+
+    @cython.final
+    cdef inline gr _undefined(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_undefined(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create undefined")
+        return res
+
+    @cython.final
+    cdef inline gr _unknown(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_unknown(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create unknown")
+        return res
+
+    @cython.final
+    cdef inline gr _gen(self):
+        cdef int err
+        cdef gr res = self.new_gr()
+        err = gr_gen(res.pval, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot create generator")
+        return res
+
+    @cython.final
+    cdef inline list _gens(self):
+        cdef int err
+        cdef gr g
+        cdef gr_vec_t gens
+        gr_vec_init(gens, 0, self.ctx_t)
+        err = gr_gens(gens, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot get generators")
+        length = gr_vec_length(gens, self.ctx_t)
+        py_gens = [None] * length
+        for 0 <= i < length:
+            g = self.new_gr()
+            err = gr_set(g.pval, gr_vec_entry_ptr(gens, i, self.ctx_t), self.ctx_t)
+            if err != GR_SUCCESS:
+                raise self._error(err, "Failed to copy generator.")
+            py_gens[i] = g
+        gr_vec_clear(gens, self.ctx_t)
+        return py_gens
+
+    @cython.final
+    cdef inline list _gens_recursive(self):
+        cdef int err
+        cdef gr g
+        cdef gr_vec_t gens
+        gr_vec_init(gens, 0, self.ctx_t)
+        err = gr_gens_recursive(gens, self.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Cannot get recursive generators")
+        length = gr_vec_length(gens, self.ctx_t)
+        py_gens = [None] * length
+        for 0 <= i < length:
+            g = self.new_gr()
+            err = gr_set(g.pval, gr_vec_entry_ptr(gens, i, self.ctx_t), self.ctx_t)
+            if err != GR_SUCCESS:
+                raise self._error(err, "Failed to copy generator.")
+            py_gens[i] = g
+        gr_vec_clear(gens, self.ctx_t)
+        return py_gens
+
+
+cdef class gr_scalar_ctx(gr_ctx):
+    pass
+
+
+cdef class gr_poly_ctx(gr_ctx):
+    pass
+
+
+cdef class gr_mpoly_ctx(gr_ctx):
+    pass
+
+
+# cdef class gr_matrix_domain_ctx(gr_ctx):
+#     pass
+
+
+# cdef class gr_matrix_space_ctx(gr_ctx):
+#     pass
+
+
+# cdef class gr_matrix_ring_ctx(gr_ctx):
+#     pass
+
+
+@cython.no_gc
+cdef class _gr_fmpz_ctx(gr_scalar_ctx):
+
+    @staticmethod
+    cdef inline _gr_fmpz_ctx _new():
+        cdef _gr_fmpz_ctx ctx
+        ctx = _gr_fmpz_ctx.__new__(_gr_fmpz_ctx)
+        gr_ctx_init_fmpz(ctx.ctx_t)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class _gr_fmpq_ctx(gr_scalar_ctx):
+
+    @staticmethod
+    cdef inline _gr_fmpq_ctx _new():
+        cdef _gr_fmpq_ctx ctx
+        ctx = _gr_fmpq_ctx.__new__(_gr_fmpq_ctx)
+        gr_ctx_init_fmpq(ctx.ctx_t)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class _gr_fmpzi_ctx(gr_scalar_ctx):
+
+    @staticmethod
+    cdef inline _gr_fmpzi_ctx _new():
+        cdef _gr_fmpzi_ctx ctx
+        ctx = _gr_fmpzi_ctx.__new__(_gr_fmpzi_ctx)
+        gr_ctx_init_fmpzi(ctx.ctx_t)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class _gr_fexpr_ctx(gr_scalar_ctx):
+
+    @staticmethod
+    cdef inline _gr_fexpr_ctx _new():
+        cdef _gr_fexpr_ctx ctx
+        ctx = _gr_fexpr_ctx.__new__(_gr_fexpr_ctx)
+        gr_ctx_init_fexpr(ctx.ctx_t)
+        ctx._init = True
+        return ctx
+
+
+# The global contexts for use in cython code:
+cdef _gr_fmpz_ctx gr_fmpz_ctx_c
+cdef _gr_fmpq_ctx gr_fmpq_ctx_c
+cdef _gr_fmpzi_ctx gr_fmpzi_ctx_c
+cdef _gr_fexpr_ctx gr_fexpr_ctx_c
+
+
+@cython.no_gc
+cdef class gr_nmod_ctx(gr_scalar_ctx):
+    cdef ulong n
+    cdef bint is_field
+
+    @staticmethod
+    cdef inline gr_nmod_ctx _new(ulong n):
+        cdef gr_nmod_ctx ctx
+        cdef bint is_prime = n_is_prime(n)
+        ctx = gr_nmod_ctx.__new__(gr_nmod_ctx)
+        ctx.n = n
+        ctx.is_field = is_prime
+        gr_ctx_init_nmod(ctx.ctx_t, n)
+        ctx._init = True
+        # if is_prime:
+        #     gr_ctx_set_is_field(ctx.ctx_t, T_TRUE)
+        # else:
+        #     gr_ctx_set_is_field(ctx.ctx_t, T_FALSE)
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_fmpz_mod_ctx(gr_scalar_ctx):
+    cdef fmpz_t n
+    cdef bint is_field
+
+    @staticmethod
+    cdef inline gr_fmpz_mod_ctx _new(fmpz n):
+        cdef gr_fmpz_mod_ctx ctx
+        cdef bint is_prime = n.is_prime()
+        ctx = gr_fmpz_mod_ctx.__new__(gr_fmpz_mod_ctx)
+        ctx.is_field = is_prime
+        fmpz_init_set(ctx.n, n.val)
+        gr_ctx_init_fmpz_mod(ctx.ctx_t, ctx.n)
+        ctx._init = True
+        # if is_prime:
+        #     gr_ctx_set_is_field(ctx.ctx_t, T_TRUE)
+        # else:
+        #     gr_ctx_set_is_field(ctx.ctx_t, T_FALSE)
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_fq_ctx(gr_scalar_ctx):
+    cdef fmpz_t p
+    cdef slong d
+
+    @staticmethod
+    cdef inline gr_fq_ctx _new(fmpz p, slong d, char* name):
+        cdef gr_fq_ctx ctx
+        ctx = gr_fq_ctx.__new__(gr_fq_ctx)
+        ctx.d = d
+        fmpz_init_set(ctx.p, p.val)
+        gr_ctx_init_fq(ctx.ctx_t, ctx.p, d, name)
+        ctx._init = True
+        # XXX: free name_c?
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_fq_nmod_ctx(gr_scalar_ctx):
+    cdef ulong p
+    cdef slong d
+
+    @staticmethod
+    cdef inline gr_fq_nmod_ctx _new(ulong p, slong d, char* name):
+        cdef gr_fq_nmod_ctx ctx
+        ctx = gr_fq_nmod_ctx.__new__(gr_fq_nmod_ctx)
+        ctx.p = p
+        ctx.d = d
+        gr_ctx_init_fq_nmod(ctx.ctx_t, p, d, name)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_fq_zech_ctx(gr_scalar_ctx):
+    cdef ulong p
+    cdef slong d
+
+    @staticmethod
+    cdef inline gr_fq_zech_ctx _new(ulong p, slong d, char* name):
+        cdef gr_fq_zech_ctx ctx
+        ctx = gr_fq_zech_ctx.__new__(gr_fq_zech_ctx)
+        ctx.p = p
+        ctx.d = d
+        gr_ctx_init_fq_zech(ctx.ctx_t, p, d, name)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_nf_ctx(gr_scalar_ctx):
+    cdef fmpq_poly_t poly
+
+    @staticmethod
+    cdef inline gr_nf_ctx _new(fmpq_poly poly):
+        cdef gr_nf_ctx ctx
+        ctx = gr_nf_ctx.__new__(gr_nf_ctx)
+        fmpq_poly_init(ctx.poly)
+        fmpq_poly_set(ctx.poly, poly.val)
+        gr_ctx_init_nf(ctx.ctx_t, ctx.poly)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_nf_fmpz_poly_ctx(gr_scalar_ctx):
+    cdef fmpz_poly_t poly
+
+    @staticmethod
+    cdef inline gr_nf_fmpz_poly_ctx _new(fmpz_poly poly):
+        cdef gr_nf_fmpz_poly_ctx ctx
+        ctx = gr_nf_fmpz_poly_ctx.__new__(gr_nf_fmpz_poly_ctx)
+        fmpz_poly_init(ctx.poly)
+        fmpz_poly_set(ctx.poly, poly.val)
+        gr_ctx_init_nf_fmpz_poly(ctx.ctx_t, ctx.poly)
+        ctx._init = True
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_real_qqbar_ctx(gr_scalar_ctx):
+    cdef slong deg_limit
+    cdef slong bits_limit
+
+    @staticmethod
+    cdef inline gr_real_qqbar_ctx _new(slong deg_limit, slong bits_limit):
+        cdef gr_real_qqbar_ctx ctx
+        ctx = gr_real_qqbar_ctx.__new__(gr_real_qqbar_ctx)
+        gr_ctx_init_real_qqbar(ctx.ctx_t)
+        ctx._init = True
+        ctx.deg_limit = deg_limit
+        ctx.bits_limit = bits_limit
+        if deg_limit != -1 or bits_limit != -1:
+            # Maybe add setters for these?
+            _gr_ctx_qqbar_set_limits(ctx.ctx_t, deg_limit, bits_limit)
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_complex_qqbar_ctx(gr_scalar_ctx):
+    cdef slong deg_limit
+    cdef slong bits_limit
+
+    @staticmethod
+    cdef inline gr_complex_qqbar_ctx _new(slong deg_limit, slong bits_limit):
+        cdef gr_complex_qqbar_ctx ctx
+        ctx = gr_complex_qqbar_ctx.__new__(gr_complex_qqbar_ctx)
+        gr_ctx_init_complex_qqbar(ctx.ctx_t)
+        ctx._init = True
+        ctx.deg_limit = deg_limit
+        ctx.bits_limit = bits_limit
+        if deg_limit != -1 or bits_limit != -1:
+            # Maybe add setters for these?
+            _gr_ctx_qqbar_set_limits(ctx.ctx_t, deg_limit, bits_limit)
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_real_ca_ctx(gr_scalar_ctx):
+    cdef dict options
+
+    @staticmethod
+    cdef inline gr_real_ca_ctx _new(dict options):
+        cdef gr_real_ca_ctx ctx
+        if options:
+            raise NotImplementedError("Options not implemented yet")
+        ctx = gr_real_ca_ctx.__new__(gr_real_ca_ctx)
+        gr_ctx_init_real_ca(ctx.ctx_t)
+        ctx._init = True
+        ctx.options = options.copy()
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_complex_ca_ctx(gr_scalar_ctx):
+    cdef dict options
+
+    @staticmethod
+    cdef inline gr_complex_ca_ctx _new(dict options):
+        cdef gr_complex_ca_ctx ctx
+        if options:
+            raise NotImplementedError("Options not implemented yet")
+        ctx = gr_complex_ca_ctx.__new__(gr_complex_ca_ctx)
+        gr_ctx_init_complex_ca(ctx.ctx_t)
+        ctx._init = True
+        ctx.options = options.copy()
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_real_algebraic_ca_ctx(gr_scalar_ctx):
+    cdef dict options
+
+    @staticmethod
+    cdef inline gr_real_algebraic_ca_ctx _new(dict options):
+        cdef gr_real_algebraic_ca_ctx ctx
+        if options:
+            raise NotImplementedError("Options not implemented yet")
+        ctx = gr_real_algebraic_ca_ctx.__new__(gr_real_algebraic_ca_ctx)
+        gr_ctx_init_real_algebraic_ca(ctx.ctx_t)
+        ctx._init = True
+        ctx.options = options.copy()
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_complex_algebraic_ca_ctx(gr_scalar_ctx):
+    cdef dict options
+
+    @staticmethod
+    cdef inline gr_complex_algebraic_ca_ctx _new(dict options):
+        cdef gr_complex_algebraic_ca_ctx ctx
+        if options:
+            raise NotImplementedError("Options not implemented yet")
+        ctx = gr_complex_algebraic_ca_ctx.__new__(gr_complex_algebraic_ca_ctx)
+        gr_ctx_init_complex_algebraic_ca(ctx.ctx_t)
+        ctx._init = True
+        ctx.options = options.copy()
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_complex_extended_ca_ctx(gr_scalar_ctx):
+    cdef dict options
+
+    @staticmethod
+    cdef inline gr_complex_extended_ca_ctx _new(dict options):
+        cdef gr_complex_extended_ca_ctx ctx
+        if options:
+            raise NotImplementedError("Options not implemented yet")
+        ctx = gr_complex_extended_ca_ctx.__new__(gr_complex_extended_ca_ctx)
+        gr_ctx_init_complex_extended_ca(ctx.ctx_t)
+        ctx._init = True
+        ctx.options = options.copy()
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_real_float_arf_ctx(gr_scalar_ctx):
+    cdef slong prec
+
+    @staticmethod
+    cdef inline gr_real_float_arf_ctx _new(slong prec):
+        cdef gr_real_float_arf_ctx ctx
+        ctx = gr_real_float_arf_ctx.__new__(gr_real_float_arf_ctx)
+        gr_ctx_init_real_float_arf(ctx.ctx_t, prec)
+        ctx._init = True
+        ctx.prec = prec
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_complex_float_acf_ctx(gr_scalar_ctx):
+    cdef slong prec
+
+    @staticmethod
+    cdef inline gr_complex_float_acf_ctx _new(slong prec):
+        cdef gr_complex_float_acf_ctx ctx
+        ctx = gr_complex_float_acf_ctx.__new__(gr_complex_float_acf_ctx)
+        gr_ctx_init_complex_float_acf(ctx.ctx_t, prec)
+        ctx._init = True
+        ctx.prec = prec
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_real_arb_ctx(gr_scalar_ctx):
+    cdef slong prec
+
+    @staticmethod
+    cdef inline gr_real_arb_ctx _new(slong prec):
+        cdef gr_real_arb_ctx ctx
+        ctx = gr_real_arb_ctx.__new__(gr_real_arb_ctx)
+        gr_ctx_init_real_arb(ctx.ctx_t, prec)
+        ctx._init = True
+        ctx.prec = prec
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_complex_acb_ctx(gr_scalar_ctx):
+    cdef slong prec
+
+    @staticmethod
+    cdef inline gr_complex_acb_ctx _new(slong prec):
+        cdef gr_complex_acb_ctx ctx
+        ctx = gr_complex_acb_ctx.__new__(gr_complex_acb_ctx)
+        gr_ctx_init_complex_acb(ctx.ctx_t, prec)
+        ctx._init = True
+        ctx.prec = prec
+        return ctx
+
+
+# @cython.no_gc
+# cdef class _gr_fmpz_poly_ctx(gr_poly_ctx):
+#
+#     @staticmethod
+#     cdef _gr_fmpz_poly_ctx _new()
+
+
+# @cython.no_gc
+# cdef class _gr_fmpq_poly_ctx(gr_poly_ctx):
+#
+#     @staticmethod
+#     cdef _gr_fmpq_poly_ctx _new()
+
+
+# @cython.no_gc
+# cdef class _gr_gr_poly_ctx(gr_poly_ctx):
+#
+#     @staticmethod
+#     cdef _gr_gr_poly_ctx _new()
+
+
+@cython.no_gc
+cdef class gr_gr_poly_ctx(gr_poly_ctx):
+    cdef gr_ctx base_ctx
+
+    @staticmethod
+    cdef inline gr_gr_poly_ctx _new(gr_ctx base_ctx):
+        cdef gr_gr_poly_ctx ctx
+        ctx = gr_gr_poly_ctx.__new__(gr_gr_poly_ctx)
+        gr_ctx_init_gr_poly(ctx.ctx_t, base_ctx.ctx_t)
+        ctx._init = True
+        ctx.base_ctx = base_ctx
+        return ctx
+
+
+@cython.no_gc
+cdef class gr_gr_mpoly_ctx(gr_mpoly_ctx):
+    cdef gr_ctx base_ctx
+    cdef ordering_t _order
+    cdef tuple _names
+    cdef slong _nvars
+
+    @staticmethod
+    cdef inline gr_gr_mpoly_ctx _new(gr_ctx base_ctx, tuple names, order):
+        cdef gr_gr_mpoly_ctx ctx
+        cdef ordering_t ord_c
+        cdef slong nvars
+        cdef const char **names_c
+        cdef int status
+
+        nvars = len(names)
+        names_b = [name.encode('utf-8') for name in names]
+
+        ord_c = ordering_py_to_c(order)
+        ctx = gr_gr_mpoly_ctx.__new__(gr_gr_mpoly_ctx)
+        gr_ctx_init_gr_mpoly(ctx.ctx_t, base_ctx.ctx_t, nvars, ord_c)
+        ctx._init = True
+        ctx.base_ctx = base_ctx
+        ctx._order = ord_c
+        ctx._nvars = nvars
+        ctx._names = names
+
+        names_c = <const char **>libc.stdlib.malloc(nvars * sizeof(const char *))
+        if names_c == NULL:
+            raise MemoryError("Failed to allocate memory for generator names")
+        try:
+            for i in range(nvars):
+                names_c[i] = names_b[i]
+            status = gr_ctx_set_gen_names(ctx.ctx_t, names_c)
+        finally:
+            libc.stdlib.free(names_c)
+        if status != GR_SUCCESS:
+            raise MemoryError("Failed to set generator names")
+
+        return ctx
+
+
+# @cython.no_gc
+# cdef class gr_fmpz_mpoly_q_ctx(gr_mpoly_ctx):
+#     cdef ordering_t _order
+#     cdef slong _nvars
+#
+#     @staticmethod
+#     cdef inline gr_fmpz_mpoly_q_ctx _new(slong nvars, order):
+#         cdef gr_fmpz_mpoly_q_ctx ctx
+#         cdef ordering_t ord_c
+#         ord_c = ordering_py_to_c(order)
+#         ctx = gr_fmpz_mpoly_q_ctx.__new__(gr_fmpz_mpoly_q_ctx)
+#         gr_ctx_init_fmpz_mpoly_q(ctx.ctx_t, nvars, ord_c)
+#         ctx._init = True
+#         ctx._order = ord_c
+#         ctx._nvars = nvars
+#         return ctx
+
+
+# @cython.no_gc
+# cdef class gr_series_mod_gr_poly_ctx(gr_ctx):
+#     cdef gr_ctx base_ctx
+#     cdef slong _n
+#
+#     @staticmethod
+#     cdef inline gr_series_mod_gr_poly_ctx _new(gr_ctx base_ctx, slong n):
+#         cdef gr_series_mod_gr_poly_ctx ctx
+#         ctx = gr_series_mod_gr_poly_ctx.__new__(gr_series_mod_gr_poly_ctx)
+#         gr_ctx_init_series_mod_gr_poly(ctx.ctx_t, base_ctx.ctx_t, n)
+#         ctx._init = True
+#         ctx.base_ctx = base_ctx
+#         ctx._n = n
+#         return ctx
+
+
+@cython.no_gc
+cdef class gr_series_ctx(gr_ctx):
+    cdef gr_ctx base_ctx
+    cdef slong _prec
+
+    @staticmethod
+    cdef inline gr_series_ctx _new(gr_ctx base_ctx, slong prec):
+        cdef gr_series_ctx ctx
+        ctx = gr_series_ctx.__new__(gr_series_ctx)
+        gr_ctx_init_gr_series(ctx.ctx_t, base_ctx.ctx_t, prec)
+        ctx._init = True
+        ctx.base_ctx = base_ctx
+        ctx._prec = prec
+        return ctx
+
+
+@cython.no_gc
+cdef class gr(flint_scalar):
+    cdef gr_ptr pval
+    cdef gr_ctx ctx
+    cdef bint _init
+
+    @cython.final
+    cdef inline _error(self, int err, str msg):
+        return self.ctx._error(err, msg)
+
+    @cython.final
+    cdef inline truth_t _equal(self, other: gr):
+        return gr_equal(self.pval, other.pval, self.ctx.ctx_t)
+
+    @cython.final
+    cdef inline truth_t _is_zero(self):
+        return gr_is_zero(self.pval, self.ctx.ctx_t)
+
+    @cython.final
+    cdef inline truth_t _is_one(self):
+        return gr_is_one(self.pval, self.ctx.ctx_t)
+
+    @cython.final
+    cdef inline truth_t _is_neg_one(self):
+        return gr_is_neg_one(self.pval, self.ctx.ctx_t)
+
+    # @cython.final
+    # cdef inline truth_t _is_integer(self):
+    #     return gr_is_integer(self.pval, self.ctx.ctx_t)
+
+    # @cython.final
+    # cdef inline truth_t _is_rational(self):
+    #     return gr_is_rational(self.pval, self.ctx.ctx_t)
+
+    @cython.final
+    cdef inline gr _neg(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_neg(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to negate gr object")
+        return res
+
+    @cython.final
+    cdef inline gr _add(self, other: gr):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_add(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to add gr objects")
+        return res
+
+    @cython.final
+    cdef inline gr _add_si(self, other: slong):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_add_si(res.pval, self.pval, other, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to add gr object and int")
+        return res
+
+    @cython.final
+    cdef inline gr _sub(self, other: gr):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_sub(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to subtract gr objects")
+        return res
+
+    @cython.final
+    cdef inline gr _sub_si(self, other: int):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_sub_si(res.pval, self.pval, other, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to subtract gr object and int")
+        return res
+
+    @cython.final
+    cdef inline gr _mul(self, other: gr):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_mul(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to multiply gr objects")
+        return res
+
+    @cython.final
+    cdef inline gr _mul_si(self, other: slong):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_mul_si(res.pval, self.pval, other, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to multiply gr object and int")
+        return res
+
+    @cython.final
+    cdef inline gr _inv(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_inv(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to invert gr object")
+        return res
+
+    @cython.final
+    cdef inline gr _div(self, other: gr):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_div(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to divide gr objects")
+        return res
+
+    @cython.final
+    cdef inline gr _div_si(self, other: slong):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_div_si(res.pval, self.pval, other, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to divide gr object and int")
+        return res
+
+    @cython.final
+    cdef inline gr _pow_si(self, other: slong):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_pow_si(res.pval, self.pval, other, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to raise gr object to an integer power")
+        return res
+
+    @cython.final
+    cdef inline truth_t _is_square(self):
+        return gr_is_square(self.pval, self.ctx.ctx_t)
+
+    @cython.final
+    cdef inline gr _sqrt(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_sqrt(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute square root of gr object")
+        return res
+
+    @cython.final
+    cdef inline gr _rsqrt(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_rsqrt(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute reciprocal square root of gr object")
+        return res
+
+    @cython.final
+    cdef inline gr _gcd(self, other: gr):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_gcd(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute gcd of gr objects")
+        return res
+
+    @cython.final
+    cdef inline gr _lcm(self, other: gr):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_lcm(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute lcm of gr objects")
+        return res
+
+    cdef inline _factor(self):
+        cdef int err, i
+        cdef slong length, exp_s
+        cdef gr c, f
+        cdef gr_ptr fac_i, exp_i
+        cdef gr_vec_t factors, exponents
+        cdef int flags = 0  # XXX: What is flags?
+        c = self.ctx.new_gr()
+        gr_vec_init(factors, 0, self.ctx.ctx_t)
+        gr_vec_init(exponents, 0, gr_fmpz_ctx_c.ctx_t)
+        err = gr_factor(c.pval, factors, exponents, self.pval, flags, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to factor gr object")
+        length = gr_vec_length(factors, self.ctx.ctx_t)
+        py_factors = [None] * length
+        for 0 <= i < length:
+            f = self.ctx.new_gr()
+            fac_i = gr_vec_entry_ptr(factors, i, self.ctx.ctx_t)
+            err = gr_set(f.pval, fac_i, self.ctx.ctx_t)
+            if err != GR_SUCCESS:
+                raise self._error(err, "Failed to copy factor.")
+            exp_i = gr_vec_entry_ptr(exponents, i, gr_fmpz_ctx_c.ctx_t)
+            err = gr_get_si(&exp_s, exp_i, gr_fmpz_ctx_c.ctx_t)
+            if err != GR_SUCCESS:
+                raise self._error(err, "Failed to get integer value of exponent.")
+            exp = exp_s
+            py_factors[i] = (f, exp)
+        gr_vec_clear(factors, self.ctx.ctx_t)
+        gr_vec_clear(exponents, gr_fmpz_ctx_c.ctx_t)
+        return c, py_factors
+
+    cdef inline gr _numerator(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_numerator(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute numerator of gr object")
+        return res
+
+    cdef inline gr _denominator(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_denominator(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute denominator of gr object")
+        return res
+
+    cdef inline gr _floor(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_floor(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute floor of gr object")
+        return res
+
+    cdef inline gr _ceil(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_ceil(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute ceil of gr object")
+        return res
+
+    cdef inline gr _trunc(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_trunc(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute trunc of gr object")
+        return res
+
+    cdef inline gr _nint(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_nint(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute nint of gr object")
+        return res
+
+    cdef inline gr _abs(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_abs(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute absolute value of gr object")
+        return res
+
+    cdef inline gr _conj(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_conj(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute conjugate of gr object")
+        return res
+
+    cdef inline gr _re(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_re(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute real part of gr object")
+        return res
+
+    cdef inline gr _im(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_im(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute imaginary part of gr object")
+        return res
+
+    cdef inline gr _sgn(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_sgn(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute sign of gr object")
+        return res
+
+    cdef inline gr _csgn(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_csgn(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute complex sign of gr object")
+        return res
+
+    cdef inline gr _arg(self):
+        cdef int err
+        cdef gr res = self.ctx.new_gr()
+        err = gr_arg(res.pval, self.pval, self.ctx.ctx_t)
+        if err != GR_SUCCESS:
+            raise self._error(err, "Failed to compute argument of gr object")
+        return res
+
+    # cdef inline truth_t _le(self, other: gr):
+    #     return gr_le(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _lt(self, other: gr):
+    #     return gr_lt(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _ge(self, other: gr):
+    #     return gr_ge(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _gt(self, other: gr):
+    #     return gr_gt(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _abs_le(self, other: gr):
+    #    return gr_abs_le(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _abs_lt(self, other: gr):
+    #     return gr_abs_lt(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _abs_ge(self, other: gr):
+    #     return gr_abs_ge(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline truth_t _abs_gt(self, other: gr):
+    #     return gr_abs_gt(self.pval, other.pval, self.ctx.ctx_t)
+
+    # cdef inline gr _min(self, other: gr):
+    #     cdef int err
+    #     cdef gr res = self.ctx.new_gr()
+    #     err = gr_min(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+    #     if err != GR_SUCCESS:
+    #         raise self._error(err, "Failed to compute minimum of gr objects")
+    #     return res
+
+    # cdef inline gr _max(self, other: gr):
+    #     cdef int err
+    #     cdef gr res = self.ctx.new_gr()
+    #     err = gr_max(res.pval, self.pval, other.pval, self.ctx.ctx_t)
+    #     if err != GR_SUCCESS:
+    #         raise self._error(err, "Failed to compute maximum of gr objects")
+    #     return res
diff --git a/src/flint/types/_gr.pyx b/src/flint/types/_gr.pyx
new file mode 100644
index 00000000..ff9859ce
--- /dev/null
+++ b/src/flint/types/_gr.pyx
@@ -0,0 +1,1605 @@
+cimport cython
+
+from flint.flintlib.types.flint cimport slong
+from flint.flintlib.types.gr cimport GR_SUCCESS
+from flint.flintlib.functions.gr cimport (
+    gr_ctx_clear,
+    gr_heap_clear,
+
+    gr_set_si,
+)
+from flint.flintlib.functions.gr_domains cimport (
+    gr_ctx_is_finite,
+    gr_ctx_is_multiplicative_group,
+    gr_ctx_is_ring,
+    gr_ctx_is_commutative_ring,
+    gr_ctx_is_integral_domain,
+    gr_ctx_is_unique_factorization_domain,
+    gr_ctx_is_field,
+    gr_ctx_is_algebraically_closed,
+    gr_ctx_is_finite_characteristic,
+    gr_ctx_is_ordered_ring,
+    gr_ctx_is_zero_ring,
+    gr_ctx_is_exact,
+    gr_ctx_is_canonical,
+    gr_ctx_has_real_prec,
+    gr_ctx_set_real_prec,
+    gr_ctx_get_real_prec,
+)
+
+from flint.flint_base.flint_base cimport ordering_c_to_py
+from flint.flint_base.flint_base import Ordering
+
+
+@cython.no_gc
+cdef class gr_ctx(flint_ctx):
+    """Base class for all gr contexts.
+
+    This class should not be instantiated directly. Instead, use one of the
+    derived classes, e.g. gr_nmod_ctx or gr_fmpz_ctx.
+    """
+    def __init__(self, *args, **kwargs):
+        raise TypeError("Cannot create gr_ctx object directly."
+                        " Use e.g. gr_nmod_ctx.new(n) instead.")
+
+    def new(self):
+        """Create a context of a given type.
+
+        This method should be called on a derived class:
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> ctx10 = gr_real_float_arf_ctx.new(10)
+        >>> ctx10
+        gr_real_float_arf_ctx(10)
+        >>> ctx10.real_prec
+        10
+        >>> ctx10(2).sqrt()  # 10 bits
+        1.414
+        >>> ctx100 = gr_real_float_arf_ctx.new(100)
+        >>> ctx100(2).sqrt()  # 100 bits
+        1.414213562373095048801688724209
+        """
+        raise NotImplementedError("Cannot create gr_ctx object directly."
+                                  " Use e.g. gr_nmod_ctx.new(n) instead.")
+
+    def __dealloc__(self):
+        if self._init:
+            self._init = False
+            gr_ctx_clear(self.ctx_t)
+
+    @property
+    def is_finite(self) -> bool | None:
+        """True if the domain is finite (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_fmpz_ctx, gr_nmod_ctx
+        >>> R1 = gr_fmpz_ctx
+        >>> R2 = gr_nmod_ctx.new(5)
+        >>> R1.is_finite, R2.is_finite
+        (False, True)
+        """
+        return truth_to_py(gr_ctx_is_finite(self.ctx_t))
+
+    @property
+    def is_multiplicative_group(self) -> bool | None:
+        """True if the domain is a multiplicative group (can be ``None`` if unknown)."""
+        return truth_to_py(gr_ctx_is_multiplicative_group(self.ctx_t))
+
+    @property
+    def is_ring(self) -> bool | None:
+        """True if the domain is a ring (can be ``None`` if unknown)."""
+        return truth_to_py(gr_ctx_is_ring(self.ctx_t))
+
+    @property
+    def is_commutative_ring(self) -> bool | None:
+        """True if the domain is a commutative ring (can be ``None`` if unknown)."""
+        return truth_to_py(gr_ctx_is_commutative_ring(self.ctx_t))
+
+    @property
+    def is_integral_domain(self) -> bool | None:
+        """True if the domain is an integral domain (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_nmod_ctx
+        >>> R1 = gr_nmod_ctx.new(5)
+        >>> R2 = gr_nmod_ctx.new(6)
+        >>> R1.is_integral_domain, R2.is_integral_domain
+        (True, False)
+        """
+        return truth_to_py(gr_ctx_is_integral_domain(self.ctx_t))
+
+    @property
+    def is_unique_factorization_domain(self) -> bool | None:
+        """True if the domain is a unique factorization domain (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_fmpz_ctx
+        >>> R1 = gr_fmpz_ctx
+        >>> R1.is_unique_factorization_domain
+        True
+        """
+        return truth_to_py(gr_ctx_is_unique_factorization_domain(self.ctx_t))
+
+    @property
+    def is_field(self) -> bool | None:
+        """True if the domain is a field (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_fmpz_ctx, gr_fmpq_ctx
+        >>> R1 = gr_fmpz_ctx
+        >>> R2 = gr_fmpq_ctx
+        >>> R1.is_field, R2.is_field
+        (False, True)
+        """
+        return truth_to_py(gr_ctx_is_field(self.ctx_t))
+
+    @property
+    def is_algebraically_closed(self) -> bool | None:
+        """True if the domain is algebraically closed (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_complex_qqbar_ctx, gr_real_qqbar_ctx
+        >>> R1 = gr_complex_qqbar_ctx.new()
+        >>> R2 = gr_real_qqbar_ctx.new()
+        >>> R1.is_algebraically_closed, R2.is_algebraically_closed
+        (True, False)
+        """
+        return truth_to_py(gr_ctx_is_algebraically_closed(self.ctx_t))
+
+    @property
+    def is_finite_characteristic(self) -> bool | None:
+        """True if the domain has finite characteristic (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_nmod_ctx, gr_fmpz_ctx
+        >>> R1 = gr_nmod_ctx.new(5)
+        >>> R2 = gr_fmpz_ctx
+        >>> R1.is_finite_characteristic, R2.is_finite_characteristic
+        (True, False)
+        """
+        return truth_to_py(gr_ctx_is_finite_characteristic(self.ctx_t))
+
+    @property
+    def is_ordered_ring(self) -> bool | None:
+        """True if the domain is an ordered ring (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_fmpz_ctx, gr_nmod_ctx
+        >>> R1 = gr_fmpz_ctx
+        >>> R2 = gr_nmod_ctx.new(5)
+        >>> R1.is_ordered_ring, R2.is_ordered_ring
+        (True, None)
+        """
+        return truth_to_py(gr_ctx_is_ordered_ring(self.ctx_t))
+
+    @property
+    def is_zero_ring(self) -> bool | None:
+        """True if the domain is a zero ring (can be ``None`` if unknown)."""
+        return truth_to_py(gr_ctx_is_zero_ring(self.ctx_t))
+
+    @property
+    def is_exact(self) -> bool | None:
+        """True if the domain is exact (can be ``None`` if unknown).
+
+        >>> from flint.types._gr import gr_fmpz_ctx, gr_real_float_arf_ctx
+        >>> R1 = gr_fmpz_ctx
+        >>> R2 = gr_real_float_arf_ctx.new(10)
+        >>> R1.is_exact, R2.is_exact
+        (True, False)
+        """
+        return truth_to_py(gr_ctx_is_exact(self.ctx_t))
+
+    @property
+    def is_canonical(self) -> bool | None:
+        """True if elements of the domain are always canonical."""
+        return truth_to_py(gr_ctx_is_canonical(self.ctx_t))
+
+    @property
+    def has_real_prec(self) -> bool | None:
+        """True if the domain has real precision (can be ``None`` if unknown).
+
+        See also :attr:`real_prec`.
+        """
+        return truth_to_py(gr_ctx_has_real_prec(self.ctx_t))
+
+    @property
+    def real_prec(self) -> int:
+        """Real precision of the domain.
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> ctx = gr_real_float_arf_ctx.new(10)
+        >>> ctx.real_prec
+        10
+        >>> ctx.real_prec = 20
+        >>> ctx.real_prec
+        20
+        """
+        cdef slong prec
+        if gr_ctx_has_real_prec(self.ctx_t) != T_TRUE:
+            raise ValueError("Real precision is not available")
+        if gr_ctx_get_real_prec(&prec, self.ctx_t) != GR_SUCCESS:
+            raise ValueError("Failed to get real precision")
+        return prec
+
+    @real_prec.setter
+    def real_prec(self, prec: slong) -> None:
+        if gr_ctx_has_real_prec(self.ctx_t) != T_TRUE:
+            raise ValueError("Real precision is not available")
+        if gr_ctx_set_real_prec(self.ctx_t, prec) != GR_SUCCESS:
+            raise ValueError("Failed to set real precision")
+
+    def __call__(self, arg) -> gr:
+        """Create a new element of the domain.
+
+        >>> from flint.types._gr import gr_fmpz_ctx
+        >>> ctx = gr_fmpz_ctx
+        >>> ctx(2)
+        2
+        """
+        return self.from_str(str(arg))
+
+    def zero(self) -> gr:
+        """Return the zero element of the domain.
+
+        >>> from flint.types._gr import gr_fmpz_ctx
+        >>> ctx = gr_fmpz_ctx
+        >>> ctx.zero()
+        0
+        """
+        return self._zero()
+
+    def one(self) -> gr:
+        """Return the one element of the domain.
+
+        >>> from flint.types._gr import gr_fmpz_ctx
+        >>> ctx = gr_fmpz_ctx
+        >>> ctx.one()
+        1
+        """
+        return self._one()
+
+    def i(self) -> gr:
+        """Return the imaginary unit of the domain (if available).
+
+        >>> from flint.types._gr import gr_fmpzi_ctx
+        >>> ctx = gr_fmpzi_ctx
+        >>> ctx.i()
+        I
+        """
+        return self._i()
+
+    def pos_inf(self) -> gr:
+        """Return the positive infinity element of the domain (if available).
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> ctx = gr_real_float_arf_ctx.new(10)
+        >>> ctx.pos_inf()
+        +inf
+        """
+        return self._pos_inf()
+
+    def neg_inf(self) -> gr:
+        """Return the negative infinity element of the domain (if available).
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> ctx = gr_real_float_arf_ctx.new(10)
+        >>> ctx.neg_inf()
+        -inf
+        """
+        return self._neg_inf()
+
+    def uinf(self) -> gr:
+        """Return the unsigned infinity element of the domain (if available).
+
+        >>> from flint.types._gr import gr_complex_extended_ca_ctx
+        >>> ctx = gr_complex_extended_ca_ctx.new()
+        >>> ctx.uinf()
+        UnsignedInfinity
+        """
+        return self._uinf()
+
+    def undefined(self) -> gr:
+        """Return the undefined element of the domain (if available).
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> ctx = gr_real_float_arf_ctx.new(10)
+        >>> ctx.undefined()
+        nan
+        """
+        return self._undefined()
+
+    def unknown(self) -> gr:
+        """Return the unknown element of the domain (if available).
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> ctx = gr_real_float_arf_ctx.new(10)
+        >>> ctx.unknown()
+        nan
+        """
+        return self._unknown()
+
+    def gen(self) -> gr:
+        """Return the generator of the domain (if available).
+
+        >>> from flint.types._gr import gr_fmpzi_ctx, gr_fq_nmod_ctx
+        >>> ctx = gr_fmpzi_ctx
+        >>> ctx.gen()
+        I
+        >>> ctx = gr_fq_nmod_ctx.new(5, 2)
+        >>> ctx.gen()
+        a
+        """
+        return self._gen()
+
+    def gens(self) -> list[gr]:
+        """Return the top-level generators of the domain
+
+        >>> from flint.types._gr import gr_fmpzi_ctx, gr_gr_mpoly_ctx
+        >>> ctx = gr_gr_mpoly_ctx.new(gr_fmpzi_ctx, ['x', 'y'])
+        >>> ctx.gens()
+        [x, y]
+        >>> gr_fmpzi_ctx.gens()
+        [I]
+        >>> ctx.gens_recursive()
+        [I, x, y]
+        """
+        return self._gens()
+
+    def gens_recursive(self) -> list[gr]:
+        """Return all generators of the domain
+
+        See :meth:`gens` for an example.
+        """
+        return self._gens_recursive()
+
+
+cdef class gr_scalar_ctx(gr_ctx):
+    """Base class for all scalar contexts."""
+    pass
+
+
+cdef class gr_poly_ctx(gr_ctx):
+    """Base class for dense univariate polynomial contexts."""
+    pass
+
+
+cdef class gr_mpoly_ctx(gr_ctx):
+    """Base class for dense multivariate polynomial contexts."""
+    pass
+
+
+# cdef class gr_matrix_domain_ctx(gr_ctx):
+#     pass
+
+
+# cdef class gr_matrix_space_ctx(gr_ctx):
+#     pass
+
+
+# cdef class gr_matrix_ring_ctx(gr_ctx):
+#     pass
+
+
+# Global contexts for Cython code
+cdef _gr_fmpz_ctx gr_fmpz_ctx_c = _gr_fmpz_ctx._new()
+cdef _gr_fmpq_ctx gr_fmpq_ctx_c = _gr_fmpq_ctx._new()
+cdef _gr_fmpzi_ctx gr_fmpzi_ctx_c = _gr_fmpzi_ctx._new()
+cdef _gr_fexpr_ctx gr_fexpr_ctx_c = _gr_fexpr_ctx._new()
+
+# Global contexts for Python code
+gr_fmpz_ctx = gr_fmpz_ctx_c
+gr_fmpq_ctx = gr_fmpq_ctx_c
+gr_fmpzi_ctx = gr_fmpzi_ctx_c
+gr_fexpr_ctx = gr_fexpr_ctx_c
+
+
+@cython.no_gc
+cdef class _gr_fmpz_ctx(gr_scalar_ctx):
+    r"""Context for arbitrary precision integers, `\mathbb{Z}`."""
+
+    @staticmethod
+    def new() -> _gr_fmpz_ctx:
+        """The global context for arbitrary precision integers.
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z
+        gr_fmpz_ctx
+        >>> Z(10)
+        10
+        >>> Z(10) + Z(20)
+        30
+        >>> Z(10) * Z(20)
+        200
+        >>> Z(10) / Z(20)
+        Traceback (most recent call last):
+            ...
+        AssertionError: Failed to divide gr objects
+        >>> Z(10).factor()
+        (1, [(2, 1), (5, 1)])
+        >>> Z(9).sqrt()
+        3
+        """
+        return gr_fmpz_ctx_c
+
+    def __repr__(self):
+        return "gr_fmpz_ctx"
+
+
+@cython.no_gc
+cdef class _gr_fmpq_ctx(gr_scalar_ctx):
+    r"""Context for arbitrary precision rationals, `\mathbb{Q}`."""
+
+    @staticmethod
+    def new() -> _gr_fmpq_ctx:
+        """The global context for arbitrary precision rationals.
+
+        >>> from flint.types._gr import gr_fmpq_ctx as Q
+        >>> Q
+        gr_fmpq_ctx
+        >>> Q(2) / 3
+        2/3
+        """
+        return gr_fmpq_ctx_c
+
+    def __repr__(self):
+        return "gr_fmpq_ctx"
+
+
+@cython.no_gc
+cdef class _gr_fmpzi_ctx(gr_scalar_ctx):
+    r"""Context for Gaussian integers, `\mathbb{Z}[i]`. """
+
+    @staticmethod
+    def new() -> _gr_fmpzi_ctx:
+        """The global context for Gaussian integers.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx as ZI
+        >>> ZI
+        gr_fmpzi_ctx
+        >>> ZI.i()
+        I
+        >>> ZI.gen()
+        I
+        >>> I = ZI.i()
+        >>> (2 + 3*I) ** 2
+        (-5+12*I)
+        """
+        return gr_fmpzi_ctx_c
+
+    def __repr__(self):
+        return "gr_fmpzi_ctx"
+
+
+@cython.no_gc
+cdef class _gr_fexpr_ctx(gr_scalar_ctx):
+    """Context for symbolic expressions. """
+
+    @staticmethod
+    def new() -> _gr_fexpr_ctx:
+        """The global context for symbolic expressions.
+
+        >>> from flint.types._gr import gr_fexpr_ctx as EX
+        >>> EX
+        gr_fexpr_ctx
+        >>> (EX(1) + EX(2)) / 3
+        Div(Add(1, 2), 3)
+        """
+        return gr_fexpr_ctx_c
+
+    def __repr__(self):
+        return "gr_fexpr_ctx"
+
+
+@cython.no_gc
+cdef class gr_nmod_ctx(gr_scalar_ctx):
+    r"""Context for integers modulo `n`, `\mathbb{Z}/n\mathbb{Z}`.
+
+    The modulus `n` must be word-sized. Use :class:`gr_fmpz_mod_ctx` for
+    arbitrary precision modulus.
+    """
+    @staticmethod
+    def new(ulong n) -> gr_nmod_ctx:
+        """Create a new context for integers modulo `n`.
+
+        >>> from flint.types._gr import gr_nmod_ctx
+        >>> Z5 = gr_nmod_ctx.new(5)
+        >>> Z5
+        gr_nmod_ctx(5)
+        >>> Z5(2) + Z5(3)
+        0
+        >>> Z5.modulus()
+        5
+        >>> Z5.is_prime
+        True
+        """
+        return gr_nmod_ctx._new(n)
+
+    @property
+    def is_prime(self) -> bool:
+        return self.is_field
+
+    def modulus(self):
+        return self.n
+
+    def __repr__(self):
+        return f"gr_nmod_ctx({self.n})"
+
+
+@cython.no_gc
+cdef class gr_fmpz_mod_ctx(gr_scalar_ctx):
+    r"""Context for integers modulo `n`, `\mathbb{Z}/n\mathbb{Z}`.
+
+    The modulus `n` can be arbitrary precision. See also :class:`gr_nmod_ctx`
+    for word-sized modulus.
+    """
+
+    @staticmethod
+    def new(n) -> gr_fmpz_mod_ctx:
+        """Create a new context for integers modulo `n`.
+
+        >>> from flint.types._gr import gr_fmpz_mod_ctx
+        >>> Z64 = gr_fmpz_mod_ctx.new(2**64 + 1)
+        >>> Z64
+        gr_fmpz_mod_ctx(18446744073709551617)
+        >>> Z64(2**64+1)
+        0
+        >>> Z64(2)**64 + 2
+        1
+        >>> Z64.is_prime
+        False
+        >>> Z64.modulus()
+        18446744073709551617
+        """
+        n = fmpz(n)
+        return gr_fmpz_mod_ctx._new(n)
+
+    @property
+    def is_prime(self) -> bool:
+        return self.is_field
+
+    def modulus(self):
+        cdef fmpz n = fmpz.__new__(fmpz)
+        fmpz_init_set(n.val, self.n)
+        return n
+
+    def __repr__(self):
+        return f"gr_fmpz_mod_ctx({self.modulus()})"
+
+
+@cython.no_gc
+cdef class gr_fq_ctx(gr_scalar_ctx):
+    r"""Context for finite fields, `\mathbb{F}_q` where `q = p^d`.
+
+    The characteristic `p` must be a prime number and the degree `d` must be
+    positive. The characteristic `p` can be arbitrary precision. Use
+    :class:`gr_fq_nmod_ctx` for word-sized characteristic. For very small
+    characteristic, use :class:`gr_fq_zech_ctx`.
+    """
+
+    @staticmethod
+    def new(p, d, name=None) -> gr_fq_ctx:
+        """Create a new context for finite fields.
+
+        >>> from flint.types._gr import gr_fq_ctx
+        >>> F9 = gr_fq_ctx.new(3, 2)
+        >>> F9
+        gr_fq_ctx(3, 2)
+        >>> F9(2) + F9(3)
+        2
+        >>> F9.characteristic()
+        3
+        >>> F9.degree()
+        2
+        >>> F9.gen()
+        a
+        >>> a = F9.gen()
+        >>> (1 + a) ** 2 + a
+        a+2
+        """
+        cdef bytes name_b
+        cdef char *name_c
+        if name is not None:
+            name_b = name.encode('utf-8')
+            name_c = name_b
+        else:
+            name_c = NULL
+        p = fmpz(p)
+        return gr_fq_ctx._new(p, d, name_c)
+
+    def characteristic(self):
+        cdef fmpz p = fmpz.__new__(fmpz)
+        fmpz_init_set(p.val, self.p)
+        return p
+
+    def degree(self):
+        return self.d
+
+    def __repr__(self):
+        return f"gr_fq_ctx({self.characteristic()}, {self.degree()})"
+
+
+@cython.no_gc
+cdef class gr_fq_nmod_ctx(gr_scalar_ctx):
+    r"""Context for finite fields, `\mathbb{F}_q` where `q = p^d`.
+
+    The characteristic `p` must be a prime number and the degree `d` must be
+    positive. The characteristic `p` must be word-sized. Use :class:`gr_fq_ctx`
+    for arbitrary precision characteristic. For very small characteristic, use
+    :class:`gr_fq_zech_ctx`.
+    """
+
+    @staticmethod
+    def new(p, d, name=None) -> gr_fq_nmod_ctx:
+        """Create a new context for finite fields.
+
+        >>> from flint.types._gr import gr_fq_nmod_ctx
+        >>> F9 = gr_fq_nmod_ctx.new(3, 2)
+        >>> F9
+        gr_fq_nmod_ctx(3, 2)
+        >>> F9(2) + F9(3)
+        2
+        >>> F9.characteristic()
+        3
+        >>> F9.degree()
+        2
+        >>> F9.gen()
+        a
+        >>> a = F9.gen()
+        >>> (1 + a) ** 2 + a
+        a+2
+        """
+        cdef bytes name_b
+        cdef char *name_c
+        if name is not None:
+            name_b = name.encode('utf-8')
+            name_c = name_b
+        else:
+            name_c = NULL
+        return gr_fq_nmod_ctx._new(p, d, name_c)
+
+    def characteristic(self):
+        return self.p
+
+    def degree(self):
+        return self.d
+
+    def __repr__(self):
+        return f"gr_fq_nmod_ctx({self.characteristic()}, {self.degree()})"
+
+
+@cython.no_gc
+cdef class gr_fq_zech_ctx(gr_scalar_ctx):
+    r"""Context for finite fields, `\mathbb{F}_q` where `q = p^d`.
+
+    The characteristic `p` must be a prime number and the degree `d` must be
+    positive. The characteristic `p` must be small. Use :class:`gr_nmod_ctx`
+    for word-sized characteristic and :class:`gr_fq_ctx` for arbitrary
+    precision characteristic.
+
+    .. warning::
+        This type is possibly not working correctly. Use with caution.
+    """
+
+    @staticmethod
+    def new(p, d, name=None) -> gr_fq_zech_ctx:
+        """Create a new context for finite fields with small characteristic.
+
+        >>> from flint.types._gr import gr_fq_zech_ctx
+        >>> F9 = gr_fq_zech_ctx.new(3, 2)
+        >>> F9
+        gr_fq_zech_ctx(3, 2)
+        >>> F9(2) + F9(3) # XXX: Is this correct?
+        a^4
+        >>> F9.characteristic()
+        3
+        >>> F9.degree()
+        2
+        >>> F9.gen()
+        a^1
+        >>> a = F9.gen()
+        >>> (1 + a) ** 2 + a  # doctest: +SKIP
+        a+2
+        """
+        cdef bytes name_b
+        cdef char *name_c
+        if name is not None:
+            name_b = name.encode('utf-8')
+            name_c = name_b
+        else:
+            name_c = NULL
+        return gr_fq_zech_ctx._new(p, d, name_c)
+
+    def characteristic(self):
+        return self.p
+
+    def degree(self):
+        return self.d
+
+    def __repr__(self):
+        return f"gr_fq_zech_ctx({self.characteristic()}, {self.degree()})"
+
+
+@cython.no_gc
+cdef class gr_nf_ctx(gr_scalar_ctx):
+    r"""Context for number fields, `\mathbb{Q}(\alpha)`.
+
+    The number field is defined by a minimal polynomial `f(x)` in
+    `\mathbb{Q}[x]` and represents the field `\mathbb{Q}(\alpha)` where
+    `\alpha` is a root of `f(x)`. The minimal polynomial `f(x)` must be
+    irreducible.
+    """
+
+    @staticmethod
+    def new(poly) -> gr_nf_ctx:
+        """Create a new context for number fields.
+
+        >>> from flint.types._gr import gr_nf_ctx
+        >>> Qa = gr_nf_ctx.new([-2, 0, 1])
+        >>> Qa
+        gr_nf_ctx(x^2 + (-2))
+        >>> Qa.modulus()
+        x^2 + (-2)
+        >>> a = Qa.gen()
+        >>> a
+        a
+        >>> a**2
+        2
+        >>> (1 + a) ** 2
+        2*a+3
+        >>> (1 + a) / 2
+        1/2*a+1/2
+        """
+        poly = fmpq_poly(poly)
+        return gr_nf_ctx._new(poly)
+
+    def modulus(self):
+        cdef fmpq_poly poly = fmpq_poly.__new__(fmpq_poly)
+        fmpq_poly_init(poly.val)
+        fmpq_poly_set(poly.val, self.poly)
+        return poly
+
+    def __repr__(self):
+        return f"gr_nf_ctx({self.modulus()})"
+
+
+@cython.no_gc
+cdef class gr_nf_fmpz_poly_ctx(gr_scalar_ctx):
+    r"""Context for number fields, `\mathbb{Q}(\alpha)`.
+
+    The number field is defined by a minimal polynomial `f(x)` in
+    `\mathbb{Z}[x]` and represents the ring `\mathbb{Z}[\alpha]` where
+    `\alpha` is a root of `f(x)`. The minimal polynomial `f(x)` must be
+    irreducible.
+    """
+
+    @staticmethod
+    def new(poly) -> gr_nf_fmpz_poly_ctx:
+        """Create a new context for number fields.
+
+        >>> from flint.types._gr import gr_nf_fmpz_poly_ctx
+        >>> Qa = gr_nf_fmpz_poly_ctx.new([-2, 0, 1])
+        >>> Qa
+        gr_nf_fmpz_poly_ctx(x^2 + (-2))
+        >>> Qa.modulus()
+        x^2 + (-2)
+        >>> a = Qa.gen()
+        >>> a
+        a
+        >>> a**2
+        2
+        >>> (1 + a) ** 2
+        2*a+3
+        """
+        poly = fmpz_poly(poly)
+        return gr_nf_fmpz_poly_ctx._new(poly)
+
+    def modulus(self):
+        cdef fmpz_poly poly = fmpz_poly.__new__(fmpz_poly)
+        fmpz_poly_init(poly.val)
+        fmpz_poly_set(poly.val, self.poly)
+        return poly
+
+    def __repr__(self):
+        return f"gr_nf_fmpz_poly_ctx({self.modulus()})"
+
+
+@cython.no_gc
+cdef class gr_real_qqbar_ctx(gr_scalar_ctx):
+    r"""Context for real algebraic numbers, `\mathbb{R} \cap \overline{\mathbb{Q}}`.
+
+    The real algebraic numbers are a subset of the real numbers that are roots
+    of non-zero polynomials with integer coefficients. See
+    :class:`gr_complex_qqbar_ctx` for the full algebraic closure of the
+    rationals. See :class:`gr_nf_ctx` for number fields with a single
+    generator.
+    """
+    @staticmethod
+    def new(deg_limit=-1, bits_limit=-1) -> gr_real_qqbar_ctx:
+        """Create a new context for real algebraic numbers.
+
+        >>> from flint.types._gr import gr_real_qqbar_ctx
+        >>> R = gr_real_qqbar_ctx.new()
+        >>> R
+        gr_real_qqbar_ctx(-1, -1)
+        >>> R(2).sqrt()
+        Root a = 1.41421 of a^2-2
+        """
+        assert deg_limit == -1 and bits_limit == -1 or deg_limit > 0 and bits_limit > 0
+        return gr_real_qqbar_ctx._new(deg_limit, bits_limit)
+
+    def limits(self):
+        return {'deg_limit': self.deg_limit, 'bits_limit': self.bits_limit}
+
+    def __repr__(self):
+        return f"gr_real_qqbar_ctx({self.deg_limit}, {self.bits_limit})"
+
+
+@cython.no_gc
+cdef class gr_complex_qqbar_ctx(gr_scalar_ctx):
+    r"""Context for algebraic numbers, `\overline{\mathbb{Q}}`.
+
+    The algebraic numbers are a field that contains the rationals and all roots
+    of non-zero polynomials with integer coefficients. See
+    :class:`gr_real_qqbar_ctx` for the real algebraic numbers,
+    :class:`gr_nf_ctx` for number fields with a single generator.
+    """
+    @staticmethod
+    def new(deg_limit=-1, bits_limit=-1) -> gr_complex_qqbar_ctx:
+        """Create a new context for algebraic numbers.
+
+        >>> from flint.types._gr import gr_complex_qqbar_ctx
+        >>> C = gr_complex_qqbar_ctx.new()
+        >>> C
+        gr_complex_qqbar_ctx(-1, -1)
+        >>> C(-2).sqrt()
+        Root a = 1.41421*I of a^2+2
+        """
+        assert deg_limit == -1 and bits_limit == -1 or deg_limit > 0 and bits_limit > 0
+        return gr_complex_qqbar_ctx._new(deg_limit, bits_limit)
+
+    def limits(self):
+        return {'deg_limit': self.deg_limit, 'bits_limit': self.bits_limit}
+
+    def __repr__(self):
+        return f"gr_complex_qqbar_ctx({self.deg_limit}, {self.bits_limit})"
+
+
+@cython.no_gc
+cdef class gr_real_ca_ctx(gr_scalar_ctx):
+    r"""Context for calcium exact real numbers, `\mathbb{R}`."""
+
+    @staticmethod
+    def new(**options) -> gr_real_ca_ctx:
+        """Create a new context for calcium exact real numbers.
+
+        >>> from flint.types._gr import gr_real_ca_ctx
+        >>> R = gr_real_ca_ctx.new()
+        >>> R
+        gr_real_ca_ctx({})
+        >>> R(2).sqrt()
+        1.41421 {a where a = 1.41421 [a^2-2=0]}
+        """
+        return gr_real_ca_ctx._new(options)
+
+    def __repr__(self):
+        return f"gr_real_ca_ctx({self.options})"
+
+
+@cython.no_gc
+cdef class gr_complex_ca_ctx(gr_scalar_ctx):
+    r"""Context for calcium exact complex numbers, `\mathbb{C}`."""
+
+    @staticmethod
+    def new(**options) -> gr_complex_ca_ctx:
+        """Create a new context for calcium exact complex numbers.
+
+        >>> from flint.types._gr import gr_complex_ca_ctx
+        >>> C = gr_complex_ca_ctx.new()
+        >>> C
+        gr_complex_ca_ctx({})
+        >>> C(2).sqrt()
+        1.41421 {a where a = 1.41421 [a^2-2=0]}
+        """
+        return gr_complex_ca_ctx._new(options)
+
+    def __repr__(self):
+        return f"gr_complex_ca_ctx({self.options})"
+
+
+@cython.no_gc
+cdef class gr_real_algebraic_ca_ctx(gr_scalar_ctx):
+    r"""Context for calcium exact real algebraic numbers, `\mathbb{R} \cap \overline{\mathbb{Q}}`."""
+
+    @staticmethod
+    def new(**options) -> gr_real_algebraic_ca_ctx:
+        """Create a new context for calcium exact real algebraic numbers.
+
+        >>> from flint.types._gr import gr_real_algebraic_ca_ctx
+        >>> R = gr_real_algebraic_ca_ctx.new()
+        >>> R
+        gr_real_algebraic_ca_ctx({})
+        >>> R(2).sqrt()
+        1.41421 {a where a = 1.41421 [a^2-2=0]}
+        """
+        return gr_real_algebraic_ca_ctx._new(options)
+
+    def __repr__(self):
+        return f"gr_real_algebraic_ca_ctx({self.options})"
+
+
+@cython.no_gc
+cdef class gr_complex_algebraic_ca_ctx(gr_scalar_ctx):
+    r"""Context for calcium exact complex algebraic numbers, `\overline{\mathbb{Q}}`."""
+
+    @staticmethod
+    def new(**options) -> gr_complex_algebraic_ca_ctx:
+        """Create a new context for calcium exact complex algebraic numbers.
+
+        >>> from flint.types._gr import gr_complex_algebraic_ca_ctx
+        >>> C = gr_complex_algebraic_ca_ctx.new()
+        >>> C
+        gr_complex_algebraic_ca_ctx({})
+        >>> C(2).sqrt()
+        1.41421 {a where a = 1.41421 [a^2-2=0]}
+        """
+        return gr_complex_algebraic_ca_ctx._new(options)
+
+    def __repr__(self):
+        return f"gr_complex_algebraic_ca_ctx({self.options})"
+
+
+@cython.no_gc
+cdef class gr_complex_extended_ca_ctx(gr_scalar_ctx):
+    r"""Context for calcium exact extended complex numbers, `\mathbb{C} \cup \{\infty\}`."""
+
+    @staticmethod
+    def new(**options) -> gr_complex_extended_ca_ctx:
+        """Create a new context for calcium exact extended complex numbers.
+
+        >>> from flint.types._gr import gr_complex_extended_ca_ctx
+        >>> C = gr_complex_extended_ca_ctx.new()
+        >>> C
+        gr_complex_extended_ca_ctx({})
+        >>> C(1) / 0
+        UnsignedInfinity
+        """
+        return gr_complex_extended_ca_ctx._new(options)
+
+    def __repr__(self):
+        return f"gr_complex_extended_ca_ctx({self.options})"
+
+
+@cython.no_gc
+cdef class gr_real_float_arf_ctx(gr_scalar_ctx):
+    r"""Context for arbitrary precision approximate real numbers, `\mathbb{R}`."""
+
+    @staticmethod
+    def new(prec) -> gr_real_float_arf_ctx:
+        """Create a new context for arbitrary precision approximate real numbers.
+
+        >>> from flint.types._gr import gr_real_float_arf_ctx
+        >>> R = gr_real_float_arf_ctx.new(10)
+        >>> R
+        gr_real_float_arf_ctx(10)
+        >>> R(2).sqrt()
+        1.414
+        >>> R.real_prec = 20
+        >>> R(2).sqrt()
+        1.414213
+        """
+        return gr_real_float_arf_ctx._new(prec)
+
+    def __repr__(self):
+        return f"gr_real_float_arf_ctx({self.prec})"
+
+
+@cython.no_gc
+cdef class gr_complex_float_acf_ctx(gr_scalar_ctx):
+    r"""Context for arbitrary precision approximate complex numbers, `\mathbb{C}`."""
+
+    @staticmethod
+    def new(prec) -> gr_complex_float_acf_ctx:
+        """Create a new context for arbitrary precision approximate complex numbers.
+
+        >>> from flint.types._gr import gr_complex_float_acf_ctx
+        >>> C = gr_complex_float_acf_ctx.new(10)
+        >>> C
+        gr_complex_float_acf_ctx(10)
+        >>> C(-2).sqrt()
+        1.414*I
+        >>> C.real_prec = 20
+        >>> C(-2).sqrt()
+        1.414213*I
+        """
+        return gr_complex_float_acf_ctx._new(prec)
+
+    def __repr__(self):
+        return f"gr_complex_float_acf_ctx({self.prec})"
+
+
+@cython.no_gc
+cdef class gr_real_arb_ctx(gr_scalar_ctx):
+    r"""Context for Arb arbitrary precision real ball arithmetic, `\mathbb{R}`."""
+
+    @staticmethod
+    def new(prec) -> gr_real_arb_ctx:
+        """Create a new context for arbitrary precision real ball arithmetic.
+
+        >>> from flint.types._gr import gr_real_arb_ctx
+        >>> R = gr_real_arb_ctx.new(10)
+        >>> R
+        gr_real_arb_ctx(10)
+        >>> R(2).sqrt()
+        [1.41 +/- 6.02e-3]
+        >>> R.real_prec = 20
+        >>> R(2).sqrt()
+        [1.41421 +/- 5.09e-6]
+        """
+        return gr_real_arb_ctx._new(prec)
+
+    def __repr__(self):
+        return f"gr_real_arb_ctx({self.prec})"
+
+
+@cython.no_gc
+cdef class gr_complex_acb_ctx(gr_scalar_ctx):
+    r"""Context for Arb arbitrary precision complex ball arithmetic, `\mathbb{C}`."""
+
+    @staticmethod
+    def new(prec) -> gr_complex_acb_ctx:
+        """Create a new context for arbitrary precision complex ball arithmetic.
+
+        >>> from flint.types._gr import gr_complex_acb_ctx
+        >>> C = gr_complex_acb_ctx.new(10)
+        >>> C
+        gr_complex_acb_ctx(10)
+        >>> C(-2).sqrt()
+        [1.41 +/- 6.02e-3]*I
+        >>> C.real_prec = 20
+        >>> C(-2).sqrt()
+        [1.41421 +/- 5.09e-6]*I
+        """
+        return gr_complex_acb_ctx._new(prec)
+
+    def __repr__(self):
+        return f"gr_complex_acb_ctx({self.prec})"
+
+
+@cython.no_gc
+cdef class gr_gr_poly_ctx(gr_poly_ctx):
+    r"""Context for dense univariate polynomial rings, `\mathbb{D}[x]`."""
+
+    @staticmethod
+    def new(base_ring) -> gr_gr_poly_ctx:
+        """Create a new context for dense univariate polynomial rings.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx, gr_gr_poly_ctx
+        >>> ZI = gr_fmpzi_ctx
+        >>> R = gr_gr_poly_ctx.new(ZI)
+        >>> R
+        gr_gr_poly_ctx(gr_fmpzi_ctx)
+        >>> R.base_ring
+        gr_fmpzi_ctx
+        >>> R.gen()
+        x
+        >>> ZI.gen()
+        I
+        >>> R.gens_recursive()
+        [I, x]
+        >>> I, x = R.gens_recursive()
+        >>> (I + x) ** 2
+        -1 + (2*I)*x + x^2
+        """
+        return gr_gr_poly_ctx._new(base_ring)
+
+    def __repr__(self):
+        return f"gr_gr_poly_ctx({self.base_ring})"
+
+    @property
+    def base_ring(self):
+        return self.base_ctx
+
+
+@cython.no_gc
+cdef class gr_gr_mpoly_ctx(gr_mpoly_ctx):
+    r"""Context for dense multivariate polynomial rings, `\mathbb{D}[x, y, ...]`."""
+
+    @staticmethod
+    def new(base_ring, names, order=None) -> gr_gr_mpoly_ctx:
+        """Create a new context for dense multivariate polynomial rings.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx, gr_gr_mpoly_ctx
+        >>> ZI = gr_fmpzi_ctx
+        >>> R = gr_gr_mpoly_ctx.new(ZI, ['x', 'y'])
+        >>> R
+        gr_gr_mpoly_ctx(gr_fmpzi_ctx, ('x', 'y'), Ordering.lex)
+        >>> R.base_ring
+        gr_fmpzi_ctx
+        >>> R.names
+        ('x', 'y')
+        >>> R.nvars
+        2
+        >>> R.order
+        <Ordering.lex: 'lex'>
+        >>> R.gens()
+        [x, y]
+        >>> R.gens_recursive()
+        [I, x, y]
+        >>> I, x, y = R.gens_recursive()
+        >>> (I + x + y) ** 2
+        x^2 + 2*x*y + (2*I)*x + y^2 + (2*I)*y - 1
+        """
+        if order is None:
+            order = Ordering.lex
+        return gr_gr_mpoly_ctx._new(base_ring, tuple(names), order)
+
+    def __repr__(self):
+        return f"gr_gr_mpoly_ctx({self.base_ring}, {self.names}, {self.order})"
+
+    @property
+    def base_ring(self):
+        return self.base_ctx
+
+    @property
+    def names(self):
+        return self._names
+
+    @property
+    def nvars(self):
+        return self._nvars
+
+    @property
+    def order(self):
+        return ordering_c_to_py(self._order)
+
+
+# @cython.no_gc
+# cdef class gr_mpoly_q_ctx(gr_mpoly_ctx):
+#
+#     @staticmethod
+#     def new(nvars, order=None) -> gr_gr_mpoly_q_ctx:
+#         if order is None:
+#             order = Ordering.lex
+#         return gr_mpoly_q_ctx._new(nvars, order)
+#
+#     def __repr__(self):
+#         return f"gr_mpoly_q_ctx({self.nvars}, {self.order})"
+#
+#     @property
+#     def nvars(self):
+#         return self._nvars
+#
+#     @property
+#     def order(self):
+#         return ordering_c_to_py(self._order)
+
+
+# @cython.no_gc
+# cdef class gr_series_mod_gr_poly_ctx(gr_ctx):
+#
+#     @staticmethod
+#     def new(base_ring, n) -> gr_series_mod_gr_poly_ctx:
+#         return gr_series_mod_gr_poly_ctx._new(base_ring, n)
+#
+#     def __repr__(self):
+#         return f"gr_series_mod_gr_poly_ctx({self.base_ring}, {self.nvars})"
+#
+#     @property
+#     def base_ring(self):
+#         return self.base_ctx
+#
+#     @property
+#     def n(self):
+#         return self._n
+
+
+@cython.no_gc
+cdef class gr_series_ctx(gr_ctx):
+    r"""Context for series with precision `n`, `\mathbb{D}[[x]]`."""
+
+    @staticmethod
+    def new(base_ring, prec) -> gr_series_ctx:
+        """Create a new context for series with precision `n`.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx, gr_series_ctx
+        >>> ZI = gr_fmpzi_ctx
+        >>> R = gr_series_ctx.new(ZI, 10)
+        >>> R
+        gr_series_ctx(gr_fmpzi_ctx, 10)
+        >>> R.base_ring
+        gr_fmpzi_ctx
+        >>> R.prec
+        10
+        >>> R.gens_recursive()
+        [I, x]
+        >>> I, x = R.gens_recursive()
+        >>> 1 / (I - x)
+        -I - x + I*x^2 + x^3 - I*x^4 - x^5 + I*x^6 + x^7 - I*x^8 - x^9 + O(x^10)
+        """
+        return gr_series_ctx._new(base_ring, prec)
+
+    def __repr__(self):
+        return f"gr_series_ctx({self.base_ring}, {self.prec})"
+
+    @property
+    def base_ring(self):
+        return self.base_ctx
+
+    @property
+    def prec(self):
+        return self._prec
+
+
+@cython.no_gc
+cdef class gr(flint_scalar):
+    """Type of elements of :class:`gr_ctx` domains."""
+
+    def __dealloc__(self):
+        if self._init:
+            self._init = False
+            gr_heap_clear(self.pval, self.ctx.ctx_t)
+
+    def __repr__(self):
+        return self.ctx.to_str(self)
+
+    def is_zero(self):
+        """Return whether the element is zero (may return ``None``).
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(0).is_zero()
+        True
+        """
+        return truth_to_py(self._is_zero())
+
+    def is_one(self):
+        """Return whether the element is one (may return ``None``).
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(1).is_one()
+        True
+        """
+        return truth_to_py(self._is_one())
+
+    def is_neg_one(self):
+        """Return whether the element is negative one (may return ``None``).
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(-1).is_neg_one()
+        True
+        """
+        return truth_to_py(self._is_neg_one())
+
+#     def is_integer(self):
+#         """Return whether the element is an integer (may return ``None``).
+#
+#         >>> from flint.types._gr import gr_fmpq_ctx as Q
+#         >>> Q(2).is_integer()
+#         True
+#         >>> Q(2, 3).is_integer()
+#         False
+#         """
+#         return truth_to_py(self._is_integer())
+#
+#     def is_rational(self):
+#         """Return whether the element is a rational number (may return ``None``).
+#
+#         >>> from flint.types._gr import gr_nf_ctx
+#         >>> Qa = gr_nf_ctx.new([-2, 0, 1])
+#         >>> a = Qa.gen()
+#         >>> a.is_rational()
+#         False
+#         >>> (a**2).is_rational()
+#         True
+#         """
+#         return truth_to_py(self._is_rational())
+
+    def __bool__(self):
+        return not truth_to_bool(self._is_zero())
+
+    def __richcmp__(self, other, int op):
+        cdef int err
+        cdef gr other_gr
+        cdef truth_t res
+
+        if isinstance(other, int):
+            other_gr = self.ctx.new_gr()
+            err = gr_set_si(other_gr.pval, other, self.ctx.ctx_t)
+            if err != GR_SUCCESS:
+                raise self.ctx._error(err, "Cannot set gr from int")
+        elif isinstance(other, gr):
+            other_gr = other
+        else:
+            raise NotImplementedError("Cannot compare gr with non-gr objects")
+
+        if self.ctx != other_gr.ctx:
+            raise NotImplementedError("Cannot compare gr with different contexts")
+
+        if op == 0:
+            raise NotImplementedError("Cannot compare gr with <")
+            # res = self._lt(other_gr)
+        elif op == 1:
+            raise NotImplementedError("Cannot compare gr with <=")
+            # res = self._le(other_gr)
+        elif op == 2:
+            res = self._equal(other_gr)
+        elif op == 3:
+            res = self._equal(other_gr)
+            if res == T_TRUE:
+                res = T_FALSE
+            elif res == T_FALSE:
+                res = T_TRUE
+        elif op == 4:
+            raise NotImplementedError("Cannot compare gr with >")
+            # res = self._gt(other_gr)
+        elif op == 5:
+            raise NotImplementedError("Cannot compare gr with >=")
+            # res = self._ge(other_gr)
+        else:
+            assert False, "Invalid rich comparison operator"
+
+        return truth_to_py(res)
+
+    def __neg__(self) -> gr:
+        return self._neg()
+
+    # XXX: Maybe +a should return a copy or for arb it should round...
+    # def __pos__(self) -> gr:
+    #     return self
+
+    def __add__(self, other) -> gr:
+        cdef gr other_gr
+        if isinstance(other, gr):
+            other_gr = other
+            if self.ctx != other_gr.ctx:
+                raise ValueError("Cannot add elements from different contexts")
+            return self._add(other_gr)
+        elif isinstance(other, int):
+            return self._add_si(other)
+        else:
+            return NotImplemented
+
+    def __radd__(self, other) -> gr:
+        if isinstance(other, int):
+            return self._add_si(other)
+        else:
+            return NotImplemented
+
+    def __sub__(self, other) -> gr:
+        cdef gr other_gr
+        if isinstance(other, gr):
+            other_gr = other
+            if self.ctx != other_gr.ctx:
+                raise ValueError("Cannot subtract elements from different contexts")
+            return self._sub(other_gr)
+        elif isinstance(other, int):
+            return self._sub_si(other)
+        else:
+            return NotImplemented
+
+    def __rsub__(self, other) -> gr:
+        if isinstance(other, int):
+            return self._neg()._add_si(other)
+        else:
+            return NotImplemented
+
+    def __mul__(self, other) -> gr:
+        cdef gr other_gr
+        if isinstance(other, gr):
+            other_gr = other
+            if self.ctx != other_gr.ctx:
+                raise ValueError("Cannot multiply elements from different contexts")
+            return self._mul(other_gr)
+        elif isinstance(other, int):
+            return self._mul_si(other)
+        else:
+            return NotImplemented
+
+    def __rmul__(self, other) -> gr:
+        if isinstance(other, int):
+            return self._mul_si(other)
+        else:
+            return NotImplemented
+
+    def __truediv__(self, other) -> gr:
+        cdef gr other_gr
+        if isinstance(other, gr):
+            other_gr = other
+            if self.ctx != other_gr.ctx:
+                raise ValueError("Cannot divide elements from different contexts")
+            return self._div(other_gr)
+        elif isinstance(other, int):
+            return self._div_si(other)
+        else:
+            return NotImplemented
+
+    def __rtruediv__(self, other) -> gr:
+        if isinstance(other, int):
+            return self._div_si(other)._inv()
+        else:
+            return NotImplemented
+
+    def __pow__(self, other) -> gr:
+        if isinstance(other, int):
+            return self._pow_si(other)
+        else:
+            return NotImplemented
+
+    def is_square(self):
+        """Return whether the element is a square (may return ``None``).
+
+        >>> from flint.types._gr import gr_fmpq_ctx as Q
+        >>> Q(2).is_square()
+        False
+        >>> Q(4).is_square()
+        True
+        >>> Q(4).sqrt()
+        2
+        """
+        return truth_to_py(self._is_square())
+
+    def sqrt(self):
+        """Return the square root of the element if it exists.
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(4).sqrt()
+        2
+        """
+        return self._sqrt()
+
+    def rsqrt(self):
+        """Return the reciprocal square root of the element if it exists.
+
+        >>> from flint.types._gr import gr_fmpq_ctx as Q
+        >>> Q(4).rsqrt()
+        1/2
+        """
+        return self._rsqrt()
+
+    def gcd(self, other):
+        """Return the greatest common divisor of two elements.
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(4).gcd(Z(6))
+        2
+        """
+        cdef gr other_gr
+        if not isinstance(other, gr):
+            raise TypeError("gcd when other is not gr.")
+        other_gr = other
+        if not self.ctx == other_gr.ctx:
+            raise TypeError("gcd of gr with different contexts.")
+        return self._gcd(other_gr)
+
+    def lcm(self, other):
+        """Return the least common multiple of two elements.
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(4).lcm(Z(6))
+        12
+        """
+        cdef gr other_gr
+        if not isinstance(other, gr):
+            raise TypeError("gcd when other is not gr.")
+        other_gr = other
+        if not self.ctx == other_gr.ctx:
+            raise TypeError("gcd of gr with different contexts.")
+        return self._lcm(other_gr)
+
+    def factor(self):
+        """Return the factorization of the element.
+
+        >>> from flint.types._gr import gr_fmpz_ctx as Z
+        >>> Z(12).factor()
+        (1, [(2, 2), (3, 1)])
+        """
+        return self._factor()
+
+    def numer(self) -> gr:
+        """Return the numerator of the element.
+
+        >>> from flint.types._gr import gr_fmpq_ctx as Q
+        >>> q = Q(2) / 3
+        >>> q
+        2/3
+        >>> q.numer()
+        2
+
+        See also :meth:`denom`.
+        """
+        return self._numerator()
+
+    def denom(self) -> gr:
+        """Return the denominator of the element.
+
+        >>> from flint.types._gr import gr_fmpq_ctx as Q
+        >>> q = Q(2) / 3
+        >>> q
+        2/3
+        >>> q.denom()
+        3
+
+        See also :meth:`numer`.
+        """
+        return self._denominator()
+
+    def __floor__(self) -> gr:
+        return self._floor()
+
+    def __ceil__(self) -> gr:
+        return self._ceil()
+
+    def __trunc__(self) -> gr:
+        return self._trunc()
+
+    def __round__(self, ndigits: int = 0) -> gr:
+        if ndigits != 0:
+            raise NotImplementedError("Rounding to a specific number of digits is not supported")
+        return self._nint()
+
+    # def __int__(self) -> int:
+    #     return self._floor().to_int()
+
+    # def __float__(self) -> float:
+    #     return ...
+
+    def __abs__(self) -> gr:
+        return self._abs()
+
+    def conjugate(self) -> gr:
+        """Return complex conjugate of the element.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx as ZI
+        >>> I = ZI.gen()
+        >>> (1 + I).conjugate()
+        (1-I)
+        """
+        return self._conj()
+
+    @property
+    def real(self) -> gr:
+        """Return the real part of the element.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx as ZI
+        >>> I = ZI.gen()
+        >>> (1 + I).real
+        1
+        """
+        return self._re()
+
+    @property
+    def imag(self) -> gr:
+        """Return the imaginary part of the element.
+
+        >>> from flint.types._gr import gr_fmpzi_ctx as ZI
+        >>> I = ZI.gen()
+        >>> (1 + I).imag
+        1
+        """
+        return self._im()
+
+    # XXX: Return -1, 0, 1 as int?
+    def sgn(self) -> gr:
+        """Return the sign of the element.
+
+        >>> from flint.types._gr import gr_fmpq_ctx as Q
+        >>> Q(2).sgn()
+        1
+        >>> Q(-2).sgn()
+        -1
+        >>> Q(0).sgn()
+        0
+        """
+        return self._sgn()
+
+    def csgn(self) -> gr:
+        """Return the complex sign of the element.
+
+        >>> from flint.types._gr import gr_complex_acb_ctx
+        >>> C = gr_complex_acb_ctx.new(10)
+        >>> (1 + C.i()).csgn()
+        1
+        """
+        return self._csgn()
+
+    def arg(self) -> gr:
+        """Return the argument of the element.
+
+        >>> from flint.types._gr import gr_complex_acb_ctx
+        >>> C = gr_complex_acb_ctx.new(10)
+        >>> (1 + C.i()).arg()
+        [0.785 +/- 6.45e-4]
+        """
+        return self._arg()
diff --git a/src/flint/types/meson.build b/src/flint/types/meson.build
index 6f32aff5..b5c2f953 100644
--- a/src/flint/types/meson.build
+++ b/src/flint/types/meson.build
@@ -47,6 +47,8 @@ exts = [
   'fmpz_mpoly',
   'fmpz_mod_mpoly',
   'fmpq_mpoly',
+
+  '_gr',
 ]
 
 if have_acb_theta
diff --git a/src/flint/utils/flint_exceptions.py b/src/flint/utils/flint_exceptions.py
index 9d9c271b..05c6e5e5 100644
--- a/src/flint/utils/flint_exceptions.py
+++ b/src/flint/utils/flint_exceptions.py
@@ -16,3 +16,18 @@ class IncompatibleContextError(FlintError):
     incompatible contexts.
     """
     pass
+
+
+class UnableError(FlintError):
+    """
+    Exception intended to be called when the implementation is unable to
+    perform the requested operation.
+    """
+    pass
+
+
+class UnknownError(FlintError):
+    """
+    Exception intended to be called when the value of a predicate is unknown.
+    """
+    pass