diff --git a/acb_elliptic.h b/acb_elliptic.h
index 626a00f4..de7d964e 100644
--- a/acb_elliptic.h
+++ b/acb_elliptic.h
@@ -23,6 +23,8 @@ void acb_elliptic_rf(acb_t res, const acb_t x, const acb_t y, const acb_t z, int
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_rg(acb_t res, const acb_t x, const acb_t y, const acb_t z, int flags, slong prec);
+
void acb_elliptic_rc1(acb_t res, const acb_t x, slong prec);
#ifdef __cplusplus
diff --git a/acb_elliptic/rg.c b/acb_elliptic/rg.c
new file mode 100644
index 00000000..35868336
--- /dev/null
+++ b/acb_elliptic/rg.c
@@ -0,0 +1,85 @@
+/*
+ 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_rg(acb_t res, const acb_t x, const acb_t y, const acb_t z,
+ int flags, slong prec)
+{
+ acb_t a, b, c, t;
+ slong wp;
+
+ wp = prec + 10;
+
+ acb_init(a);
+ acb_init(b);
+ acb_init(c);
+ acb_init(t);
+
+ acb_elliptic_rf(a, x, y, z, 0, wp);
+ acb_mul(a, a, z, wp);
+
+ acb_elliptic_rj(b, x, y, z, z, 0, wp);
+ acb_sub(c, x, z, wp);
+ acb_mul(b, b, c, wp);
+ acb_sub(c, z, y, wp);
+ acb_mul(b, b, c, wp);
+ acb_div_ui(b, b, 3, wp);
+
+ acb_sqrt(c, x, wp);
+ acb_sqrt(t, y, wp);
+ acb_mul(c, c, t, wp);
+ acb_rsqrt(t, z, wp);
+ acb_mul(c, c, t, wp);
+
+ acb_add(res, a, b, wp);
+ acb_add(res, res, c, prec);
+ acb_mul_2exp_si(res, res, -1);
+
+ acb_clear(a);
+ acb_clear(b);
+ acb_clear(c);
+ acb_clear(t);
+}
+
+void
+acb_elliptic_rg(acb_t res, const acb_t x, const acb_t y, const acb_t z,
+ int flags, slong prec)
+{
+ if (acb_is_zero(x) && acb_is_zero(y))
+ {
+ acb_sqrt(res, z, prec);
+ acb_mul_2exp_si(res, res, -1);
+ }
+ else if (acb_is_zero(x) && acb_is_zero(z))
+ {
+ acb_sqrt(res, y, prec);
+ acb_mul_2exp_si(res, res, -1);
+ }
+ else if (acb_is_zero(y) && acb_is_zero(z))
+ {
+ acb_sqrt(res, x, prec);
+ acb_mul_2exp_si(res, res, -1);
+ }
+ else if (acb_contains_zero(z))
+ {
+ if (acb_contains_zero(y))
+ _acb_elliptic_rg(res, y, z, x, flags, prec);
+ else
+ _acb_elliptic_rg(res, x, z, y, flags, prec);
+ }
+ else
+ {
+ _acb_elliptic_rg(res, x, y, z, flags, prec);
+ }
+}
+
diff --git a/acb_elliptic/test/t-rg.c b/acb_elliptic/test/t-rg.c
new file mode 100644
index 00000000..3eec2464
--- /dev/null
+++ b/acb_elliptic/test/t-rg.c
@@ -0,0 +1,134 @@
+/*
+ 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_rg[7][8] = {
+ {0.0, 0.0, 16.0, 0.0, 16.0, 0.0, 3.1415926535897932385, 0.0},
+ {2.0, 0.0, 3.0, 0.0, 4.0, 0.0, 1.7255030280692277601, 0.0},
+ {0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.4236065423969895433, 0.0},
+ {-1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.44660591677018372657, 0.70768352357515390073},
+ {0.0, -1.0, -1.0, 1.0, 0.0, 1.0, 0.36023392184473309034, 0.40348623401722113741},
+ {0.0, 0.0, 0.0796, 0.0, 4.0, 0.0, 1.0284758090288040022, 0.0},
+ /* more tests */
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+};
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("rg....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+ {
+ acb_t x, y, z, 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(r1);
+ acb_init(r2);
+
+ if (iter == 0)
+ {
+ slong k;
+
+ for (k = 0; k < 7; k++)
+ {
+ acb_set_d_d(x, testdata_rg[k][0], testdata_rg[k][1]);
+ acb_set_d_d(y, testdata_rg[k][2], testdata_rg[k][3]);
+ acb_set_d_d(z, testdata_rg[k][4], testdata_rg[k][5]);
+ acb_set_d_d(r2, testdata_rg[k][6], testdata_rg[k][7]);
+ mag_set_d(arb_radref(acb_realref(r2)), 1e-14 * fabs(testdata_rg[k][6]));
+ mag_set_d(arb_radref(acb_imagref(r2)), 1e-14 * fabs(testdata_rg[k][7]));
+
+ for (prec1 = 16; prec1 <= 256; prec1 *= 2)
+ {
+ acb_elliptic_rg(r1, x, y, z, 0, prec1);
+
+ if (!acb_overlaps(r1, r2) || acb_rel_accuracy_bits(r1) < prec1 * 0.9 - 10)
+ {
+ flint_printf("FAIL: overlap (testdata rg)\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_elliptic_rg(r1, x, y, z, 0, prec1);
+
+ switch (n_randint(state, 6))
+ {
+ case 0:
+ acb_elliptic_rg(r2, x, y, z, 0, prec2);
+ break;
+ case 1:
+ acb_elliptic_rg(r2, x, z, y, 0, prec2);
+ break;
+ case 2:
+ acb_elliptic_rg(r2, y, x, z, 0, prec2);
+ break;
+ case 3:
+ acb_elliptic_rg(r2, y, z, x, 0, prec2);
+ break;
+ case 4:
+ acb_elliptic_rg(r2, z, x, y, 0, prec2);
+ break;
+ default:
+ acb_elliptic_rg(r2, z, y, x, 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("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(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 2505f3ff..7b22e0ef 100644
--- a/doc/source/acb_elliptic.rst
+++ b/doc/source/acb_elliptic.rst
@@ -33,8 +33,10 @@ in [Car1995]_ and chapter 19 in [NIST2012]_.
\int_0^{\infty} \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}
where the square root extends continuously from positive infinity.
- The function is well-defined for `x,y,z \notin (-\infty,0)`, and with
+ The integral is well-defined for `x,y,z \notin (-\infty,0)`, and with
at most one of `x,y,z` being zero.
+ When some parameters are negative real numbers, the function is
+ still defined by analytic continuation.
In general, one or more duplication steps are applied until
`x,y,z` are close enough to use a multivariate Taylor polynomial
@@ -74,13 +76,26 @@ in [Car1995]_ and chapter 19 in [NIST2012]_.
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.
+ may be computed by setting *z* and *p* 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_rg(acb_t res, const acb_t x, const acb_t y, const acb_t z, int flags, slong prec)
+
+ Evaluates the Carlson symmetric elliptic integral of the second kind
+
+ .. math ::
+
+ R_G(x,y,z) = \frac{1}{4} \int_0^{\infty}
+ \frac{t}{\sqrt{(t+x)(t+y)(t+z)}}
+ \left( \frac{x}{t+x} + \frac{y}{t+y} + \frac{z}{t+z}\right) dt.
+
+ The evaluation is done by expressing `R_G` in terms of `R_F` and `R_D`.
+ There are no restrictions on the variables.
+
.. function:: void acb_elliptic_rc1(acb_t res, const acb_t x, slong prec)
This helper function computes the special case