diff --git a/acb_elliptic.h b/acb_elliptic.h
index a3227a2a..626a00f4 100644
--- a/acb_elliptic.h
+++ b/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
diff --git a/acb_elliptic/rc1.c b/acb_elliptic/rc1.c
new file mode 100644
index 00000000..b48d3ea1
--- /dev/null
+++ b/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 .
+*/
+
+#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);
+}
+
diff --git a/acb_elliptic/rj.c b/acb_elliptic/rj.c
new file mode 100644
index 00000000..ce3e08ac
--- /dev/null
+++ b/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 .
+*/
+
+#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);
+}
+
diff --git a/acb_elliptic/test/t-rc1.c b/acb_elliptic/test/t-rc1.c
new file mode 100644
index 00000000..e32cee96
--- /dev/null
+++ b/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 .
+*/
+
+#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;
+}
+
diff --git a/acb_elliptic/test/t-rj.c b/acb_elliptic/test/t-rj.c
new file mode 100644
index 00000000..dc769b64
--- /dev/null
+++ b/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 .
+*/
+
+#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;
+}
+
diff --git a/doc/source/acb_elliptic.rst b/doc/source/acb_elliptic.rst
index 06c08a94..2505f3ff 100644
--- a/doc/source/acb_elliptic.rst
+++ b/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`.
+
+