From 303913f631c29864c2ff65891fb62a5d208c507d Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 17 Jun 2015 18:37:53 +0200 Subject: [PATCH] add arf_root and improve arb_root --- arb/fmprb_wrappers.c | 10 ---- arb/pow_fmpq.c | 31 +++++++---- arb/root.c | 125 +++++++++++++++++++++++++++++++++++++++++ arb/test/t-root.c | 7 ++- arf.h | 2 + arf/root.c | 105 +++++++++++++++++++++++++++++++++++ arf/test/t-root.c | 129 +++++++++++++++++++++++++++++++++++++++++++ doc/source/arb.rst | 19 ++++++- doc/source/arf.rst | 7 +++ 9 files changed, 409 insertions(+), 26 deletions(-) create mode 100644 arb/root.c create mode 100644 arf/root.c create mode 100644 arf/test/t-root.c diff --git a/arb/fmprb_wrappers.c b/arb/fmprb_wrappers.c index 1c191259..4645347e 100644 --- a/arb/fmprb_wrappers.c +++ b/arb/fmprb_wrappers.c @@ -35,13 +35,3 @@ void arb_div_2expm1_ui(arb_t z, const arb_t x, ulong n, long prec) fmprb_clear(t); } -void arb_root(arb_t z, const arb_t x, ulong k, long prec) -{ - fmprb_t t; - fmprb_init(t); - arb_get_fmprb(t, x); - fmprb_root(t, t, k, prec); - arb_set_fmprb(z, t); - fmprb_clear(t); -} - diff --git a/arb/pow_fmpq.c b/arb/pow_fmpq.c index 563e51a9..887b4062 100644 --- a/arb/pow_fmpq.c +++ b/arb/pow_fmpq.c @@ -32,21 +32,30 @@ arb_pow_fmpq(arb_t y, const arb_t x, const fmpq_t a, long prec) { arb_pow_fmpz(y, x, fmpq_numref(a), prec); } - else if (fmpz_cmp_ui(fmpq_denref(a), 36) <= 0) - { - arb_root(y, x, *fmpq_denref(a), prec); - arb_pow_fmpz(y, y, fmpq_numref(a), prec); - } else { - long wp; + int use_exp; + long k = *fmpq_denref(a); - wp = prec + 10; + if (k == 2 || k == 4) + use_exp = 0; + else if (k < 50) + use_exp = prec < (1L << ((k / 8) + 8)); + else + use_exp = 1; - arb_log(y, x, wp); - arb_mul_fmpz(y, y, fmpq_numref(a), wp); - arb_div_fmpz(y, y, fmpq_denref(a), wp); - arb_exp(y, y, prec); + if (use_exp) + { + arb_log(y, x, prec + 10); + arb_mul_fmpz(y, y, fmpq_numref(a), prec + 10); + arb_div_fmpz(y, y, fmpq_denref(a), prec + 10); + arb_exp(y, y, prec); + } + else + { + arb_root(y, x, *fmpq_denref(a), prec); + arb_pow_fmpz(y, y, fmpq_numref(a), prec); + } } } diff --git a/arb/root.c b/arb/root.c new file mode 100644 index 00000000..e9e96ecb --- /dev/null +++ b/arb/root.c @@ -0,0 +1,125 @@ +/*============================================================================= + + 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" + +static void +arb_root_arf(arb_t z, const arf_t x, ulong k, long prec) +{ + int inexact = arf_root(arb_midref(z), x, k, prec, ARB_RND); + + if (inexact) + arf_mag_set_ulp(arb_radref(z), arb_midref(z), prec); + else + mag_zero(arb_radref(z)); +} + +void +arb_root_algebraic(arb_t res, const arb_t x, ulong k, long prec) +{ + mag_t r, msubr, m1k, t; + + if (arb_is_exact(x)) + { + arb_root_arf(res, arb_midref(x), k, prec); + return; + } + + if (!arb_is_nonnegative(x)) + { + arb_indeterminate(res); + return; + } + + mag_init(r); + mag_init(msubr); + mag_init(m1k); + mag_init(t); + + /* x = [m-r, m+r] */ + mag_set(r, arb_radref(x)); + /* m - r */ + arb_get_mag_lower(msubr, x); + + /* m^(1/k) */ + arb_root_arf(res, arb_midref(x), k, prec); + + /* bound for m^(1/k) */ + arb_get_mag(m1k, res); + + /* C = min(1, log(1+r/(m-r))/k) */ + mag_div(t, r, msubr); + mag_log1p(t, t); + mag_div_ui(t, t, k); + if (mag_cmp_2exp_si(t, 0) > 0) + mag_one(t); + + /* C m^(1/k) */ + mag_mul(t, m1k, t); + mag_add(arb_radref(res), arb_radref(res), t); + + mag_clear(r); + mag_clear(msubr); + mag_clear(m1k); + mag_clear(t); +} + +void +arb_root_exp(arb_t res, const arb_t x, ulong k, long prec) +{ + arb_log(res, x, prec + 4); + arb_div_ui(res, res, k, prec + 4); + arb_exp(res, res, prec); +} + +void +arb_root(arb_t res, const arb_t x, ulong k, long prec) +{ + if (k == 0) + { + arb_indeterminate(res); + } + else if (k == 1) + { + arb_set_round(res, x, prec); + } + else if (k == 2) + { + arb_sqrt(res, x, prec); + } + else if (k == 4) + { + arb_sqrt(res, x, prec + 2); + arb_sqrt(res, res, prec); + } + else + { + if (k > 50 || prec < (1L << ((k / 8) + 8))) + arb_root_exp(res, x, k, prec); + else + arb_root_algebraic(res, x, k, prec); + } +} + diff --git a/arb/test/t-root.c b/arb/test/t-root.c index 97036188..9f179df2 100644 --- a/arb/test/t-root.c +++ b/arb/test/t-root.c @@ -41,14 +41,15 @@ int main() ulong k; long prec; - prec = 2 + n_randint(state, 200); - k = 1 + n_randint(state, 10); + prec = 2 + n_randint(state, 2000); + k = n_randtest_not_zero(state); arb_init(a); arb_init(b); arb_init(c); - arb_randtest(a, state, 1 + n_randint(state, 200), 10); + arb_randtest(a, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100)); + arb_randtest(b, state, 1 + n_randint(state, 2000), 1 + n_randint(state, 100)); arb_root(b, a, k, prec); arb_pow_ui(c, b, k, prec); diff --git a/arf.h b/arf.h index 1dcf4948..18ef1326 100644 --- a/arf.h +++ b/arf.h @@ -1190,6 +1190,8 @@ int arf_sqrt_fmpz(arf_t z, const fmpz_t x, long prec, arf_rnd_t rnd); int arf_rsqrt(arf_ptr z, arf_srcptr x, long prec, arf_rnd_t rnd); +int arf_root(arf_t z, const arf_t x, ulong k, long prec, arf_rnd_t rnd); + /* Magnitude bounds */ void arf_get_mag(mag_t y, const arf_t x); diff --git a/arf/root.c b/arf/root.c new file mode 100644 index 00000000..cf306d80 --- /dev/null +++ b/arf/root.c @@ -0,0 +1,105 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_root(arf_ptr z, arf_srcptr x, ulong k, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, zn, val; + mp_srcptr xptr; + mp_ptr tmp, zptr; + mpfr_t xf, zf; + fmpz_t q, r; + int inexact; + + if (k == 0) + { + arf_nan(z); + return 0; + } + + if (k == 1) + return arf_set_round(z, x, prec, rnd); + + if (k == 2) + return arf_sqrt(z, x, prec, rnd); + + if (arf_is_special(x)) + { + if (arf_is_neg_inf(x)) + arf_nan(z); + else + arf_set(z, x); + return 0; + } + + if (ARF_SGNBIT(x)) + { + arf_nan(z); + return 0; + } + + fmpz_init(q); + fmpz_init(r); + + /* x = m * 2^e where e = qk + r */ + /* x^(1/k) = (m * 2^(qk+r))^(1/k) */ + /* x^(1/k) = (m * 2^r)^(1/k) * 2^q */ + fmpz_set_ui(r, k); + fmpz_fdiv_qr(q, r, ARF_EXPREF(x), r); + + ARF_GET_MPN_READONLY(xptr, xn, x); + zn = (prec + FLINT_BITS - 1) / FLINT_BITS; + + zf->_mpfr_d = tmp = flint_malloc(zn * sizeof(mp_limb_t)); + zf->_mpfr_prec = prec; + zf->_mpfr_sign = 1; + zf->_mpfr_exp = 0; + + xf->_mpfr_d = (mp_ptr) xptr; + xf->_mpfr_prec = xn * FLINT_BITS; + xf->_mpfr_sign = 1; + xf->_mpfr_exp = fmpz_get_ui(r); + + inexact = mpfr_root(zf, xf, k, arf_rnd_to_mpfr(rnd)); + inexact = (inexact != 0); + + val = 0; + while (tmp[val] == 0) + val++; + + ARF_GET_MPN_WRITE(zptr, zn - val, z); + flint_mpn_copyi(zptr, tmp + val, zn - val); + + fmpz_add_si(ARF_EXPREF(z), q, zf->_mpfr_exp); + + flint_free(tmp); + fmpz_clear(q); + fmpz_clear(r); + + return inexact; +} + diff --git a/arf/test/t-root.c b/arf/test/t-root.c new file mode 100644 index 00000000..c08555cf --- /dev/null +++ b/arf/test/t-root.c @@ -0,0 +1,129 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_root_naive(arf_t z, const arf_t x, ulong k, long prec, arf_rnd_t rnd) +{ + fmpr_t a; + long r; + + fmpr_init(a); + + arf_get_fmpr(a, x); + + r = fmpr_root(a, a, k, prec, rnd); + + arf_set_fmpr(z, a); + + fmpr_clear(a); + + return (r == FMPR_RESULT_EXACT) ? 0 : 1; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("root...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + long prec, r1, r2; + ulong k; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 10; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + prec = 2 + n_randint(state, 2000); + k = n_randint(state, 50); + + if (n_randint(state, 20) == 0) + arf_mul(x, x, x, prec, ARF_RND_DOWN); + else if (n_randint(state, 20) == 0) + arf_mul(x, x, x, prec, ARF_RND_UP); + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_root(z, x, k, prec, rnd); + r2 = arf_root_naive(v, x, k, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("k = %lu, prec = %ld, rnd = %d\n\n", k, prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_root_naive(v, x, k, prec, rnd); + r1 = arf_root(x, x, k, prec, rnd); + if (!arf_equal(v, x) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("k = %lu, prec = %ld, rnd = %d\n\n", k, prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/doc/source/arb.rst b/doc/source/arb.rst index 32f629eb..2e95393c 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -680,8 +680,23 @@ Powers and roots .. function:: void arb_root(arb_t z, const arb_t x, ulong k, long prec) Sets *z* to the *k*-th root of *x*, rounded to *prec* bits. - As currently implemented, this function is only fast for small *k*. - For large *k* it is better to use :func:`arb_pow_fmpq` or :func:`arb_pow`. + This function selects between different algorithms. For large *k*, + it evaluates `\exp(\log(x)/k)`. For small *k*, it uses :func:`arf_root` + at the midpoint and computes a propagated error bound as follows: + if input interval is `[m-r, m+r]` with `r \le m`, the error is largest at + `m-r` where it satisfies + + .. math :: + + m^{1/k} - (m-r)^{1/k} = m^{1/k} [1 - (1-r/m)^{1/k}] + + = m^{1/k} [1 - \exp(\log(1-r/m)/k)] + + \le m^{1/k} \min(1, -\log(1-r/m)/k) + + = m^{1/k} \min(1, \log(1+r/(m-r))/k). + + This is evaluated using :func:`mag_log1p`. .. function:: void arb_pow_fmpz_binexp(arb_t y, const arb_t b, const fmpz_t e, long prec) diff --git a/doc/source/arf.rst b/doc/source/arf.rst index 0c17ac88..a4f09775 100644 --- a/doc/source/arf.rst +++ b/doc/source/arf.rst @@ -595,6 +595,13 @@ Square roots returning nonzero iff the operation is inexact. The result is NaN if *x* is negative, and `+\infty` if *x* is zero. +.. function:: int arf_root(arf_t z, const arf_t x, ulong k, long prec, arf_rnd_t rnd) + + Sets `z = x^{1/k}`, rounded to *prec* bits in the direction specified by *rnd*, + returning nonzero iff the operation is inexact. The result is NaN if *x* is negative. + Warning: this function is a wrapper around the MPFR root function. + It gets slow and uses much memory for large *k*. + Complex arithmetic -------------------------------------------------------------------------------