diff --git a/doc/source/ca.rst b/doc/source/ca.rst
index 891b511b4a..d6cad232e7 100644
--- a/doc/source/ca.rst
+++ b/doc/source/ca.rst
@@ -1231,6 +1231,37 @@ Special functions
 
     Sets *res* to the imaginary error function of *x*.
 
+.. function:: void ca_elliptic_k(ca_t res, const ca_t m, ca_ctx_t ctx)
+
+    Sets *res* to the complete elliptic integral of the first kind `K(m)` of *m*, defined by
+    
+    .. math::
+
+        K(m) = \int_0^{\pi/2} \frac{dt}{\sqrt{1-m \sin^2 t}}
+                  = \int_0^1
+            \frac{dt}{\left(\sqrt{1-t^2}\right)\left(\sqrt{1-mt^2}\right)}.
+
+.. function:: void ca_elliptic_e(ca_t res, const ca_t m, ca_ctx_t ctx)
+
+    Sets *res* to the complete elliptic integral of the second kind `E(m)` of *m*, defined by
+    
+    .. math::
+
+        E(m) = \int_0^{\pi/2} \sqrt{1-m \sin^2 t} \, dt =
+                    \int_0^1
+                    \frac{\sqrt{1-mt^2}}{\sqrt{1-t^2}} \, dt.
+
+.. function:: void ca_elliptic_pi(ca_t res, const ca_t n, const ca_t m, ca_ctx_t ctx)
+
+    Sets *res* to the complete elliptic integral of the third kind `\Pi(n,m)` of *n* and *m*, defined by
+    
+    .. math::
+
+        \Pi(n, m) = \int_0^{\pi/2}
+            \frac{dt}{(1-n \sin^2 t) \sqrt{1-m \sin^2 t}} =
+            \int_0^1
+            \frac{dt}{(1-nt^2) \sqrt{1-t^2} \sqrt{1-mt^2}}.
+
 
 Numerical evaluation
 -------------------------------------------------------------------------------
diff --git a/src/ca.h b/src/ca.h
index ee12508ee6..cdc6244c7b 100644
--- a/src/ca.h
+++ b/src/ca.h
@@ -480,6 +480,10 @@ void ca_erfi(ca_t res, const ca_t x, ca_ctx_t ctx);
 
 void ca_gamma(ca_t res, const ca_t x, ca_ctx_t ctx);
 
+void ca_elliptic_k(ca_t res, const ca_t m, ca_ctx_t ctx);
+void ca_elliptic_e(ca_t res, const ca_t m, ca_ctx_t ctx);
+void ca_elliptic_pi(ca_t res, const ca_t n, const ca_t m, ca_ctx_t ctx);
+
 /* Numerical evaluation */
 
 void ca_get_acb_raw(acb_t res, const ca_t x, slong prec, ca_ctx_t ctx);
diff --git a/src/ca/elliptic.c b/src/ca/elliptic.c
new file mode 100644
index 0000000000..1dee10ae6b
--- /dev/null
+++ b/src/ca/elliptic.c
@@ -0,0 +1,156 @@
+/*
+    Copyright (C) 2020 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "ca.h"
+
+void
+ca_elliptic_k(ca_t res, const ca_t m, ca_ctx_t ctx)
+{
+    if (CA_IS_SPECIAL(m))
+    {
+        if (ca_check_is_pos_inf(m, ctx) == T_TRUE || ca_check_is_neg_inf(m, ctx) == T_TRUE ||
+            ca_check_is_pos_i_inf(m, ctx) == T_TRUE || ca_check_is_neg_i_inf(m, ctx) == T_TRUE)
+            ca_zero(res, ctx);
+        else if (ca_check_is_undefined(m, ctx) == T_TRUE || ca_check_is_uinf(m, ctx) == T_TRUE)
+            ca_undefined(res, ctx);
+        else
+            ca_unknown(res, ctx);
+        return;
+    }
+
+    if (ca_check_is_zero(m, ctx) == T_TRUE)
+    {
+        ca_pi(res, ctx);
+        ca_div_ui(res, res, 2, ctx);
+        return;
+    }
+    else if (ca_check_is_one(m, ctx) == T_TRUE)
+    {
+        ca_uinf(res, ctx);
+        return;
+    }
+
+    _ca_make_field_element(res, _ca_ctx_get_field_fx(ctx, CA_EllipticK, m), ctx);
+    fmpz_mpoly_q_gen(CA_MPOLY_Q(res), 0, CA_MCTX_1(ctx));
+}
+
+void
+ca_elliptic_e(ca_t res, const ca_t m, ca_ctx_t ctx)
+{
+    if (CA_IS_SPECIAL(m))
+    {
+        if (ca_check_is_pos_inf(m, ctx) == T_TRUE)
+            ca_pos_i_inf(res, ctx);
+        else if (ca_check_is_neg_inf(m, ctx) == T_TRUE)
+            ca_pos_inf(res, ctx);
+        else if (ca_check_is_pos_i_inf(m, ctx) == T_TRUE)
+        {
+            ca_pos_inf(res, ctx);
+            ca_t phase;
+            ca_init(phase, ctx);
+            ca_i(phase, ctx);
+            ca_sqrt(phase, phase, ctx);
+            ca_pow_ui(phase, phase, 3, ctx);
+            ca_neg(phase, phase, ctx);
+            ca_mul(res, res, phase, ctx);
+            ca_clear(phase, ctx);
+        }
+        else if (ca_check_is_neg_i_inf(m, ctx) == T_TRUE)
+        {
+            ca_pos_inf(res, ctx);
+            ca_t phase;
+            ca_init(phase, ctx);
+            ca_i(phase, ctx);
+            ca_sqrt(phase, phase, ctx);
+            ca_mul(res, res, phase, ctx);
+            ca_clear(phase, ctx);
+        }
+        else if (ca_check_is_undefined(m, ctx) == T_TRUE || ca_check_is_uinf(m, ctx) == T_TRUE)
+            ca_undefined(res, ctx);
+        else
+            ca_unknown(res, ctx);
+        return;
+    }
+
+    if (ca_check_is_zero(m, ctx) == T_TRUE)
+    {
+        ca_pi(res, ctx);
+        ca_div_ui(res, res, 2, ctx);
+        return;
+    }
+    else if (ca_check_is_one(m, ctx) == T_TRUE)
+    {
+        ca_one(res, ctx);
+        return;
+    }
+
+    _ca_make_field_element(res, _ca_ctx_get_field_fx(ctx, CA_EllipticE, m), ctx);
+    fmpz_mpoly_q_gen(CA_MPOLY_Q(res), 0, CA_MCTX_1(ctx));
+}
+
+void
+ca_elliptic_pi(ca_t res, const ca_t n, const ca_t m, ca_ctx_t ctx)
+{
+    if (CA_IS_SPECIAL(n))
+    {
+        if (ca_check_is_pos_inf(n, ctx) == T_TRUE || ca_check_is_neg_inf(n, ctx) == T_TRUE ||
+            ca_check_is_pos_i_inf(n, ctx) == T_TRUE || ca_check_is_neg_i_inf(n, ctx) == T_TRUE)
+            ca_zero(res, ctx);
+        else if (ca_check_is_undefined(n, ctx) == T_TRUE || ca_check_is_uinf(n, ctx) == T_TRUE)
+            ca_undefined(res, ctx);
+        else
+            ca_unknown(res, ctx);
+        return;
+    }
+    
+    if (CA_IS_SPECIAL(m))
+    {
+        if (ca_check_is_pos_inf(m, ctx) == T_TRUE || ca_check_is_neg_inf(m, ctx) == T_TRUE ||
+            ca_check_is_pos_i_inf(m, ctx) == T_TRUE || ca_check_is_neg_i_inf(m, ctx) == T_TRUE)
+            ca_zero(res, ctx);
+        else if (ca_check_is_undefined(m, ctx) == T_TRUE || ca_check_is_uinf(m, ctx) == T_TRUE)
+            ca_undefined(res, ctx);
+        else
+            ca_unknown(res, ctx);
+        return;
+    }
+    
+    if (ca_check_is_one(n, ctx) == T_TRUE || ca_check_is_one(m, ctx) == T_TRUE)
+    {
+        ca_uinf(res, ctx);
+        return;
+    }
+
+    if (ca_check_is_zero(n, ctx) == T_TRUE)
+    {
+        _ca_make_field_element(res, _ca_ctx_get_field_fx(ctx, CA_EllipticK, m), ctx);
+        fmpz_mpoly_q_gen(CA_MPOLY_Q(res), 0, CA_MCTX_1(ctx));
+        return;
+    }
+    else if (ca_check_is_zero(m, ctx) == T_TRUE)
+    {
+        ca_t nc;
+        ca_init(nc, ctx);
+        ca_one(nc, ctx);
+        ca_sub(nc, nc, n, ctx);
+        ca_sqrt(nc, nc, ctx);
+        ca_inv(nc, nc, ctx);
+        ca_div_ui(nc, nc, 2, ctx);
+        ca_pi(res, ctx);
+        ca_mul(res, res, nc, ctx);
+        ca_clear(nc, ctx);
+        return;
+    }
+
+    _ca_make_field_element(res, _ca_ctx_get_field_fxy(ctx, CA_EllipticPi, n, m), ctx);
+    fmpz_mpoly_q_gen(CA_MPOLY_Q(res), 0, CA_MCTX_1(ctx));
+}
+
diff --git a/src/ca/test/main.c b/src/ca/test/main.c
index 13dfc46875..c9512b2d55 100644
--- a/src/ca/test/main.c
+++ b/src/ca/test/main.c
@@ -19,6 +19,7 @@
 #include "t-ctx_init_clear.c"
 #include "t-div.c"
 #include "t-erf.c"
+#include "t-elliptic.c"
 #include "t-exp.c"
 #include "t-field_init_clear.c"
 #include "t-fmpz_mpoly_evaluate.c"
@@ -53,6 +54,7 @@ test_struct tests[] =
     TEST_FUNCTION(ca_ctx_init_clear),
     TEST_FUNCTION(ca_div),
     TEST_FUNCTION(ca_erf),
+    TEST_FUNCTION(ca_elliptic),
     TEST_FUNCTION(ca_exp),
     TEST_FUNCTION(ca_field_init_clear),
     TEST_FUNCTION(ca_fmpz_mpoly_evaluate),
diff --git a/src/ca/test/t-elliptic.c b/src/ca/test/t-elliptic.c
new file mode 100644
index 0000000000..bf239a3f0b
--- /dev/null
+++ b/src/ca/test/t-elliptic.c
@@ -0,0 +1,153 @@
+/*
+    Copyright (C) 2020 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "test_helpers.h"
+#include "acb_elliptic.h"
+#include "ca.h"
+
+TEST_FUNCTION_START(ca_elliptic, state)
+{
+    slong iter;
+
+    for (iter = 0; iter < 1000 * 0.1 * flint_test_multiplier(); iter++)
+    {
+        ca_ctx_t ctx;
+        ca_t x, y, ekx, eex, epxy;
+        acb_t ay, ekax, eeax, epaxy, aekx, aeex, aepxy;
+
+        ca_ctx_init(ctx);
+        ca_init(x, ctx);
+        ca_init(y, ctx);
+        ca_init(ekx, ctx);
+        ca_init(eex, ctx);
+        ca_init(epxy, ctx);
+        acb_init(ekax);
+        acb_init(eeax);
+        acb_init(epaxy);
+        acb_init(aekx);
+        acb_init(aeex);
+        acb_init(aepxy);
+        acb_init(ay);
+
+        ca_randtest(x, state, 3, 5, ctx);
+        if (n_randint(state, 2))
+            ca_randtest(y, state, 3, 5, ctx);
+        else
+            ca_randtest_rational(y, state, 5, ctx);
+        if (n_randint(state, 2))
+            ca_add(y, y, x, ctx);
+
+        ca_elliptic_k(ekx, x, ctx);
+        ca_elliptic_e(eex, x, ctx);
+        ca_elliptic_pi(epxy, x, y, ctx);
+
+        ca_get_acb(aekx, ekx, 53, ctx);
+        ca_get_acb(aeex, eex, 53, ctx);
+        ca_get_acb(aepxy, epxy, 53, ctx);
+
+        ca_get_acb(ekax, x, 53, ctx);
+        acb_elliptic_k(ekax, ekax, 53);
+        ca_get_acb(eeax, x, 53, ctx);
+        acb_elliptic_e(eeax, eeax, 53);
+        ca_get_acb(epaxy, x, 53, ctx);
+        ca_get_acb(ay, y, 53, ctx);
+        acb_elliptic_pi(epaxy, epaxy, ay, 53);
+
+        if (!acb_overlaps(aekx, ekax) || !acb_overlaps(aeex, eeax) ||
+            !acb_overlaps(aepxy, epaxy))
+        {
+            flint_printf("FAIL (Test 1)\n");
+            flint_printf("x = "); ca_print(x, ctx); printf("\n\n");
+            flint_printf("y = "); ca_print(y, ctx); printf("\n\n");
+            flint_printf("ekx = "); ca_print(ekx, ctx); printf("\n\n");
+            flint_printf("eex = "); ca_print(eex, ctx); printf("\n\n");
+            flint_printf("epxy = "); ca_print(epxy, ctx); printf("\n\n");
+            flint_abort();
+        }
+        
+        /* Test the Legendre relation. */
+        ca_si_sub(y, 1, x, ctx);
+        ca_t eky, eey, p;
+        acb_t ekay, eeay, ap, az;
+        ca_init(eky, ctx);
+        ca_init(eey, ctx);
+        ca_init(p, ctx);
+        acb_init(ekay);
+        acb_init(eeay);
+        acb_init(ap);
+        acb_init(az);
+        acb_zero(az);
+        
+        ca_pi(p, ctx);
+        ca_div_si(p, p, 2, ctx);
+        ca_get_acb(ap, p, 53, ctx);
+        ca_elliptic_k(eky, y, ctx);
+        ca_elliptic_e(eey, y, ctx);
+        
+        ca_mul(eey, eey, ekx, ctx);
+        ca_mul(eex, eex, eky, ctx);
+        ca_mul(ekx, ekx, eky, ctx);
+        ca_sub(p, p, eey, ctx);
+        ca_sub(p, p, eex, ctx);
+        ca_add(p, p, ekx, ctx);
+        ca_get_acb(ay, p, 53, ctx);
+        
+        ca_get_acb(ekay, y, 53, ctx);
+        acb_elliptic_k(ekay, ekay, 53);
+        ca_get_acb(eeay, y, 53, ctx);
+        acb_elliptic_e(eeay, eeay, 53);
+        
+        acb_mul(eeay, eeay, ekax, 53);
+        acb_mul(eeax, eeax, ekay, 53);
+        acb_mul(ekax, ekax, ekay, 53);
+        acb_sub(ap, ap, eeay, 53);
+        acb_sub(ap, ap, eeax, 53);
+        acb_add(ap, ap, ekax, 53);
+        
+        if (!acb_overlaps(ap, ay) || !acb_contains(ay, az))
+        {
+            flint_printf("FAIL (Test 2)\n");
+            flint_printf("x = "); ca_print(x, ctx); printf("\n\n");
+            flint_printf("y = "); ca_print(y, ctx); printf("\n\n");
+            flint_printf("ekx = "); ca_print(ekx, ctx); printf("\n\n");
+            flint_printf("eex = "); ca_print(eex, ctx); printf("\n\n");
+            flint_printf("eky = "); ca_print(eky, ctx); printf("\n\n");
+            flint_printf("eey = "); ca_print(eey, ctx); printf("\n\n");
+            flint_printf("ap = %{acb} \n", ap);
+            flint_printf("ay = %{acb} \n", ay);
+            
+            flint_abort();
+        }
+        
+        ca_clear(x, ctx);
+        ca_clear(y, ctx);
+        ca_clear(ekx, ctx);
+        ca_clear(eex, ctx);
+        ca_clear(eky, ctx);
+        ca_clear(eey, ctx);
+        ca_clear(p, ctx);
+        ca_clear(epxy, ctx);
+        acb_clear(ay);
+        acb_clear(ap);
+        acb_clear(az);
+        acb_clear(ekax);
+        acb_clear(eeax);
+        acb_clear(ekay);
+        acb_clear(eeay);
+        acb_clear(epaxy);
+        acb_clear(aekx);
+        acb_clear(aeex);
+        acb_clear(aepxy);
+        ca_ctx_clear(ctx);
+    }
+
+    TEST_FUNCTION_END(state);
+}
diff --git a/src/ca_ext/get_acb_raw.c b/src/ca_ext/get_acb_raw.c
index 3342f5b07e..415a3cf052 100644
--- a/src/ca_ext/get_acb_raw.c
+++ b/src/ca_ext/get_acb_raw.c
@@ -10,6 +10,7 @@
 */
 
 #include "acb_hypgeom.h"
+#include "acb_elliptic.h"
 #include "ca_ext.h"
 
 #define ARB_CONST(f) \
@@ -122,6 +123,10 @@ ca_ext_get_acb_raw(acb_t res, ca_ext_t x, slong prec, ca_ctx_t ctx)
         case CA_Erfi:         ACB_UNARY(acb_hypgeom_erfi)
         case CA_RiemannZeta:  ACB_UNARY(acb_zeta)
         case CA_HurwitzZeta:  ACB_BINARY(acb_hurwitz_zeta)
+        /* Complete elliptic integrals */
+        case CA_EllipticK:    ACB_UNARY(acb_elliptic_k)
+        case CA_EllipticE:    ACB_UNARY(acb_elliptic_e)
+        case CA_EllipticPi:   ACB_BINARY(acb_elliptic_pi)
         default:
             flint_throw(FLINT_ERROR, "ca_ext_get_acb_raw: unknown function\n");
     }
diff --git a/src/ca_field.h b/src/ca_field.h
index c6e8e1b09e..6989b71430 100644
--- a/src/ca_field.h
+++ b/src/ca_field.h
@@ -41,6 +41,7 @@ int ca_field_cmp(const ca_field_t K1, const ca_field_t K2, ca_ctx_t ctx);
 void ca_field_build_ideal(ca_field_t K, ca_ctx_t ctx);
 void ca_field_build_ideal_erf(ca_field_t K, ca_ctx_t ctx);
 void ca_field_build_ideal_gamma(ca_field_t K, ca_ctx_t ctx);
+void ca_field_build_ideal_elliptic(ca_field_t K, ca_ctx_t ctx);
 
 void ca_field_cache_init(ca_field_cache_t cache, ca_ctx_t ctx);
 void ca_field_cache_clear(ca_field_cache_t cache, ca_ctx_t ctx);
diff --git a/src/ca_field/build_ideal.c b/src/ca_field/build_ideal.c
index e12ff12a32..30c708f82e 100644
--- a/src/ca_field/build_ideal.c
+++ b/src/ca_field/build_ideal.c
@@ -1111,6 +1111,7 @@ ca_field_build_ideal(ca_field_t K, ca_ctx_t ctx)
     /* ca_field_build_ideal_sin_cos(K, ctx); */
     ca_field_build_ideal_erf(K, ctx);
     ca_field_build_ideal_gamma(K, ctx);
+    ca_field_build_ideal_elliptic(K, ctx);
 
     if (ctx->options[CA_OPT_USE_GROEBNER])
     {
diff --git a/src/ca_field/build_ideal_elliptic.c b/src/ca_field/build_ideal_elliptic.c
new file mode 100644
index 0000000000..5186acaf64
--- /dev/null
+++ b/src/ca_field/build_ideal_elliptic.c
@@ -0,0 +1,156 @@
+/*
+    Copyright (C) 2020 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "gr.h"
+#include "ca.h"
+#include "ca_ext.h"
+#include "ca_field.h"
+#include "fmpz_mpoly.h"
+
+void _ca_field_ideal_insert_clear_mpoly(ca_field_t K, fmpz_mpoly_t poly, fmpz_mpoly_ctx_t mctx, ca_ctx_t ctx);
+
+/* This merely implements the Legendre relation between the generators i, j (elliptic K-functions), k, l (elliptic E-functions) and pi_index (pi), if it exists. */
+void
+ca_field_build_elliptic_relation(ca_field_t K, slong i, slong j, slong k, slong l, slong pi_index, ca_ctx_t ctx)
+{
+    ca_ptr xi, xj, xk, xl;
+    ca_t t, u;
+    truth_t e;
+    
+    ca_init(t, ctx);
+    ca_init(u, ctx);
+    
+    xi = CA_EXT_FUNC_ARGS(CA_FIELD_EXT_ELEM(K, i));
+    xj = CA_EXT_FUNC_ARGS(CA_FIELD_EXT_ELEM(K, j));
+    xk = CA_EXT_FUNC_ARGS(CA_FIELD_EXT_ELEM(K, k));
+    xl = CA_EXT_FUNC_ARGS(CA_FIELD_EXT_ELEM(K, l));
+    
+    ca_add(t, xi, xj, ctx);
+    ca_add(u, xk, xl, ctx);
+    
+    /* swap indices if necessary */
+    if (ca_check_equal(xi, xl, ctx) == T_TRUE){
+        ca_swap(xk, xl, ctx);
+        slong d = k;
+        k = l;
+        l = d;
+    }
+    
+    e = truth_and(truth_and(ca_check_equal(xi, xk, ctx),
+        ca_check_is_one(t, ctx)), ca_check_is_one(u, ctx));
+    
+    if (e == T_TRUE)
+    {
+        fmpz_mpoly_t poly;
+        fmpz_mpoly_init(poly, CA_FIELD_MCTX(K, ctx));
+        
+        slong len = CA_FIELD_LENGTH(K);
+        ulong* exp = (ulong*) GR_TMP_ALLOC(len * sizeof(ulong));
+        for (slong v = 0; v < len; v++){
+            exp[v] = 0;
+        }
+        
+        exp[i] = 1;
+        exp[l] = 1;
+        fmpz_mpoly_set_coeff_si_ui(poly, 2, exp, CA_FIELD_MCTX(K, ctx));
+        exp[i] = 0;
+        exp[l] = 0;
+        exp[j] = 1;
+        exp[k] = 1;
+        fmpz_mpoly_set_coeff_si_ui(poly, 2, exp, CA_FIELD_MCTX(K, ctx));
+        exp[k] = 0;
+        exp[i] = 1;
+        fmpz_mpoly_set_coeff_si_ui(poly, -2, exp, CA_FIELD_MCTX(K, ctx));
+        exp[i] = 0;
+        exp[j] = 0;
+        exp[pi_index] = 1;
+        fmpz_mpoly_set_coeff_si_ui(poly, -1, exp, CA_FIELD_MCTX(K, ctx));
+        
+        _ca_field_ideal_insert_clear_mpoly(K, poly, CA_FIELD_MCTX(K, ctx), ctx);
+        
+        GR_TMP_FREE(exp, len * sizeof(ulong));
+    }
+    
+    ca_clear(t, ctx);
+    ca_clear(u, ctx);
+}
+
+void
+ca_field_build_ideal_elliptic(ca_field_t K, ca_ctx_t ctx)
+{
+    calcium_func_code Fi, Fj, Fk, Fl;
+    slong i, j, k, l, len, num_ellK, num_ellE, pi_index;
+
+    len = CA_FIELD_LENGTH(K);
+
+    if (len < 5)
+    {
+        return;
+    }
+
+    num_ellK = 0;
+    num_ellE = 0;
+    pi_index = -1;
+    for (i = 0; i < len; i++)
+    {
+        Fi = CA_EXT_HEAD(CA_FIELD_EXT_ELEM(K, i));
+
+        if (Fi == CA_EllipticK)
+        {
+            num_ellK++;
+        }
+        else if (Fi == CA_EllipticE)
+        {
+            num_ellE++;
+        }
+        else if (Fi == CA_Pi)
+        {
+            pi_index = i;
+        }
+    }
+
+    if (num_ellK < 2 || num_ellE < 2 || pi_index == -1)
+        return;
+
+    for (i = 0; i < len; i++)
+    {
+        Fi = CA_EXT_HEAD(CA_FIELD_EXT_ELEM(K, i));
+
+        if (Fi == CA_EllipticK)
+        {
+            for (j = i + 1; j < len; j++)
+            {
+                Fj = CA_EXT_HEAD(CA_FIELD_EXT_ELEM(K, j));
+
+                if (Fj == CA_EllipticK)
+                {
+                    for (k = 0; k < len; k++)
+                    {
+                        Fk = CA_EXT_HEAD(CA_FIELD_EXT_ELEM(K, k));
+
+                        if (Fk == CA_EllipticE)
+                        {
+                            for (l = k + 1; l < len; l++)
+                            {
+                                Fl = CA_EXT_HEAD(CA_FIELD_EXT_ELEM(K, l));
+
+                                if (Fl == CA_EllipticE)
+                                {
+                                    ca_field_build_elliptic_relation(K, i, j, k, l, pi_index, ctx);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
diff --git a/src/ca_types.h b/src/ca_types.h
index 33a33d70fa..bb3fff91cc 100644
--- a/src/ca_types.h
+++ b/src/ca_types.h
@@ -87,6 +87,10 @@ typedef enum
     CA_Erfi,
     CA_RiemannZeta,
     CA_HurwitzZeta,
+    /* Complete elliptic integrals */
+    CA_EllipticK,
+    CA_EllipticE,
+    CA_EllipticPi,
     CA_FUNC_CODE_LENGTH
 } calcium_func_code;
 
diff --git a/src/calcium/func_name.c b/src/calcium/func_name.c
index 44541f7ff3..72bf7e30d5 100644
--- a/src/calcium/func_name.c
+++ b/src/calcium/func_name.c
@@ -63,6 +63,10 @@ const char * calcium_func_name(calcium_func_code func)
         case CA_Erfi: return "Erfi";
         case CA_RiemannZeta: return "RiemannZeta";
         case CA_HurwitzZeta: return "HurwitzZeta";
+        /* Complete elliptic integrals */
+        case CA_EllipticK: return "EllipticK";
+        case CA_EllipticE: return "EllipticE";
+        case CA_EllipticPi: return "EllipticPi";
         default: return "<unknown function>";
     }
 }