Skip to content

Add Calcium support for complete elliptic integrals #2302

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions doc/source/ca.rst
Original file line number Diff line number Diff line change
@@ -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
-------------------------------------------------------------------------------
4 changes: 4 additions & 0 deletions src/ca.h
Original file line number Diff line number Diff line change
@@ -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);
156 changes: 156 additions & 0 deletions src/ca/elliptic.c
Original file line number Diff line number Diff line change
@@ -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));
}

2 changes: 2 additions & 0 deletions src/ca/test/main.c
Original file line number Diff line number Diff line change
@@ -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),
153 changes: 153 additions & 0 deletions src/ca/test/t-elliptic.c
Original file line number Diff line number Diff line change
@@ -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);
}
5 changes: 5 additions & 0 deletions src/ca_ext/get_acb_raw.c
Original file line number Diff line number Diff line change
@@ -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");
}
1 change: 1 addition & 0 deletions src/ca_field.h
Original file line number Diff line number Diff line change
@@ -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);
1 change: 1 addition & 0 deletions src/ca_field/build_ideal.c
Original file line number Diff line number Diff line change
@@ -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])
{
156 changes: 156 additions & 0 deletions src/ca_field/build_ideal_elliptic.c
Original file line number Diff line number Diff line change
@@ -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);
}
}
}
}
}
}
}
}
}
4 changes: 4 additions & 0 deletions src/ca_types.h
Original file line number Diff line number Diff line change
@@ -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;

4 changes: 4 additions & 0 deletions src/calcium/func_name.c
Original file line number Diff line number Diff line change
@@ -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>";
}
}