summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorzimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4>2021-02-05 16:47:26 +0000
committerzimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4>2021-02-05 16:47:26 +0000
commit4f2a30795c961d88cddaa848691e46eedfa5e458 (patch)
tree257df99e5a29c0a44b9b3444911fda0b7e251ef6 /src
parent88db8636c893e759bda2c3df26e0b50618aed402 (diff)
downloadmpfr-4f2a30795c961d88cddaa848691e46eedfa5e458.tar.gz
added new function mpfr_exp10m1
git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/trunk@14388 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am2
-rw-r--r--src/exp10m1.c179
-rw-r--r--src/mpfr.h5
3 files changed, 183 insertions, 3 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index bfc68cd20..fea49b27a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -68,7 +68,7 @@ random_deviate.h random_deviate.c erandom.c mpfr-mini-gmp.c \
mpfr-mini-gmp.h fmma.c log_ui.c gamma_inc.c ubf.c invert_limb.h \
invsqrt_limb.h beta.c odd_p.c get_q.c pool.c total_order.c set_d128.c \
get_d128.c nbits_ulong.c cmpabs_ui.c sinu.c cosu.c tanu.c fmod_ui.c \
-acosu.c asinu.c atanu.c compound.c exp2m1.c
+acosu.c asinu.c atanu.c compound.c exp2m1.c exp10m1.c
nodist_libmpfr_la_SOURCES = $(BUILT_SOURCES)
diff --git a/src/exp10m1.c b/src/exp10m1.c
new file mode 100644
index 000000000..d5a7663b4
--- /dev/null
+++ b/src/exp10m1.c
@@ -0,0 +1,179 @@
+/* mpfr_exp10m1 -- Compute 10^x-1
+
+Copyright 2001-2021 Free Software Foundation, Inc.
+Contributed by the AriC and Caramba projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library 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 Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
+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
+#include "mpfr-impl.h"
+
+/* The computation of exp10m1 is done by expm1(x) = 10^x-1 */
+
+/* In case x is small in absolute value, 10^x - 1 ~ x*log(10).
+ If this is enough to deduce correct rounding, put in the auxiliary variable
+ t the approximation that will be rounded to get y, and return non-zero.
+ If we put 0 in t, it means underflow.
+ Otherwise return 0. */
+static int
+mpfr_exp10m1_small (mpfr_ptr y, mpfr_srcptr x , mpfr_rnd_t rnd_mode,
+ mpfr_ptr t)
+{
+ mpfr_prec_t prec;
+ mpfr_exp_t e;
+
+ /* for |x| < 0.25, we have |10^x-1-x*log(10)| < 4*x^2 */
+ if (MPFR_EXP(x) > -2)
+ return 0;
+ /* now EXP(x) <= -2, thus x < 0.25 */
+ prec = MPFR_PREC(t);
+ mpfr_log_ui (t, 10, MPFR_RNDN);
+ /* t = log(10)*(1 + theta) with |theta| <= 2^(-prec) */
+ mpfr_mul (t, t, x, MPFR_RNDN);
+ /* no underflow can occur, since log(10) > 1 */
+ /* t = x*log(10)*(1 + theta)^2 with |theta| <= 2^(-prec) */
+ /* |t - x*log(10)| <= ((1 + theta)^2 - 1) * |t| <= 3*2^(-prec)*|t| */
+ /* |t - x*log(10)| < 3*2^(EXP(t)-prec) */
+ e = 2 * MPFR_GET_EXP (x) + 2 + prec - MPFR_GET_EXP(t);
+ /* |4*x^2| < 2^e*2^(EXP(t)-prec) thus
+ |t - exp10m1(x)| < (3+2^e)*2^(EXP(t)-prec) */
+ e = (e <= 1) ? 2 + (e == 1) : e + 1;
+ /* now |t - exp10m1(x)| < 2^e*2^(EXP(t)-prec) */
+ return MPFR_CAN_ROUND (t, prec - e, MPFR_PREC(y), rnd_mode);
+}
+
+int
+mpfr_exp10m1 (mpfr_ptr y, mpfr_srcptr x , mpfr_rnd_t rnd_mode)
+{
+ int inexact, nloop;
+ mpfr_t t;
+ mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */
+ mpfr_prec_t Nt; /* working precision */
+ mpfr_exp_t err, exp_te; /* error */
+ MPFR_ZIV_DECL (loop);
+ MPFR_SAVE_EXPO_DECL (expo);
+
+ MPFR_LOG_FUNC
+ (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
+ ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
+ inexact));
+
+ if (MPFR_IS_SINGULAR (x))
+ return mpfr_expm1 (y, x, rnd_mode); /* singular cases are identical */
+
+ MPFR_ASSERTN(!MPFR_IS_ZERO(x));
+
+ MPFR_SAVE_EXPO_MARK (expo);
+
+ /* Check case where result is -1 or nextabove(-1) because x is a huge
+ negative number. */
+ if (MPFR_IS_NEG(x) && mpfr_cmpabs_ui (x, 2 + (MPFR_PREC(y) - 1) / 3) > 0)
+ {
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_INEXACT);
+ /* 1/2*ulp(-1) = 2^(-PREC(y)):
+ since |x| >= PREC(y)/3 + 1, then 3*|x| >= PREC(y) + 3,
+ thus 10^x < 8^x <= 2^(-PREC(y)-3) <= 1/2*ulp(-1), thus the
+ result is -1 for RNDA,RNDD,RNDN, and nextabove(-1) for RNDZ,RNDU */
+ mpfr_set_si (y, -1, MPFR_RNDZ);
+ if (!MPFR_IS_LIKE_RNDZ(rnd_mode,1))
+ inexact = -1;
+ else
+ {
+ mpfr_nextabove (y);
+ inexact = 1;
+ }
+ goto end;
+ }
+
+ Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6;
+
+ mpfr_init2 (t, Nt);
+
+ MPFR_ZIV_INIT (loop, Nt);
+ for (nloop = 0;; nloop++)
+ {
+ int inex1;
+
+ MPFR_BLOCK_DECL (flags);
+
+ /* 10^x may overflow and underflow */
+ MPFR_BLOCK (flags, inex1 = mpfr_exp10 (t, x, MPFR_RNDN));
+
+ if (MPFR_OVERFLOW (flags)) /* overflow case */
+ {
+ inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN_POS);
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
+ goto clear;
+ }
+
+ /* integer case */
+ if (inex1 == 0)
+ {
+ inexact = mpfr_sub_ui (y, t, 1, rnd_mode);
+ goto clear;
+ }
+
+ /* The case of underflow in 10^x (huge negative x)
+ was already detected before Ziv's loop. */
+ MPFR_ASSERTD(!MPFR_UNDERFLOW (flags));
+
+ MPFR_ASSERTN(!MPFR_IS_ZERO(t));
+ exp_te = MPFR_GET_EXP (t);
+ mpfr_sub_ui (t, t, 1, MPFR_RNDN); /* 10^x-1 */
+
+ /* error estimate */
+ /* err = __gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))) */
+ if (!MPFR_IS_ZERO(t))
+ {
+ err = MAX (exp_te - MPFR_GET_EXP (t), 0) + 1;
+ /* if inex1=0, this means that t=o(10^x) is exact, thus the correct
+ rounding is simply o(t-1) */
+ if (inex1 == 0 ||
+ MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Ny, rnd_mode)))
+ break;
+ }
+
+ /* check small case: we need to do it at each step of Ziv's loop,
+ since the multiplication x*log(10) might not enable correct
+ rounding at the first loop */
+ if (mpfr_exp10m1_small (y, x, rnd_mode, t))
+ {
+ if (MPFR_IS_ZERO(t)) /* underflow */
+ {
+ mpfr_clear (t);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_underflow (y, (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ
+ : rnd_mode, 1);
+ }
+ break;
+ }
+
+ /* increase the precision */
+ MPFR_ZIV_NEXT (loop, Nt);
+ mpfr_set_prec (t, Nt);
+ }
+
+ inexact = mpfr_set (y, t, rnd_mode);
+ clear:
+ MPFR_ZIV_FREE (loop);
+ mpfr_clear (t);
+
+ end:
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (y, inexact, rnd_mode);
+}
diff --git a/src/mpfr.h b/src/mpfr.h
index 083e4f918..e98c6530b 100644
--- a/src/mpfr.h
+++ b/src/mpfr.h
@@ -618,10 +618,11 @@ __MPFR_DECLSPEC int mpfr_log10p1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
__MPFR_DECLSPEC int mpfr_log_ui (mpfr_ptr, unsigned long, mpfr_rnd_t);
__MPFR_DECLSPEC int mpfr_exp (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
-__MPFR_DECLSPEC int mpfr_expm1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
__MPFR_DECLSPEC int mpfr_exp2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
-__MPFR_DECLSPEC int mpfr_exp2m1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
__MPFR_DECLSPEC int mpfr_exp10 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_expm1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_exp2m1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
+__MPFR_DECLSPEC int mpfr_exp10m1 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
__MPFR_DECLSPEC int mpfr_eint (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
__MPFR_DECLSPEC int mpfr_li2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);