1 /* mpfr_exp2 -- power of 2 function 2^y
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of y = 2^z is done by *
27 * y = exp(z*log(2)). The result is exact iff z is an integer. */
30 mpfr_exp2 (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
35 MPFR_SAVE_EXPO_DECL (expo
);
38 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x
), mpfr_log_prec
, x
, rnd_mode
),
39 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y
), mpfr_log_prec
, y
,
42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
49 else if (MPFR_IS_INF (x
))
60 MPFR_ASSERTD (MPFR_IS_ZERO(x
));
61 return mpfr_set_ui (y
, 1, rnd_mode
);
65 /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin,
66 if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */
67 MPFR_ASSERTN (MPFR_EMIN_MIN
>= LONG_MIN
+ 2);
68 if (MPFR_UNLIKELY (mpfr_cmp_si (x
, __gmpfr_emin
- 1) < 0))
70 mpfr_rnd_t rnd2
= rnd_mode
;
71 /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */
72 if (rnd_mode
== MPFR_RNDN
&&
73 mpfr_cmp_si_2exp (x
, __gmpfr_emin
- 2, 0) <= 0)
75 return mpfr_underflow (y
, rnd2
, 1);
78 MPFR_ASSERTN (MPFR_EMAX_MAX
<= LONG_MAX
);
79 if (MPFR_UNLIKELY (mpfr_cmp_si (x
, __gmpfr_emax
) >= 0))
80 return mpfr_overflow (y
, rnd_mode
, 1);
82 /* We now know that emin - 1 <= x < emax. */
84 MPFR_SAVE_EXPO_MARK (expo
);
86 /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have
87 |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1);
88 if x < 0 we must round toward 0 (dir=0). */
89 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y
, __gmpfr_one
, - MPFR_GET_EXP (x
), 0,
90 MPFR_SIGN(x
) > 0, rnd_mode
, expo
, {});
92 xint
= mpfr_get_si (x
, MPFR_RNDZ
);
93 mpfr_init2 (xfrac
, MPFR_PREC (x
));
94 mpfr_sub_si (xfrac
, x
, xint
, MPFR_RNDN
); /* exact */
96 if (MPFR_IS_ZERO (xfrac
))
98 mpfr_set_ui (y
, 1, MPFR_RNDN
);
103 /* Declaration of the intermediary variable */
106 /* Declaration of the size variable */
107 mpfr_prec_t Ny
= MPFR_PREC(y
); /* target precision */
108 mpfr_prec_t Nt
; /* working precision */
109 mpfr_exp_t err
; /* error */
110 MPFR_ZIV_DECL (loop
);
112 /* compute the precision of intermediary variable */
113 /* the optimal number of bits : see algorithms.tex */
114 Nt
= Ny
+ 5 + MPFR_INT_CEIL_LOG2 (Ny
);
116 /* initialise of intermediary variable */
119 /* First computation */
120 MPFR_ZIV_INIT (loop
, Nt
);
123 /* compute exp(x*ln(2))*/
124 mpfr_const_log2 (t
, MPFR_RNDU
); /* ln(2) */
125 mpfr_mul (t
, xfrac
, t
, MPFR_RNDU
); /* xfrac * ln(2) */
126 err
= Nt
- (MPFR_GET_EXP (t
) + 2); /* Estimate of the error */
127 mpfr_exp (t
, t
, MPFR_RNDN
); /* exp(xfrac * ln(2)) */
129 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, err
, Ny
, rnd_mode
)))
132 /* Actualisation of the precision */
133 MPFR_ZIV_NEXT (loop
, Nt
);
134 mpfr_set_prec (t
, Nt
);
136 MPFR_ZIV_FREE (loop
);
138 inexact
= mpfr_set (y
, t
, rnd_mode
);
145 mpfr_mul_2si (y
, y
, xint
, MPFR_RNDN
); /* exact or overflow */
146 /* Note: We can have an overflow only when t was rounded up to 2. */
147 MPFR_ASSERTD (MPFR_IS_PURE_FP (y
) || inexact
> 0);
148 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, __gmpfr_flags
);
149 MPFR_SAVE_EXPO_FREE (expo
);
150 return mpfr_check_range (y
, inexact
, rnd_mode
);