diff options
| -rw-r--r-- | NEWS | 4 | ||||
| -rw-r--r-- | doc/mpfr.texi | 11 | ||||
| -rw-r--r-- | src/Makefile.am | 3 | ||||
| -rw-r--r-- | src/mpfr-impl.h | 3 | ||||
| -rw-r--r-- | src/mpfr.h | 10 | ||||
| -rw-r--r-- | src/nbits_ulong.c | 54 | ||||
| -rw-r--r-- | src/pow_si.c | 40 | ||||
| -rw-r--r-- | src/pow_ui.c | 22 | ||||
| -rw-r--r-- | tests/tpow.c | 3 |
9 files changed, 118 insertions, 32 deletions
@@ -26,9 +26,9 @@ Changes from versions 4.1.* to version 4.2.0: mpfr_atanu and mpfr_atan2u. - New functions mpfr_cospi, mpfr_sinpi, mpfr_tanpi, mpfr_acospi, mpfr_asinpi, mpfr_atanpi and mpfr_atan2pi. -- New function mpfr_log2p1, mpfr_log10p1, mpfr_exp2m1, mpfr_exp10m1 and +- New functions mpfr_log2p1, mpfr_log10p1, mpfr_exp2m1, mpfr_exp10m1 and mpfr_compound. -- New function mpfr_fmod_ui. +- New functions mpfr_fmod_ui, mpfr_pow_uj and mpfr_pow_sj. - Bug fixes. In particular, for the formatted output functions (mpfr_printf, etc.), the case where the precision consists only of a period has been fixed diff --git a/doc/mpfr.texi b/doc/mpfr.texi index 023afc7a5..87a031c82 100644 --- a/doc/mpfr.texi +++ b/doc/mpfr.texi @@ -2242,15 +2242,16 @@ rounded in the direction @var{rnd}. @deftypefun int mpfr_pow (mpfr_t @var{rop}, mpfr_t @var{op1}, mpfr_t @var{op2}, mpfr_rnd_t @var{rnd}) @deftypefunx int mpfr_pow_ui (mpfr_t @var{rop}, mpfr_t @var{op1}, unsigned long int @var{op2}, mpfr_rnd_t @var{rnd}) @deftypefunx int mpfr_pow_si (mpfr_t @var{rop}, mpfr_t @var{op1}, long int @var{op2}, mpfr_rnd_t @var{rnd}) -@deftypefunx int mpfr_pow_n (mpfr_t @var{rop}, mpfr_t @var{op1}, long int @var{op2}, mpfr_rnd_t @var{rnd}) +@deftypefunx int mpfr_pow_uj (mpfr_t @var{rop}, mpfr_t @var{op1}, uintmax_t @var{op2}, mpfr_rnd_t @var{rnd}) +@deftypefunx int mpfr_pow_sj (mpfr_t @var{rop}, mpfr_t @var{op1}, intmax_t @var{op2}, mpfr_rnd_t @var{rnd}) +@deftypefunx int mpfr_pow_n (mpfr_t @var{rop}, mpfr_t @var{op1}, intmax_t @var{op2}, mpfr_rnd_t @var{rnd}) @deftypefunx int mpfr_pow_z (mpfr_t @var{rop}, mpfr_t @var{op1}, mpz_t @var{op2}, mpfr_rnd_t @var{rnd}) @deftypefunx int mpfr_ui_pow_ui (mpfr_t @var{rop}, unsigned long int @var{op1}, unsigned long int @var{op2}, mpfr_rnd_t @var{rnd}) @deftypefunx int mpfr_ui_pow (mpfr_t @var{rop}, unsigned long int @var{op1}, mpfr_t @var{op2}, mpfr_rnd_t @var{rnd}) Set @var{rop} to @m{@var{op1}^{@var{op2}}, @var{op1} raised to @var{op2}}, rounded in the direction @var{rnd}. -@c FIXME: mpfr_pown should be an alias of mpfr_pow_sj to follow ISO C2x. -@c The function @code{mpfr_pow_n} is just an alias for @code{mpfr_pow_si}, -@c %to mimic the IEEE@tie{}754-2019 function @code{pown}. +(The @code{mpfr_pow_n} function is just an alias for @code{mpfr_pow_sj}, +to follow the C2x function @code{pown}.) Special values are handled as described in the ISO C99 and IEEE@tie{}754-2008 standards for the @code{pow} function: @itemize @bullet @@ -4279,6 +4280,8 @@ use @code{mpfr_get_z_exp}. @item @code{mpfr_nrandom} in MPFR@tie{}4.0. +@item @code{mpfr_pow_sj} and @code{mpfr_pow_uj} in MPFR@tie{}4.2. + @item @code{mpfr_printf} in MPFR@tie{}2.4. @item @code{mpfr_rec_sqrt} in MPFR@tie{}2.4. diff --git a/src/Makefile.am b/src/Makefile.am index fea49b27a..935f54990 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -38,7 +38,8 @@ div_2ui.c div.c div_ui.c dot.c dump.c eq.c exp10.c exp2.c exp3.c exp.c \ frac.c frexp.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c \ isinf.c isnan.c isnum.c const_log2.c log.c modf.c mul_2exp.c mul_2si.c \ mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c printf.c vasprintf.c \ -const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c \ +const_pi.c pow.c pow_si.c pow_ui.c pow_sj.c pow_uj.c print_raw.c \ +print_rnd_mode.c \ reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c \ set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c \ set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c \ diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h index 5707eec10..b7e7d1792 100644 --- a/src/mpfr-impl.h +++ b/src/mpfr-impl.h @@ -2489,6 +2489,9 @@ __MPFR_DECLSPEC void mpfr_mpz_clear (mpz_ptr); __MPFR_DECLSPEC int mpfr_odd_p (mpfr_srcptr); __MPFR_DECLSPEC int mpfr_nbits_ulong (unsigned long); +#ifdef _MPFR_H_HAVE_INTMAX_T +__MPFR_DECLSPEC int mpfr_nbits_uj (uintmax_t); +#endif #ifdef _MPFR_H_HAVE_VA_LIST /* Declared only if <stdarg.h> has been included. */ diff --git a/src/mpfr.h b/src/mpfr.h index 220d6af05..5a4d3cfa1 100644 --- a/src/mpfr.h +++ b/src/mpfr.h @@ -563,10 +563,6 @@ __MPFR_DECLSPEC int mpfr_snprintf (char*, size_t, const char*, ...); __MPFR_DECLSPEC int mpfr_pow (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_pow_si (mpfr_ptr, mpfr_srcptr, long, mpfr_rnd_t); -/* FIXME: mpfr_pown should be an alias of mpfr_pow_sj to follow ISO C2x, - where pown is defined with intmax_t n. */ -/* define mpfr_pown (defined in IEEE 754-2019) as an alias for mpfr_pow_si */ -#define mpfr_pown mpfr_pow_si __MPFR_DECLSPEC int mpfr_compound (mpfr_ptr, mpfr_srcptr, long, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_pow_ui (mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t); @@ -1119,6 +1115,8 @@ extern "C" { #define mpfr_set_uj_2exp __gmpfr_set_uj_2exp #define mpfr_get_sj __gmpfr_mpfr_get_sj #define mpfr_get_uj __gmpfr_mpfr_get_uj +#define mpfr_pow_uj __gmpfr_mpfr_pow_uj +#define mpfr_pow_sj __gmpfr_mpfr_pow_sj __MPFR_DECLSPEC int mpfr_set_sj (mpfr_ptr, intmax_t, mpfr_rnd_t); __MPFR_DECLSPEC int mpfr_set_sj_2exp (mpfr_ptr, intmax_t, intmax_t, mpfr_rnd_t); @@ -1127,6 +1125,10 @@ __MPFR_DECLSPEC int mpfr_set_uj_2exp (mpfr_ptr, uintmax_t, intmax_t, mpfr_rnd_t); __MPFR_DECLSPEC intmax_t mpfr_get_sj (mpfr_srcptr, mpfr_rnd_t); __MPFR_DECLSPEC uintmax_t mpfr_get_uj (mpfr_srcptr, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_pow_uj (mpfr_ptr, mpfr_srcptr, uintmax_t, mpfr_rnd_t); +__MPFR_DECLSPEC int mpfr_pow_sj (mpfr_ptr, mpfr_srcptr, intmax_t, mpfr_rnd_t); +/* define mpfr_pown (defined in IEEE 754-2019) as an alias for mpfr_pow_sj */ +#define mpfr_pown mpfr_pow_sj #if defined (__cplusplus) } diff --git a/src/nbits_ulong.c b/src/nbits_ulong.c index 13bf11cd3..4f6d35110 100644 --- a/src/nbits_ulong.c +++ b/src/nbits_ulong.c @@ -1,4 +1,5 @@ /* mpfr_nbits_ulong -- number of significant bits in an unsigned long + mpfr_nbits_uj -- number of significant bits in an uintmax_t Copyright 2018-2021 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. @@ -21,6 +22,7 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ #define MPFR_NEED_LONGLONG_H /* for count_leading_zeros */ +#define MPFR_NEED_INTMAX_H #include "mpfr-impl.h" /* count the number of significant bits of n, i.e., @@ -81,3 +83,55 @@ mpfr_nbits_ulong (unsigned long n) MPFR_ASSERTD (cnt >= 0); return cnt; } + +#ifdef _MPFR_H_HAVE_INTMAX_T +/* count the number of significant bits of n, i.e., + nbits(uintmax_t) - count_leading_zeros (n) */ +int +mpfr_nbits_uj (uintmax_t n) +{ + int cnt; + + MPFR_ASSERTD (n > 0); + + cnt = 0; + + while (n >= 0x10000) + { + n >>= 16; + cnt += 16; + } + + MPFR_ASSERTD (n <= 0xffff); + + if (n >= 0x100) + { + n >>= 8; + cnt += 8; + } + + MPFR_ASSERTD (n <= 0xff); + + if (n >= 0x10) + { + n >>= 4; + cnt += 4; + } + + MPFR_ASSERTD (n <= 0xf); + + if (n >= 4) + { + n >>= 2; + cnt += 2; + } + + MPFR_ASSERTD (n <= 3); + + /* now n = 1, 2, or 3 */ + cnt += 1 + (n >= 2); + + MPFR_ASSERTD (cnt >= 0); + return cnt; +} +#endif diff --git a/src/pow_si.c b/src/pow_si.c index c17a19fa1..5bfd74775 100644 --- a/src/pow_si.c +++ b/src/pow_si.c @@ -23,21 +23,34 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_NEED_LONGLONG_H /* for MPFR_INT_CEIL_LOG2 */ #include "mpfr-impl.h" -/* The computation of y = pow_si(x,n) is done by - * y = pow_ui(x,n) if n >= 0 - * y = 1 / pow_ui(x,-n) if n < 0 +/* The computation of y = pow_si/sj(x,n) is done by + * y = pow_ui/uj(x,n) if n >= 0 + * y = 1 / pow_ui/uj(x,-n) if n < 0 */ +#ifndef POW_S +#define POW_S mpfr_pow_si +#define POW_U mpfr_pow_ui +#define SET_S mpfr_set_si +#define SET_S_2EXP mpfr_set_si_2exp +#define NBITS_UTYPE mpfr_nbits_ulong +#define TYPE long int +#define UTYPE unsigned long +#define LONG_TYPE 1 +#endif + int -mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) +POW_S (mpfr_ptr y, mpfr_srcptr x, TYPE n, mpfr_rnd_t rnd) { +#ifdef LONG_TYPE /* FIXME: what printf format for intmax_t? */ MPFR_LOG_FUNC (("x[%Pu]=%.*Rg n=%ld rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, n, rnd), ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y)); +#endif if (n >= 0) - return mpfr_pow_ui (y, x, n, rnd); + return POW_U (y, x, n, rnd); else { if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) @@ -49,7 +62,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) } else { - int positive = MPFR_IS_POS (x) || ((unsigned long) n & 1) == 0; + int positive = MPFR_IS_POS (x) || ((UTYPE) n & 1) == 0; if (MPFR_IS_INF (x)) MPFR_SET_ZERO (y); else /* x is zero */ @@ -130,8 +143,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) MPFR_EMIN_MIN - 2 /* Underflow */ : n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ? MPFR_EMAX_MAX /* Overflow */ : n * expx; - return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1, - expy, rnd); + return SET_S_2EXP (y, n % 2 ? MPFR_INT_SIGN (x) : 1, expy, rnd); } /* General case */ @@ -144,12 +156,12 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) mpfr_rnd_t rnd1; int size_n; int inexact; - unsigned long abs_n; + UTYPE abs_n; MPFR_SAVE_EXPO_DECL (expo); MPFR_ZIV_DECL (loop); - abs_n = - (unsigned long) n; - size_n = mpfr_nbits_ulong (abs_n); + abs_n = - (UTYPE) n; + size_n = NBITS_UTYPE (abs_n); /* initial working precision */ Ny = MPFR_PREC (y); @@ -177,7 +189,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */ if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags))) goto overflow; - MPFR_BLOCK (flags, mpfr_pow_ui (t, t, abs_n, rnd)); + MPFR_BLOCK (flags, POW_U (t, t, abs_n, rnd)); /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt). If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat} @@ -210,8 +222,8 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mpfr_rnd_t rnd) MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x), MPFR_EXP (x) - 1) != 0); mpfr_init2 (y2, 2); - mpfr_init2 (nn, sizeof (long) * CHAR_BIT); - inexact = mpfr_set_si (nn, n, MPFR_RNDN); + mpfr_init2 (nn, sizeof (TYPE) * CHAR_BIT); + inexact = SET_S (nn, n, MPFR_RNDN); MPFR_ASSERTN (inexact == 0); inexact = mpfr_pow_general (y2, x, nn, rnd, 1, (mpfr_save_expo_t *) NULL); diff --git a/src/pow_ui.c b/src/pow_ui.c index fba1e44ab..20e4de1c7 100644 --- a/src/pow_ui.c +++ b/src/pow_ui.c @@ -1,5 +1,4 @@ -/* mpfr_pow_ui-- compute the power of a floating-point - by a machine integer +/* mpfr_pow_ui -- compute the power of a floating-point by a machine integer Copyright 1999-2021 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. @@ -24,11 +23,18 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" +#ifndef POW_U +#define POW_U mpfr_pow_ui +#define MPZ_SET_U mpz_set_ui +#define UTYPE unsigned long int +#define LONG_TYPE 1 +#endif + /* sets y to x^n, and return 0 if exact, non-zero otherwise */ int -mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mpfr_rnd_t rnd) +POW_U (mpfr_ptr y, mpfr_srcptr x, UTYPE n, mpfr_rnd_t rnd) { - unsigned long m; + UTYPE m; mpfr_t res; mpfr_prec_t prec, err; int inexact; @@ -37,11 +43,13 @@ mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mpfr_rnd_t rnd) MPFR_ZIV_DECL (loop); MPFR_BLOCK_DECL (flags); +#ifdef LONG_TYPE /* FIXME: how to print a uintmax_t? */ MPFR_LOG_FUNC (("x[%Pu]=%.*Rg n=%lu rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, n, rnd), ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inexact)); +#endif /* x^0 = 1 for any x, even a NaN */ if (MPFR_UNLIKELY (n == 0)) @@ -110,12 +118,12 @@ mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mpfr_rnd_t rnd) MPFR_BLOCK (flags, inexact = mpfr_sqr (res, x, MPFR_RNDU); MPFR_ASSERTD (i >= 2); - if (n & (1UL << (i-2))) + if (n & ((UTYPE) 1 << (i-2))) inexact |= mpfr_mul (res, res, x, rnd1); for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--) { inexact |= mpfr_sqr (res, res, MPFR_RNDU); - if (n & (1UL << i)) + if (n & ((UTYPE) 1 << i)) inexact |= mpfr_mul (res, res, x, rnd1); }); /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2, @@ -150,7 +158,7 @@ mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mpfr_rnd_t rnd) mpfr_clear (res); MPFR_SAVE_EXPO_FREE (expo); mpz_init (z); - mpz_set_ui (z, n); + MPZ_SET_U (z, n); inexact = mpfr_pow_z (y, x, z, rnd); mpz_clear (z); return inexact; diff --git a/tests/tpow.c b/tests/tpow.c index 6a3192b21..50d04e4c7 100644 --- a/tests/tpow.c +++ b/tests/tpow.c @@ -21,6 +21,7 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ #define _MPFR_NO_DEPRECATED_ROOT +#define MPFR_NEED_INTMAX_H #include "mpfr-test.h" #ifdef CHECK_EXTERNAL @@ -360,6 +361,7 @@ check_pow_si (void) static void check_pown_ieee754_2019 (void) { +#ifdef _MPFR_H_HAVE_INTMAX_T mpfr_t x; mpfr_init2 (x, 5); /* ensures 17 is exact */ @@ -458,6 +460,7 @@ check_pown_ieee754_2019 (void) MPFR_ASSERTN(mpfr_zero_p (x) && mpfr_signbit (x) == 0); mpfr_clear (x); +#endif } static void |