diff --git a/doc/source/mag.rst b/doc/source/mag.rst index ede01d3b..95421082 100644 --- a/doc/source/mag.rst +++ b/doc/source/mag.rst @@ -285,6 +285,12 @@ Powers and logarithms Sets *z* to an upper bound for `\sqrt{x^2 + y^2}`. +.. function:: void mag_root(mag_t z, const mag_t x, ulong n) + + Sets *z* to an upper bound for `x^{1/n}`. + We evaluate `\exp(\log(1+2^{kn}x)/n) 2^{-n}`, where *k* is chosen + so that `2^{kn}x \approx 2^{30}`. + .. function:: void mag_log1p(mag_t z, const mag_t x) Sets *z* to an upper bound for `\log(1+x)`. The bound is computed diff --git a/mag.h b/mag.h index 098a8654..332c4bc3 100644 --- a/mag.h +++ b/mag.h @@ -581,6 +581,8 @@ 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_rsqrt(mag_t y, const mag_t x); +void mag_root(mag_t y, const mag_t x, ulong n); + void mag_hypot(mag_t z, const mag_t x, const mag_t y); void mag_binpow_uiui(mag_t b, ulong m, ulong n); diff --git a/mag/root.c b/mag/root.c new file mode 100644 index 00000000..7a085996 --- /dev/null +++ b/mag/root.c @@ -0,0 +1,69 @@ +/*============================================================================= + + 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 "mag.h" + +void +mag_root(mag_t y, const mag_t x, ulong n) +{ + if (n == 0) + { + mag_inf(y); + } + else if (n == 1 || mag_is_special(x)) + { + mag_set(y, x); + } + else if (n == 2) + { + mag_sqrt(y, x); + } + else if (n == 4) + { + mag_sqrt(y, x); + mag_sqrt(y, y); + } + else + { + fmpz_t e, f; + + fmpz_init_set_ui(e, MAG_BITS); + fmpz_init(f); + + fmpz_sub(e, e, MAG_EXPREF(x)); + fmpz_cdiv_q_ui(e, e, n); + fmpz_mul_ui(f, e, n); + mag_mul_2exp_fmpz(y, x, f); + mag_log1p(y, y); + mag_div_ui(y, y, n); + mag_exp(y, y); + fmpz_neg(e, e); + mag_mul_2exp_fmpz(y, y, e); + + fmpz_clear(e); + fmpz_clear(f); + } +} + diff --git a/mag/test/t-root.c b/mag/test/t-root.c new file mode 100644 index 00000000..7282c7e8 --- /dev/null +++ b/mag/test/t-root.c @@ -0,0 +1,99 @@ +/*============================================================================= + + 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("root...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, z, z2; + mag_t xb, yb; + ulong n; + + 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)); + n = n_randint(state, 30); + + mag_root(yb, xb, n); + + mag_get_fmpr(x, xb); + mag_get_fmpr(y, yb); + + fmpr_root(z, x, n, MAG_BITS, FMPR_RND_UP); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(yb) + + if (!(fmpr_cmpabs(z, y) <= 0 && fmpr_cmpabs(y, z2) <= 0)) + { + printf("FAIL\n\n"); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("y = "); fmpr_print(y); printf("\n\n"); + printf("z = "); fmpr_print(z); printf("\n\n"); + printf("z2 = "); fmpr_print(z2); printf("\n\n"); + abort(); + } + + mag_root(xb, xb, n); + + if (!mag_equal(xb, yb)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(z); + fmpr_clear(z2); + + mag_clear(xb); + mag_clear(yb); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}