hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

add arf_root and improve arb_root

Browse source
This commit is contained in:
Fredrik Johansson 2015-06-17 18:37:53 +02:00
parent e9f7d19561
commit 303913f631
9 changed files with 409 additions and 26 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
arb/fmprb_wrappers.c
View file

@ -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);
}

31
arb/pow_fmpq.c
View file

@ -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);
}
}
}

125
arb/root.c Normal file
View file

@ -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);
}
}

7
arb/test/t-root.c
View file

@ -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);

2
arf.h
View file

@ -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);

105
arf/root.c Normal file
View file

@ -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;
}

129
arf/test/t-root.c Normal file
View file

@ -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;
}

19
doc/source/arb.rst
View file

@ -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)

7
doc/source/arf.rst
View file

@ -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
-------------------------------------------------------------------------------