Carlson elliptic integral of the second kind

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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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);
+    }
+}
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

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