diff --git a/arb/rsqrt.c b/arb/rsqrt.c
index 50cd2bc2..e206dffd 100644
--- a/arb/rsqrt.c
+++ b/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 .
*/
-#include
#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);
}
diff --git a/doc/source/mag.rst b/doc/source/mag.rst
index 8656a062..9bfec1d2 100644
--- a/doc/source/mag.rst
+++ b/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)
diff --git a/mag.h b/mag.h
index 1c582e2b..965e8625 100644
--- a/mag.h
+++ b/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);
diff --git a/mag/rsqrt.c b/mag/rsqrt.c
index 623ebbe8..2814d9df 100644
--- a/mag/rsqrt.c
+++ b/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));
+ }
+ }
+}
+
diff --git a/mag/test/t-rsqrt_lower.c b/mag/test/t-rsqrt_lower.c
new file mode 100644
index 00000000..b7bbf352
--- /dev/null
+++ b/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 .
+*/
+
+#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;
+}