Carlson elliptic integral of the third kind

acb_elliptic.h

@@ -21,6 +21,10 @@ extern "C" {
 
 void acb_elliptic_rf(acb_t res, const acb_t x, const acb_t y, const acb_t z, int flags, slong prec);
 
+void acb_elliptic_rj(acb_t res, const acb_t x, const acb_t y, const acb_t z, const acb_t p, int flags, slong prec);
+
+void acb_elliptic_rc1(acb_t res, const acb_t x, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

acb_elliptic/rc1.c

@@ -0,0 +1,114 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+
+/* Compute R_C(1,1+x) = atan(sqrt(x))/sqrt(x) = 2F1(1,1/2,3/2,-x) */
+static void
+_acb_elliptic_rc1(acb_t res, const acb_t x, slong prec)
+{
+    acb_t t;
+    acb_init(t);
+    acb_sqrt(t, x, prec + 2);
+    acb_atan(res, t, prec + 2);
+    acb_div(res, res, t, prec);
+    acb_clear(t);
+}
+
+void
+acb_elliptic_rc1(acb_t res, const acb_t x, slong prec)
+{
+    mag_t xm;
+
+    mag_init(xm);
+    acb_get_mag(xm, x);
+
+    if (mag_cmp_2exp_si(xm, 0) < 0)
+    {
+        slong k, n;
+
+        for (n = 1; n <= 6; n++)
+        {
+            if (mag_cmp_2exp_si(xm, -prec / n) < 0)
+                break;
+        }
+
+        /* Use Taylor series: 1 - x/3 + x^2/5 - x^3/7 + x^4/9 + ... */
+        if (n <= 6)
+        {
+            const short coeffs[] = {3465, -1155, 693, -495, 385, -315};
+
+            acb_t s;
+            acb_init(s);
+
+            for (k = n - 1; k >= 0; k--)
+            {
+                acb_mul(s, s, x, prec);
+                acb_add_si(s, s, coeffs[k], prec);
+            }
+
+            acb_div_si(s, s, coeffs[0], prec);
+            mag_geom_series(xm, xm, n);
+            if (acb_is_real(x))
+                arb_add_error_mag(acb_realref(s), xm);
+            else
+                acb_add_error_mag(s, xm);
+
+            acb_set(res, s);
+            acb_clear(s);
+        }
+        else if (acb_is_exact(x))
+        {
+            _acb_elliptic_rc1(res, x, prec);
+        }
+        else
+        {
+            acb_t w;
+            mag_t err, rad;
+
+            acb_init(w);
+            mag_init(err);
+            mag_init(rad);
+
+            /* On the unit disc, |f'(x)| <= 0.5 / |1+x| */
+            acb_add_ui(w, x, 1, MAG_BITS);
+            acb_get_mag_lower(err, w);
+            mag_one(rad);
+            mag_mul_2exp_si(rad, rad, -1);
+            mag_div(err, rad, err);
+            mag_hypot(rad, arb_radref(acb_realref(x)), arb_radref(acb_imagref(x)));
+            mag_mul(err, err, rad);
+
+            acb_set(w, x);
+            mag_zero(arb_radref(acb_realref(w)));
+            mag_zero(arb_radref(acb_imagref(w)));
+            _acb_elliptic_rc1(w, w, prec);
+
+            if (acb_is_real(x))
+                arb_add_error_mag(acb_realref(w), err);
+            else
+                acb_add_error_mag(w, err);
+
+            acb_set(res, w);
+
+            acb_clear(w);
+            mag_clear(err);
+            mag_clear(rad);
+        }
+    }
+    else
+    {
+        _acb_elliptic_rc1(res, x, prec);
+    }
+
+    mag_clear(xm);
+}
+

acb_elliptic/rj.c

@@ -0,0 +1,270 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+
+void
+acb_elliptic_rj(acb_t res, const acb_t x, const acb_t y,
+            const acb_t z, const acb_t p, int flags, slong prec)
+{
+    acb_t xx, yy, zz, pp, sx, sy, sz, sp, t, d, delta, S;
+    acb_t A, AA, X, Y, Z, P, E2, E3, E4, E5;
+    mag_t err, err2;
+    slong k;
+    int rd;
+
+    if (!acb_is_finite(x) || !acb_is_finite(y) || !acb_is_finite(z) ||
+        !acb_is_finite(p))
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    if ((acb_contains_zero(x) + acb_contains_zero(y) + acb_contains_zero(z) > 1)
+        || acb_contains_zero(p))
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    /* Special case computing R_D(x,y,z) */
+    rd = (z == p);
+
+    acb_init(xx); acb_init(yy); acb_init(zz); acb_init(pp);
+    acb_init(sx); acb_init(sy); acb_init(sz); acb_init(sp);
+    acb_init(S); acb_init(A); acb_init(AA);
+    acb_init(X); acb_init(Y); acb_init(Z); acb_init(P);
+    acb_init(E2); acb_init(E3); acb_init(E4); acb_init(E5);
+    acb_init(t); acb_init(d); acb_init(delta);
+    mag_init(err);
+    mag_init(err2);
+
+    acb_set(xx, x);
+    acb_set(yy, y);
+    acb_set(zz, z);
+    acb_set(pp, p);
+    acb_zero(S);
+
+    if (!rd)
+    {
+        acb_mul_2exp_si(A, p, 1);
+        acb_add(A, A, z, prec);
+    }
+    else
+    {
+        acb_mul_ui(A, z, 3, prec);
+    }
+    acb_add(A, A, x, prec);
+    acb_add(A, A, y, prec);
+    acb_div_ui(A, A, 5, prec);
+    acb_set(AA, A);
+
+    if (!rd)
+    {
+        acb_sub(delta, p, x, prec);
+        acb_sub(t, p, y, prec);
+        acb_mul(delta, delta, t, prec);
+        acb_sub(t, p, z, prec);
+        acb_mul(delta, delta, t, prec);
+    }
+
+    /* must do at least one iteration */
+    for (k = 0; k < prec; k++)
+    {
+        acb_sqrt(sx, xx, prec);
+        acb_sqrt(sy, yy, prec);
+        acb_sqrt(sz, zz, prec);
+        if (!rd) acb_sqrt(sp, pp, prec);
+
+        acb_add(t, sy, sz, prec);
+        acb_mul(t, t, sx, prec);
+        acb_addmul(t, sy, sz, prec);
+
+        acb_add(xx, xx, t, prec);
+        acb_add(yy, yy, t, prec);
+        acb_add(zz, zz, t, prec);
+        if (!rd) acb_add(pp, pp, t, prec);
+        acb_add(AA, AA, t, prec);
+
+        acb_mul_2exp_si(xx, xx, -2);
+        acb_mul_2exp_si(yy, yy, -2);
+        acb_mul_2exp_si(zz, zz, -2);
+        if (!rd) acb_mul_2exp_si(pp, pp, -2);
+        acb_mul_2exp_si(AA, AA, -2);
+
+        if (!rd)
+        {
+            /* d = (sp+sx)(sp+sy)(sp+sz) */
+            /* e = 4^(-3k) delta / d^2 */
+            /* S += 4^(-k) RC(1, 1+e) / d */
+            acb_add(d, sp, sx, prec);
+            acb_add(t, sp, sy, prec);
+            acb_mul(d, d, t, prec);
+            acb_add(t, sp, sz, prec);
+            acb_mul(d, d, t, prec);
+
+            /* E2 = e */
+            acb_mul(E2, d, d, prec);
+            acb_div(E2, delta, E2, prec);
+            acb_mul_2exp_si(E2, E2, -6 * k);
+
+            acb_elliptic_rc1(E4, E2, prec);
+            acb_div(E4, E4, d, prec);
+            acb_mul_2exp_si(E4, E4, -2 * k);
+
+            acb_add(S, S, E4, prec);
+        }
+        else
+        {
+            acb_mul(t, sz, zz, prec);
+            acb_mul_2exp_si(t, t, 2);
+            acb_inv(t, t, prec);
+            acb_mul_2exp_si(t, t, -2 * k);
+            acb_mul_2exp_si(t, t, -1);
+
+            acb_add(S, S, t, prec);
+        }
+
+        /* Close enough? */
+        acb_sub(t, xx, yy, prec);
+        acb_get_mag(err, t);
+        acb_sub(t, xx, zz, prec);
+        acb_get_mag(err2, t);
+        mag_max(err, err, err2);
+        if (!rd)
+        {
+            acb_sub(t, xx, pp, prec);
+            acb_get_mag(err2, t);
+            mag_max(err, err, err2);
+        }
+        acb_get_mag_lower(err2, xx);
+        mag_div(err, err, err2);
+
+        mag_pow_ui(err, err, 8);
+
+        if (mag_cmp_2exp_si(err, -prec) < 0)
+        {
+            k++;
+            break;
+        }
+    }
+
+    /* X = (A-x)/(4^k AA) */
+    /* Y = (A-y)/(4^k AA) */
+    /* Z = (A-z)/(4^k AA) */
+    /* P = (-X-Y-Z)/2 */
+    acb_mul_2exp_si(t, AA, 2 * k);
+    acb_inv(t, t, prec);
+    acb_sub(X, A, x, prec);
+    acb_mul(X, X, t, prec);
+    acb_sub(Y, A, y, prec);
+    acb_mul(Y, Y, t, prec);
+    acb_sub(Z, A, z, prec);
+    acb_mul(Z, Z, t, prec);
+    acb_add(P, X, Y, prec);
+    acb_add(P, P, Z, prec);
+    acb_neg(P, P);
+    acb_mul_2exp_si(P, P, -1);
+
+    /* todo: improve for R_D */
+    /* E2 = XY + XZ + YZ - 3 P^2 */
+    /* E3 = XYZ + 2 E2 P + 4 P^3 */
+    /* E4 = (2 XYZ + E2 P + 3 P^3) P */
+    /* E5 = XYZP^2 */
+
+    acb_mul(t, P, P, prec); /* t = P^2 */
+
+    acb_mul(E2, X, Y, prec);
+    acb_mul(E3, E2, Z, prec);
+    acb_mul_2exp_si(E4, E3, 1);
+    acb_mul(E5, E3, t, prec);
+
+    acb_add(sx, X, Y, prec);
+    acb_addmul(E2, sx, Z, prec);
+    acb_submul_ui(E2, t, 3, prec);
+
+    acb_mul(sx, E2, P, prec);
+    acb_add(E4, E4, sx, prec);
+    acb_mul_2exp_si(sx, sx, 1);
+    acb_add(E3, E3, sx, prec);
+
+    acb_mul(t, t, P, prec); /* t = P^3 */
+    acb_addmul_ui(E3, t, 4, prec);
+    acb_addmul_ui(E4, t, 3, prec);
+    acb_mul(E4, E4, P, prec);
+
+    /* Error bound. */
+    acb_get_mag(err, X);
+    acb_get_mag(err2, Y);
+    mag_max(err, err, err2);
+    acb_get_mag(err2, Z);
+    mag_max(err, err, err2);
+    acb_get_mag(err2, P);
+    mag_max(err, err, err2);
+    mag_mul_ui(err, err, 9);
+    mag_mul_2exp_si(err, err, -3);
+    mag_geom_series(err, err, 8);
+    mag_mul_2exp_si(err, err, 1);
+
+    /*
+    1 + (-255255*E2**3 + 675675*E2**2*E3 + 417690*E2**2 - 706860*E2*E3
+        + 612612*E2*E4 - 540540*E2*E5 - 875160*E2 + 306306*E3**2 - 540540*E3*E4
+        + 680680*E3 - 556920*E4 + 471240*E5)/4084080
+    =
+    1 + (E2*(E2*(-255255*E2 + 675675*E3 + 417690) - 706860*E3
+        + 612612*E4 - 540540*E5 - 875160) + E3*(306306*E3 - 540540*E4
+        + 680680) - 556920*E4 + 471240*E5)/4084080
+    */
+    acb_set_ui(sx, 417690);
+    acb_addmul_ui(sx, E3, 675675, prec);
+    acb_submul_ui(sx, E2, 255255, prec);
+    acb_mul(sx, sx, E2, prec);
+    acb_submul_ui(sx, E3, 706860, prec);
+    acb_addmul_ui(sx, E4, 612612, prec);
+    acb_submul_ui(sx, E5, 540540, prec);
+    acb_sub_ui(sx, sx, 875160, prec);
+    acb_mul(sx, sx, E2, prec);
+
+    acb_set_ui(sy, 680680);
+    acb_submul_ui(sy, E4, 540540, prec);
+    acb_addmul_ui(sy, E3, 306306, prec);
+    acb_addmul(sx, sy, E3, prec);
+
+    acb_addmul_ui(sx, E5, 471240, prec);
+    acb_submul_ui(sx, E4, 556920, prec);
+    acb_div_ui(sx, sx, 4084080, prec);
+
+    acb_add_ui(sx, sx, 1, prec);
+
+    if (acb_is_real(X) && acb_is_real(Y) && acb_is_real(Z))
+        arb_add_error_mag(acb_realref(sx), err);
+    else
+        acb_add_error_mag(sx, err);
+
+    acb_rsqrt(t, AA, prec);
+    acb_div(t, t, AA, prec);
+    acb_mul_2exp_si(t, t, -2 * k);
+    acb_mul(t, t, sx, prec);
+
+    acb_addmul_ui(t, S, 6, prec);
+
+    acb_set(res, t);
+
+    acb_clear(xx); acb_clear(yy); acb_clear(zz); acb_clear(pp);
+    acb_clear(sx); acb_clear(sy); acb_clear(sz); acb_clear(sp);
+    acb_clear(S); acb_clear(A); acb_clear(AA);
+    acb_clear(X); acb_clear(Y); acb_clear(Z); acb_clear(P);
+    acb_clear(E2); acb_clear(E3); acb_clear(E4); acb_clear(E5);
+    acb_clear(t); acb_clear(d); acb_clear(delta);
+    mag_clear(err);
+    mag_clear(err2);
+}
+

acb_elliptic/test/t-rc1.c

@@ -0,0 +1,73 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+#include "acb_hypgeom.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("rc1....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
+    {
+        acb_t x, r1, r2, a, b, c;
+
+        acb_init(x);
+        acb_init(r1);
+        acb_init(r2);
+        acb_init(a);
+        acb_init(b);
+        acb_init(c);
+
+        acb_randtest(x, state, 1 + n_randint(state, 300), 10);
+        acb_randtest(a, state, 1 + n_randint(state, 300), 10);
+        acb_add(r2, x, a, 200);
+        acb_sub(r2, r2, a, 200);
+        acb_neg(r2, r2);
+        if (n_randint(state, 2))
+            acb_swap(x, r2);
+
+        acb_elliptic_rc1(r1, x, 2 + n_randint(state, 200));
+
+        acb_one(a);
+        acb_set_d(b, 0.5);
+        acb_set_d(c, 1.5);
+        acb_hypgeom_2f1(r2, a, b, c, r2, 0, 20 + n_randint(state, 200));
+
+        if (!acb_overlaps(r1, r2))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("x = "); acb_printd(x, 30); flint_printf("\n\n");
+            flint_printf("r1 = "); acb_printd(r1, 30); flint_printf("\n\n");
+            flint_printf("r2 = "); acb_printd(r2, 30); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_clear(x);
+        acb_clear(r1);
+        acb_clear(r2);
+        acb_clear(a);
+        acb_clear(b);
+        acb_clear(c);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_elliptic/test/t-rj.c

@@ -0,0 +1,192 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_elliptic.h"
+
+/* Test input from Carlson's paper and checked with mpmath. */
+
+static const double testdata_rj[8][10] = {
+    {0.0, 0.0, 1.0, 0.0, 2.0, 0.0, 3.0, 0.0, 0.77688623778582332014, 0.0},
+    {2.0, 0.0, 3.0, 0.0, 4.0, 0.0, 5.0, 0.0, 0.14297579667156753833, 0.0},
+    {2.0, 0.0, 3.0, 0.0, 4.0, 0.0, -1.0, 1.0, 0.13613945827770535204, -0.3820756162442716425},
+    {0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 2.0, 0.0, 1.6490011662710884518, 0.0},
+    {-1.0, 1.0, -1.0, -1.0, 1.0, 0.0, 2.0, 0.0, 0.94148358841220238083, 0.0},
+    {0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 1.0, -1.0, 1.8260115229009316249, 1.22906619086434715},
+    {-1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -3.0, 1.0, -0.61127970812028172124, -1.068403839000680788},
+    {-1.0, 1.0, -2.0, -1.0, 0.0, -1.0, -1.0, 1.0, 1.8249027393703805305, -1.2218475784827035855},
+};
+
+static const double testdata_rd[6][8] = {
+    {0.0, 0.0, 2.0, 0.0, 1.0, 0.0, 1.7972103521033883112, 0.0},
+    {2.0, 0.0, 3.0, 0.0, 4.0, 0.0, 0.16510527294261053349, 0.0},
+    {0.0, 1.0, 0.0, -1.0, 2.0, 0.0, 0.65933854154219768919, 0.0},
+    {0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 1.2708196271909686299, 2.7811120159520578776},
+    {0.0, 0.0, -1.0, 1.0, 0.0, 1.0, -1.8577235439239060056, -0.96193450888838559989},
+    {-2.0, -1.0, 0.0, -1.0, -1.0, 1.0, 1.8249027393703805305, -1.2218475784827035855},
+};
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("rj....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t x, y, z, p, r1, r2;
+        slong prec1, prec2;
+
+        prec1 = 2 + n_randint(state, 300);
+        prec2 = 2 + n_randint(state, 300);
+
+        acb_init(x);
+        acb_init(y);
+        acb_init(z);
+        acb_init(p);
+        acb_init(r1);
+        acb_init(r2);
+
+        if (iter == 0)
+        {
+            slong k;
+
+            for (k = 0; k < 8; k++)
+            {
+                acb_set_d_d(x, testdata_rj[k][0], testdata_rj[k][1]);
+                acb_set_d_d(y, testdata_rj[k][2], testdata_rj[k][3]);
+                acb_set_d_d(z, testdata_rj[k][4], testdata_rj[k][5]);
+                acb_set_d_d(p, testdata_rj[k][6], testdata_rj[k][7]);
+                acb_set_d_d(r2, testdata_rj[k][8], testdata_rj[k][9]);
+                mag_set_d(arb_radref(acb_realref(r2)), 1e-14 * fabs(testdata_rj[k][8]));
+                mag_set_d(arb_radref(acb_imagref(r2)), 1e-14 * fabs(testdata_rj[k][9]));
+
+                for (prec1 = 16; prec1 <= 256; prec1 *= 2)
+                {
+                    acb_elliptic_rj(r1, x, y, z, p, 0, prec1);
+
+                    if (!acb_overlaps(r1, r2) || acb_rel_accuracy_bits(r1) < prec1 - 12)
+                    {
+                        flint_printf("FAIL: overlap (testdata rj)\n\n");
+                        flint_printf("prec = %wd, accuracy = %wd\n\n", prec1, acb_rel_accuracy_bits(r1));
+                        flint_printf("x = "); acb_printd(x, 30); flint_printf("\n\n");
+                        flint_printf("y = "); acb_printd(y, 30); flint_printf("\n\n");
+                        flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
+                        flint_printf("p = "); acb_printd(p, 30); flint_printf("\n\n");
+                        flint_printf("r1 = "); acb_printd(r1, 30); flint_printf("\n\n");
+                        flint_printf("r2 = "); acb_printd(r2, 30); flint_printf("\n\n");
+                        abort();
+                    }
+                }
+            }
+
+            for (k = 0; k < 6; k++)
+            {
+                acb_set_d_d(x, testdata_rd[k][0], testdata_rd[k][1]);
+                acb_set_d_d(y, testdata_rd[k][2], testdata_rd[k][3]);
+                acb_set_d_d(z, testdata_rd[k][4], testdata_rd[k][5]);
+                acb_set_d_d(r2, testdata_rd[k][6], testdata_rd[k][7]);
+                mag_set_d(arb_radref(acb_realref(r2)), 1e-14 * fabs(testdata_rd[k][6]));
+                mag_set_d(arb_radref(acb_imagref(r2)), 1e-14 * fabs(testdata_rd[k][7]));
+
+                for (prec1 = 16; prec1 <= 256; prec1 *= 2)
+                {
+                    acb_elliptic_rj(r1, x, y, z, z, 0, prec1);
+
+                    if (!acb_overlaps(r1, r2) || acb_rel_accuracy_bits(r1) < prec1 - 12)
+                    {
+                        flint_printf("FAIL: overlap (testdata rd)\n\n");
+                        flint_printf("prec = %wd, accuracy = %wd\n\n", prec1, acb_rel_accuracy_bits(r1));
+                        flint_printf("x = "); acb_printd(x, 30); flint_printf("\n\n");
+                        flint_printf("y = "); acb_printd(y, 30); flint_printf("\n\n");
+                        flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
+                        flint_printf("r1 = "); acb_printd(r1, 30); flint_printf("\n\n");
+                        flint_printf("r2 = "); acb_printd(r2, 30); flint_printf("\n\n");
+                        abort();
+                    }
+                }
+            }
+        }
+
+        acb_randtest(x, state, 1 + n_randint(state, 300), 1 + n_randint(state, 30));
+        acb_randtest(y, state, 1 + n_randint(state, 300), 1 + n_randint(state, 30));
+        acb_randtest(z, state, 1 + n_randint(state, 300), 1 + n_randint(state, 30));
+        acb_randtest(p, state, 1 + n_randint(state, 300), 1 + n_randint(state, 30));
+
+        acb_elliptic_rj(r1, x, y, z, p, 0, prec1);
+
+        switch (n_randint(state, 6))
+        {
+            case 0:
+                acb_elliptic_rj(r2, x, y, z, p, 0, prec2);
+                break;
+            case 1:
+                acb_elliptic_rj(r2, x, z, y, p, 0, prec2);
+                break;
+            case 2:
+                acb_elliptic_rj(r2, y, x, z, p, 0, prec2);
+                break;
+            case 3:
+                acb_elliptic_rj(r2, y, z, x, p, 0, prec2);
+                break;
+            case 4:
+                acb_elliptic_rj(r2, z, x, y, p, 0, prec2);
+                break;
+            default:
+                acb_elliptic_rj(r2, z, y, x, p, 0, prec2);
+                break;
+        }
+
+        if (!acb_overlaps(r1, r2))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("x = "); acb_printd(x, 30); flint_printf("\n\n");
+            flint_printf("y = "); acb_printd(y, 30); flint_printf("\n\n");
+            flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
+            flint_printf("p = "); acb_printd(p, 30); flint_printf("\n\n");
+            flint_printf("r1 = "); acb_printd(r1, 30); flint_printf("\n\n");
+            flint_printf("r2 = "); acb_printd(r2, 30); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_elliptic_rj(r1, x, y, z, z, 0, prec1);
+        acb_set(p, z);
+        acb_elliptic_rj(r2, x, y, z, p, 0, prec2);
+
+        if (!acb_overlaps(r1, r2))
+        {
+            flint_printf("FAIL: overlap R_D\n\n");
+            flint_printf("x = "); acb_printd(x, 30); flint_printf("\n\n");
+            flint_printf("y = "); acb_printd(y, 30); flint_printf("\n\n");
+            flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n");
+            flint_printf("p = "); acb_printd(p, 30); flint_printf("\n\n");
+            flint_printf("r1 = "); acb_printd(r1, 30); flint_printf("\n\n");
+            flint_printf("r2 = "); acb_printd(r2, 30); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_clear(x);
+        acb_clear(y);
+        acb_clear(z);
+        acb_clear(p);
+        acb_clear(r1);
+        acb_clear(r2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_elliptic.rst

@@ -3,12 +3,22 @@
 **acb_elliptic.h** -- elliptic integrals and functions of complex variables
 ===============================================================================
 
+Warning: incomplete elliptic integrals have very complicated
+branch structure when extended to complex variables.
+For some functions in this module, branch cuts may be
+artifacts of the evaluation algorithm rather than having
+a natural mathematical justification.
+The user should, accordingly, watch out for edge cases where the functions
+implemented here may differ from other systems or literature.
+There may also exist points where a function should be well-defined
+but the implemented algorithm
+fails to produce a finite result due to artificial internal singularities.
 
 Carlson symmetric elliptic integrals
 -------------------------------------------------------------------------------
 
 Carlson symmetric forms are the preferred form of incomplete elliptic
-integrals, due to their neat analytic structure and relatively
+integrals, due to their neat properties and relatively
 simple computation based on duplication theorems.
 We largely follow the definitions and algorithms
 in [Car1995]_ and chapter 19 in [NIST2012]_.
@@ -22,7 +32,7 @@ in [Car1995]_ and chapter 19 in [NIST2012]_.
         R_F(x,y,z) = \frac{1}{2}
             \int_0^{\infty} \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}
 
-    where the square root extends continuously from the positive real axis.
+    where the square root extends continuously from positive infinity.
     The function is well-defined for `x,y,z \notin (-\infty,0)`, and with
     at most one of `x,y,z` being zero.
 
@@ -33,6 +43,48 @@ in [Car1995]_ and chapter 19 in [NIST2012]_.
     The special case `R_C(x, y) = R_F(x, y, y) = \frac{1}{2} \int_0^{\infty} (t+x)^{-1/2} (t+y)^{-1} dt`
     may be computed by
     setting *y* and *z* to the same variable.
+    (This case is not yet handled specially, but might be optimized in
+    the future.)
 
     The *flags* parameter is reserved for future use and currently
     does nothing. Passing 0 results in default behavior.
+
+.. function:: void acb_elliptic_rj(acb_t res, const acb_t x, const acb_t y, const acb_t z, const acb_t p, int flags, slong prec)
+
+    Evaluates the Carlson symmetric elliptic integral of the third kind
+
+    .. math ::
+
+        R_J(x,y,z,p) = \frac{3}{2}
+            \int_0^{\infty} \frac{dt}{(t+p)\sqrt{(t+x)(t+y)(t+z)}}.
+
+    where the square root is taken continuously as in `R_J`.
+
+    In general, one or more duplication steps are applied until
+    `x,y,z,p` are close enough to use a multivariate Taylor polynomial
+    of total degree 7.
+
+    The duplication algorithm might not be correct for all possible
+    combinations of complex variables, i.e. taking square roots
+    during the computation might introduce spurious branch cuts.
+    According to [Car1995]_, a sufficient (but not necessary) condition
+    for correctness is that *x*, *y*, *z* have nonnegative
+    real part and that *p* has positive real part.
+    In other cases, the algorithm *may* still be correct, but the user
+    should verify the results.
+
+    The special case `R_D(x, y, z) = R_J(x, y, z, z)`
+    may be computed by setting *y* and *z* to the same variable.
+    This case is handled specially to avoid redundant arithmetic operations.
+    In this case, the algorithm is correct for all *x*, *y* and *z*.
+
+    The *flags* parameter is reserved for future use and currently
+    does nothing. Passing 0 results in default behavior.
+
+.. function:: void acb_elliptic_rc1(acb_t res, const acb_t x, slong prec)
+
+    This helper function computes the special case
+    `R_C(1, 1+x) = \operatorname{atan}(\sqrt{x})/\sqrt{x} = {}_2F_1(1,1/2,3/2,-x)`,
+    which is needed in the evaluation of `R_J`.
+
+