Add Calcium support for complete elliptic integrals

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
 -------------------------------------------------------------------------------

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);

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));
+}
+

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),

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);
+}

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");
     }