diff --git a/doc/source/mag.rst b/doc/source/mag.rst index bf0c5757..fcd3e24e 100644 --- a/doc/source/mag.rst +++ b/doc/source/mag.rst @@ -164,6 +164,14 @@ Arithmetic Sets *z* to a lower bound for `x^e`. +.. function:: void mag_sqrt(mag_t z, const mag_t x) + + Sets *z* to an upper 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}`. + Fast, unsafe versions ------------------------------------------------------------------------------- diff --git a/mag.h b/mag.h index 68f053c7..e3d3dcd5 100644 --- a/mag.h +++ b/mag.h @@ -649,7 +649,6 @@ double mag_d_log_lower_bound(double x); void mag_log1p(mag_t z, const mag_t x); -/* TODO: test */ void mag_pow_ui(mag_t z, const mag_t x, ulong e); void mag_pow_ui_lower(mag_t z, const mag_t x, ulong e); @@ -663,7 +662,6 @@ void mag_bernoulli_div_fac_ui(mag_t z, ulong n); /* TODO: test */ void mag_set_fmpz_2exp_fmpz_lower(mag_t z, const fmpz_t man, const fmpz_t exp); -/* TODO: document */ void mag_sqrt(mag_t y, const mag_t x); void mag_rsqrt(mag_t y, const mag_t x); diff --git a/mag/pow_ui.c b/mag/pow_ui.c index 96095986..974b8939 100644 --- a/mag/pow_ui.c +++ b/mag/pow_ui.c @@ -28,7 +28,7 @@ void mag_pow_ui(mag_t z, const mag_t x, ulong e) { - if (mag_is_inf(z)) + if (mag_is_inf(x)) { mag_inf(z); } @@ -74,7 +74,7 @@ mag_pow_ui_lower(mag_t z, const mag_t x, ulong e) else mag_mul_lower(z, x, x); } - else if (mag_is_inf(z)) + else if (mag_is_inf(x)) { mag_inf(z); } diff --git a/mag/test/t-pow_ui.c b/mag/test/t-pow_ui.c new file mode 100644 index 00000000..cef36b0e --- /dev/null +++ b/mag/test/t-pow_ui.c @@ -0,0 +1,101 @@ +/*============================================================================= + + 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) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("pow_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, w; + mag_t xb, yb; + ulong e; + + fmpr_init(x); + fmpr_init(y); + fmpr_init(w); + + mag_init(xb); + mag_init(yb); + + mag_randtest_special(xb, state, 80); + mag_randtest_special(yb, state, 80); + e = n_randtest(state); + + mag_get_fmpr(x, xb); + + fmpr_pow_sloppy_ui(y, x, e, 2 * MAG_BITS, FMPR_RND_UP); + if (fmpr_is_nan(y)) + fmpr_pos_inf(y); + + mag_pow_ui(yb, xb, e); + mag_get_fmpr(w, yb); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(yb) + + if (!(fmpr_cmpabs(y, w) <= 0)) + { + printf("FAIL\n\n"); + printf("e = %lu\n\n", e); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("y = "); fmpr_print(y); printf("\n\n"); + printf("w = "); fmpr_print(w); printf("\n\n"); + abort(); + } + + mag_pow_ui(xb, xb, e); + + if (!mag_equal(xb, yb)) + { + printf("FAIL (aliasing)\n\n"); + printf("e = %lu\n\n", e); + mag_print(xb); printf("\n\n"); + mag_print(yb); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(w); + + mag_clear(xb); + mag_clear(yb); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/mag/test/t-pow_ui_lower.c b/mag/test/t-pow_ui_lower.c new file mode 100644 index 00000000..1aa75914 --- /dev/null +++ b/mag/test/t-pow_ui_lower.c @@ -0,0 +1,101 @@ +/*============================================================================= + + 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) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("pow_ui_lower...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, w; + mag_t xb, yb; + ulong e; + + fmpr_init(x); + fmpr_init(y); + fmpr_init(w); + + mag_init(xb); + mag_init(yb); + + mag_randtest_special(xb, state, 80); + mag_randtest_special(yb, state, 80); + e = n_randtest(state); + + mag_get_fmpr(x, xb); + + fmpr_pow_sloppy_ui(y, x, e, 2 * MAG_BITS, FMPR_RND_DOWN); + if (fmpr_is_nan(y)) + fmpr_pos_inf(y); + + mag_pow_ui_lower(yb, xb, e); + mag_get_fmpr(w, yb); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(yb) + + if (!(fmpr_cmpabs(w, y) <= 0)) + { + printf("FAIL\n\n"); + printf("e = %lu\n\n", e); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("y = "); fmpr_print(y); printf("\n\n"); + printf("w = "); fmpr_print(w); printf("\n\n"); + abort(); + } + + mag_pow_ui_lower(xb, xb, e); + + if (!mag_equal(xb, yb)) + { + printf("FAIL (aliasing)\n\n"); + printf("e = %lu\n\n", e); + mag_print(xb); printf("\n\n"); + mag_print(yb); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(w); + + mag_clear(xb); + mag_clear(yb); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +