add mag_rsqrt and improve arb_rsqrt for wide input

arb/rsqrt.c

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2012-2014 Fredrik Johansson
+    Copyright (C) 2012-2018 Fredrik Johansson
 
     This file is part of Arb.
 
@@ -9,56 +9,23 @@
     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
 */
 
-#include <math.h>
 #include "arb.h"
 
-void
-mag_pow_minus_three_half(mag_t z, const mag_t x)
-{
-    double t;
-
-    if (mag_is_zero(x))
-    {
-        mag_inf(z);
-    }
-    else if (mag_is_inf(x))
-    {
-        mag_zero(z);
-    }
-    else
-    {
-        if (fmpz_is_even(MAG_EXPREF(x)))
-        {
-            fmpz_mul_si(MAG_EXPREF(z), MAG_EXPREF(x), -3);
-            t = MAG_MAN(x) * (1.0 / (LIMB_ONE << MAG_BITS));
-        }
-        else
-        {
-            fmpz_add_ui(MAG_EXPREF(z), MAG_EXPREF(x), 1);
-            fmpz_mul_si(MAG_EXPREF(z), MAG_EXPREF(z), -3);
-            t = MAG_MAN(x) * (0.5 / (LIMB_ONE << MAG_BITS));
-        }
-
-        fmpz_tdiv_q_2exp(MAG_EXPREF(z), MAG_EXPREF(z), 1);
-        t = 1.0 / (t * sqrt(t));
-        t = t * (1.0 + 1e-10);
-
-        mag_set_d_2exp_fmpz(z, t, MAG_EXPREF(z));
-    }
-}
-
 void
 arb_rsqrt(arb_t z, const arb_t x, slong prec)
 {
+    mag_t t, u;
+    slong acc;
     int inexact;
 
-    if (arb_contains_nonpositive(x))
+    if (!arb_is_finite(x) || arf_sgn(arb_midref(x)) <= 0)
     {
-        arb_indeterminate(z);
-        return;
+        if (arf_is_pos_inf(arb_midref(x)) && mag_is_finite(arb_radref(x)))
+            arb_zero(z);
+        else
+            arb_indeterminate(z);
     }
-
-    if (arb_is_exact(x))
+    else if (mag_is_zero(arb_radref(x)))
     {
         inexact = arf_rsqrt(arb_midref(z), arb_midref(x), prec, ARB_RND);
 
@@ -69,21 +36,53 @@ arb_rsqrt(arb_t z, const arb_t x, slong prec)
     }
     else
     {
-        mag_t t;
         mag_init(t);
 
-        /* error bound: (1/2) (x-r)^(-3/2) * r */
         arb_get_mag_lower(t, x);
-        mag_pow_minus_three_half(t, t);
-        mag_mul(t, t, arb_radref(x));
-        mag_mul_2exp_si(t, t, -1);
 
-        inexact = arf_rsqrt(arb_midref(z), arb_midref(x), prec, ARB_RND);
+        acc = arb_rel_accuracy_bits(x);
+        acc = FLINT_MIN(acc, prec);
+        prec = FLINT_MIN(prec, acc + MAG_BITS);
+        prec = FLINT_MAX(prec, 2);
 
-        if (inexact)
-            arf_mag_add_ulp(arb_radref(z), t, arb_midref(z), prec);
+        if (acc <= 20)
+        {
+            if (mag_is_zero(t))
+            {
+                arb_indeterminate(z);
+            }
+            else
+            {
+                mag_init(u);
+                arb_get_mag(u, x);
+
+                mag_rsqrt(t, t);
+                mag_rsqrt_lower(u, u);
+
+                arb_set_interval_mag(z, u, t, prec);
+
+                mag_clear(u);
+            }
+        }
         else
-            mag_swap(arb_radref(z), t);
+        {
+            /* error bound: (1/2) (x-r)^(-3/2) * r */
+            mag_init(u);
+
+            mag_rsqrt(u, t);
+            mag_div(t, u, t);
+            mag_mul(t, t, arb_radref(x));
+            mag_mul_2exp_si(t, t, -1);
+
+            inexact = arf_rsqrt(arb_midref(z), arb_midref(x), prec, ARB_RND);
+
+            if (inexact)
+                arf_mag_add_ulp(arb_radref(z), t, arb_midref(z), prec);
+            else
+                mag_swap(arb_radref(z), t);
+
+            mag_clear(u);
+        }
 
         mag_clear(t);
     }

doc/source/mag.rst

@@ -327,13 +327,17 @@ Powers and logarithms
 
     Sets *z* to an upper bound for `\sqrt{x}`.
 
+.. function:: void mag_sqrt_lower(mag_t z, const mag_t x)
+
+    Sets *z* to a lower bound for `\sqrt{x}`.
+
 .. function:: void mag_rsqrt(mag_t z, const mag_t x)
 
     Sets *z* to an upper bound for `1/\sqrt{x}`.
 
-.. function:: void mag_sqrt_lower(mag_t z, const mag_t x)
+.. function:: void mag_rsqrt_lower(mag_t z, const mag_t x)
 
-    Sets *z* to a lower bound for `\sqrt{x}`.
+    Sets *z* to an lower bound for `1/\sqrt{x}`.
 
 .. function:: void mag_hypot(mag_t z, const mag_t x, const mag_t y)
 

mag.h

@@ -648,6 +648,7 @@ void mag_set_fmpz_2exp_fmpz_lower(mag_t z, const fmpz_t man, const fmpz_t exp);
 void mag_sqrt(mag_t y, const mag_t x);
 void mag_sqrt_lower(mag_t y, const mag_t x);
 void mag_rsqrt(mag_t y, const mag_t x);
+void mag_rsqrt_lower(mag_t y, const mag_t x);
 
 void mag_root(mag_t y, const mag_t x, ulong n);
 

mag/rsqrt.c

@@ -27,10 +27,11 @@ mag_rsqrt(mag_t y, const mag_t x)
         fmpz e;
 
         t = MAG_MAN(x) * ldexp(1.0, -MAG_BITS);
-        e = MAG_EXP(x);
 
-        if (!COEFF_IS_MPZ(e))
+        if (MAG_IS_LAGOM(x))
         {
+            e = MAG_EXP(x);
+
             if (e % 2 != 0)
             {
                 e = (1 - e) >> 1;
@@ -42,7 +43,9 @@ mag_rsqrt(mag_t y, const mag_t x)
             }
 
             t = (1.0 / sqrt(t)) * (1 + 1e-13);
-            mag_set_d_2exp_fmpz(y, t, &e);
+            _fmpz_demote(MAG_EXPREF(y));
+            MAG_SET_D_2EXP(MAG_MAN(y), MAG_EXP(y), t, e);
+
         }
         else
         {
@@ -55,3 +58,52 @@ mag_rsqrt(mag_t y, const mag_t x)
         }
     }
 }
+
+void
+mag_rsqrt_lower(mag_t y, const mag_t x)
+{
+    if (mag_is_special(x))
+    {
+        if (mag_is_zero(x))
+            mag_inf(y);
+        else
+            mag_zero(y);
+    }
+    else
+    {
+        double t;
+        fmpz e;
+
+        t = MAG_MAN(x) * ldexp(1.0, -MAG_BITS);
+
+        if (MAG_IS_LAGOM(x))
+        {
+            e = MAG_EXP(x);
+
+            if (e % 2 != 0)
+            {
+                e = (1 - e) >> 1;
+                t *= 2.0;
+            }
+            else
+            {
+                e = (-e) >> 1;
+            }
+
+            t = (1.0 / sqrt(t)) * (1 - 1e-13);
+            _fmpz_demote(MAG_EXPREF(y));
+            MAG_SET_D_2EXP_LOWER(MAG_MAN(y), MAG_EXP(y), t, e);
+
+        }
+        else
+        {
+            if (fmpz_is_odd(MAG_EXPREF(x)))
+                t *= 2.0;
+            fmpz_fdiv_q_2exp(MAG_EXPREF(y), MAG_EXPREF(x), 1);
+            fmpz_neg(MAG_EXPREF(y), MAG_EXPREF(y));
+            t = (1.0 / sqrt(t)) * (1 - 1e-13);
+            mag_set_d_2exp_fmpz_lower(y, t, MAG_EXPREF(y));
+        }
+    }
+}
+

mag/test/t-rsqrt_lower.c

@@ -0,0 +1,83 @@
+/*
+    Copyright (C) 2014 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 "mag.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("rsqrt_lower....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++)
+    {
+        fmpr_t x, y, z, z2;
+        mag_t xb, yb;
+
+        fmpr_init(x);
+        fmpr_init(y);
+        fmpr_init(z);
+        fmpr_init(z2);
+
+        mag_init(xb);
+        mag_init(yb);
+
+        mag_randtest_special(xb, state, 1 + n_randint(state, 200));
+        mag_randtest_special(yb, state, 1 + n_randint(state, 200));
+
+        mag_rsqrt_lower(yb, xb);
+
+        mag_get_fmpr(x, xb);
+        mag_get_fmpr(y, yb);
+
+        fmpr_rsqrt(z, x, MAG_BITS, FMPR_RND_DOWN);
+        fmpr_mul_ui(z2, z, 1023, MAG_BITS, FMPR_RND_DOWN);
+        fmpr_mul_2exp_si(z2, z2, -10);
+
+        MAG_CHECK_BITS(xb)
+        MAG_CHECK_BITS(yb)
+
+        if (!(fmpr_cmpabs(z2, y) <= 0 && fmpr_cmpabs(y, z) <= 0))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("x = "); fmpr_print(x); flint_printf("\n\n");
+            flint_printf("y = "); fmpr_print(y); flint_printf("\n\n");
+            flint_printf("z = "); fmpr_print(z); flint_printf("\n\n");
+            flint_printf("z2 = "); fmpr_print(z2); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        mag_rsqrt_lower(xb, xb);
+
+        if (!mag_equal(xb, yb))
+        {
+            flint_printf("FAIL (aliasing)\n\n");
+            flint_abort();
+        }
+
+        fmpr_clear(x);
+        fmpr_clear(y);
+        fmpr_clear(z);
+        fmpr_clear(z2);
+
+        mag_clear(xb);
+        mag_clear(yb);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}