add acb_rel_accuracy_bits and improve the real version

acb.h

@@ -834,6 +834,15 @@ void acb_randtest_special(acb_t z, flint_rand_t state, long prec, long mag_bits)
 
 void acb_randtest_precise(acb_t z, flint_rand_t state, long prec, long mag_bits);
 
+long acb_rel_error_bits(const acb_t x);
+
+ARB_INLINE long
+acb_rel_accuracy_bits(const acb_t x)
+{
+    return -acb_rel_error_bits(x);
+}
+
+
 ACB_INLINE long
 acb_bits(const acb_t x)
 {

acb/rel_error_bits.c

@@ -0,0 +1,80 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb.h"
+
+long
+acb_rel_error_bits(const acb_t x)
+{
+    int am, ar, bm, br;
+    long result;
+    const fmpz * radmag;
+    const fmpz * midmag;
+    fmpz_t t;
+
+    if (!acb_is_finite(x))
+        return ARF_PREC_EXACT;
+
+    am = !arf_is_zero(arb_midref(acb_realref(x)));
+    ar = !mag_is_zero(arb_radref(acb_realref(x)));
+
+    bm = !arf_is_zero(arb_midref(acb_imagref(x)));
+    br = !mag_is_zero(arb_radref(acb_imagref(x)));
+
+    /* no radius -- exact */
+    if (!ar && !br)
+        return -ARF_PREC_EXACT;
+
+    /* no midpoint -- infinite relative error */
+    if (!am && !bm)
+        return ARF_PREC_EXACT;
+
+#define ame ARF_EXPREF(arb_midref(acb_realref(x)))
+#define are MAG_EXPREF(arb_radref(acb_realref(x)))
+#define bme ARF_EXPREF(arb_midref(acb_imagref(x)))
+#define bre MAG_EXPREF(arb_radref(acb_imagref(x)))
+
+    if (am && bm)
+        midmag = fmpz_cmp(ame, bme) >= 0 ? ame : bme;
+    else if (am)
+        midmag = ame;
+    else
+        midmag = bme;
+
+    if (ar && br)
+        radmag = fmpz_cmp(are, bre) >= 0 ? are : bre;
+    else if (ar)
+        radmag = are;
+    else
+        radmag = bre;
+
+    fmpz_init(t);
+    fmpz_add_ui(t, radmag, 1);
+    result = _fmpz_sub_small(t, midmag);
+    fmpz_clear(t);
+
+    return result;
+}
+

acb/test/t-rel_accuracy_bits.c

@@ -0,0 +1,116 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2015 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("rel_accuracy_bits....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* test aliasing of c and a */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arb_t x;
+        acb_t z;
+        long a1, a2;
+
+        arb_init(x);
+        acb_init(z);
+
+        arb_randtest_special(x, state, 1 + n_randint(state, 200), 1 + n_randint(state, 200));
+        acb_set_arb(z, x);
+
+        a1 = arb_rel_accuracy_bits(x);
+        a2 = acb_rel_accuracy_bits(z);
+
+        if (a1 != a2)
+        {
+            printf("FAIL: acb != arb\n\n");
+            printf("x = "); arb_print(x); printf("\n\n");
+            printf("z = "); acb_print(z); printf("\n\n");
+            printf("a1 = %ld, a2 = %ld\n\n", a1, a2);
+            abort();
+        }
+
+        acb_randtest_special(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 200));
+
+        a1 = acb_rel_accuracy_bits(z);
+
+        if (n_randint(state, 2))
+            arf_swap(arb_midref(acb_realref(z)), arb_midref(acb_imagref(z)));
+
+        if (n_randint(state, 2))
+            mag_swap(arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
+
+        a2 = acb_rel_accuracy_bits(z);
+
+        if (a1 != a2)
+        {
+            printf("FAIL: swapping\n\n");
+            printf("z = "); acb_print(z); printf("\n\n");
+            printf("a1 = %ld, a2 = %ld\n\n", a1, a2);
+            abort();
+        }
+
+        acb_randtest(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 200));
+
+        if (arf_cmpabs(arb_midref(acb_realref(z)), arb_midref(acb_imagref(z))) >= 0)
+            arf_set(arb_midref(x), arb_midref(acb_realref(z)));
+        else
+            arf_set(arb_midref(x), arb_midref(acb_imagref(z)));
+
+        if (mag_cmp(arb_radref(acb_realref(z)), arb_radref(acb_imagref(z))) >= 0)
+            mag_set(arb_radref(x), arb_radref(acb_realref(z)));
+        else
+            mag_set(arb_radref(x), arb_radref(acb_imagref(z)));
+
+        a1 = acb_rel_accuracy_bits(z);
+        a2 = arb_rel_accuracy_bits(x);
+
+        if (a1 != a2)
+        {
+            printf("FAIL: acb != arb (2)\n\n");
+            printf("x = "); arb_print(x); printf("\n\n");
+            printf("z = "); acb_print(z); printf("\n\n");
+            printf("a1 = %ld, a2 = %ld\n\n", a1, a2);
+            abort();
+        }
+
+        arb_clear(x);
+        acb_clear(z);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arb/rel_error_bits.c

@@ -28,28 +28,24 @@
 long
 arb_rel_error_bits(const arb_t x)
 {
-    fmpz_t midmag, radmag;
-    arf_t t;
+    fmpz_t t;
     long result;
 
     if (mag_is_zero(arb_radref(x)))
-        return -ARF_PREC_EXACT;
+    {
+        if (arf_is_finite(arb_midref(x)))
+            return -ARF_PREC_EXACT;
+        else
+            return ARF_PREC_EXACT;
+    }
+
     if (arf_is_special(arb_midref(x)) || mag_is_inf(arb_radref(x)))
         return ARF_PREC_EXACT;
 
-    fmpz_init(midmag);
-    fmpz_init(radmag);
-
-    arf_abs_bound_lt_2exp_fmpz(midmag, arb_midref(x));
-
-    arf_init_set_mag_shallow(t, arb_radref(x)); /* no need to free */
-    arf_abs_bound_lt_2exp_fmpz(radmag, t);
-    fmpz_add_ui(radmag, radmag, 1);
-
-    result = _fmpz_sub_small(radmag, midmag);
-
-    fmpz_clear(midmag);
-    fmpz_clear(radmag);
+    fmpz_init(t);
+    fmpz_add_ui(t, MAG_EXPREF(arb_radref(x)), 1);
+    result = _fmpz_sub_small(t, ARF_EXPREF(arb_midref(x)));
+    fmpz_clear(t);
 
     return result;
 }

doc/source/acb.rst

@@ -210,6 +210,19 @@ Precision and comparisons
 
     Returns nonzero iff zero is contained in *x*.
 
+.. function:: long acb_rel_error_bits(const acb_t x)
+
+    Returns the effective relative error of *x* measured in bits.
+    This is computed as if calling :func:`arb_rel_error_bits` on the
+    real ball whose midpoint is the larger out of the real and imaginary
+    midpoints of *x*, and whose radius is the larger out of the real
+    and imaginary radiuses of *x*.
+
+.. function:: long acb_rel_accuracy_bits(const arb_t x)
+
+    Returns the effective relative accuracy of *x* measured in bits,
+    equal to the negative of the return value from :func:`acb_rel_error_bits`.
+
 .. function:: long acb_bits(const acb_t x)
 
     Returns the maximum of *arb_bits* applied to the real