add helper function for binary to decimal conversion

arb.h

@@ -456,6 +456,8 @@ int arb_contains(const arb_t x, const arb_t y);
 void arb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const arb_t x);
 int arb_get_unique_fmpz(fmpz_t z, const arb_t x);
 
+void arb_get_fmpz_mid_rad_10exp(fmpz_t mid, fmpz_t rad, fmpz_t exp, const arb_t x, long n);
+
 void arb_floor(arb_t z, const arb_t x, long prec);
 void arb_ceil(arb_t z, const arb_t x, long prec);
 

arb/get_fmpz_mid_rad_10exp.c

@@ -0,0 +1,102 @@
+/*=============================================================================
+
+    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 "arb.h"
+
+void
+arb_get_fmpz_mid_rad_10exp(fmpz_t mid, fmpz_t rad, fmpz_t exp, const arb_t x, long n)
+{
+    fmpz_t e, m;
+    arb_t t, u;
+    arf_t r;
+    long prec;
+    int roundmid, roundrad;
+
+    if (!arb_is_finite(x) || arb_is_zero(x))
+    {
+        fmpz_zero(mid);
+        fmpz_zero(rad);
+        fmpz_zero(exp);
+        return;
+    }
+
+    /*
+      We compute m such that x * 10^m ~= 10^(n+5).
+      If x = 2^e then m = (n+5) - e*log(2)/log(10).
+    */
+    fmpz_init(e);
+    fmpz_init(m);
+    arb_init(t);
+    arb_init(u);
+    arf_init(r);
+
+    if (arf_cmpabs_mag(arb_midref(x), arb_radref(x)) > 0)
+        fmpz_set(e, ARF_EXPREF(arb_midref(x)));
+    else
+        fmpz_set(e, ARF_EXPREF(arb_radref(x)));
+
+    prec = fmpz_bits(e) + 15;
+
+    arb_const_log2(t, prec);
+    arb_const_log10(u, prec);
+    arb_div(t, t, u, prec);
+    arb_mul_fmpz(t, t, e, prec);
+    arb_neg(t, t);
+    arb_add_ui(t, t, n + 5, prec);
+
+    arf_get_fmpz(m, arb_midref(t), ARF_RND_FLOOR);
+
+    arb_set_ui(t, 10);
+
+    fmpz_neg(exp, m);
+
+    prec = n * 3.32192809488736 + 30;
+
+    if (fmpz_sgn(m) >= 0)
+    {
+        arb_pow_fmpz_binexp(t, t, m, prec + 2 * fmpz_bits(m));
+        arb_mul(t, x, t, prec);
+    }
+    else
+    {
+        fmpz_neg(m, m);
+        arb_pow_fmpz_binexp(t, t, m, prec + 2 * fmpz_bits(m));
+        arb_div(t, x, t, prec);
+    }
+
+    roundmid = arf_get_fmpz_fixed_si(mid, arb_midref(t), 0);
+
+    arf_set_mag(r, arb_radref(t));
+    roundrad = arf_get_fmpz_fixed_si(rad, r, 0);
+
+    fmpz_add_ui(rad, rad, roundmid + roundrad);
+
+    fmpz_clear(e);
+    fmpz_clear(m);
+    arb_clear(t);
+    arb_clear(u);
+    arf_clear(r);
+}
+

arb/test/t-get_fmpz_mid_rad_10exp.c

@@ -0,0 +1,97 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("get_fmpz_mid_rad_10exp....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 5000; iter++)
+    {
+        arb_t x, y, t;
+        fmpz_t mid, rad, exp;
+        long n, prec;
+
+        arb_init(x);
+        arb_init(y);
+        arb_init(t);
+        fmpz_init(mid);
+        fmpz_init(rad);
+        fmpz_init(exp);
+
+        arb_randtest_special(x, state, 1 + n_randint(state, 500), 1 + n_randint(state, 500));
+        n = 1 + n_randint(state, 500);
+        prec = 2 + n_randint(state, 1500);
+
+        arb_get_fmpz_mid_rad_10exp(mid, rad, exp, x, n);
+
+        arf_set_fmpz(arb_midref(y), mid);
+        mag_set_fmpz(arb_radref(y), rad);
+        arb_set_ui(t, 10);
+        arb_pow_fmpz(t, t, exp, prec);
+        arb_mul(y, y, t, prec);
+
+        if (arb_is_finite(x) && !arb_is_zero(x) &&
+            fmpz_sizeinbase(mid, 10) < n && fmpz_sizeinbase(rad, 10) < n)
+        {
+            printf("FAIL (too few digits):\n\n");
+            printf("x = "); arb_printd(x, 50); printf("\n\n");
+            printf("y = "); arb_printd(y, 50); printf("\n\n");
+            printf("mid = "); fmpz_print(mid); printf("\n\n");
+            printf("rad = "); fmpz_print(rad); printf("\n\n");
+            printf("exp = "); fmpz_print(exp); printf("\n\n");
+            abort();
+        }
+
+        if (arb_is_finite(x) && !arb_contains(y, x))
+        {
+            printf("FAIL (containment):\n\n");
+            printf("x = "); arb_printd(x, 50); printf("\n\n");
+            printf("y = "); arb_printd(y, 50); printf("\n\n");
+            printf("mid = "); fmpz_print(mid); printf("\n\n");
+            printf("rad = "); fmpz_print(rad); printf("\n\n");
+            printf("exp = "); fmpz_print(exp); printf("\n\n");
+            abort();
+        }
+
+        arb_clear(x);
+        arb_clear(y);
+        arb_clear(t);
+        fmpz_clear(mid);
+        fmpz_clear(rad);
+        fmpz_clear(exp);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb.rst

@@ -296,6 +296,14 @@ Radius and interval operations
     Sets *y* to a ball containing `\lfloor x \rfloor` and `\lceil x \rceil`
     respectively, with the midpoint of *y* rounded to at most *prec* bits.
 
+.. function:: void arb_get_fmpz_mid_rad_10exp(fmpz_t mid, fmpz_t rad, fmpz_t exp, const arb_t x, long n)
+
+    Assuming that *x* is finite and not exactly zero, computes integers *mid*,
+    *rad*, *exp* such that `x \in [m-r, m+r] \times 10^e` and such that the
+    larger out of *mid* and *rad* has at least *n* digits plus a few guard
+    digits. If *x* is infinite or exactly zero, the outputs are all set
+    to zero.
+
 Comparisons
 -------------------------------------------------------------------------------