blob: 8db34ada536595cca3570ec1ef694c816e931bd7 [file] [log] [blame]
lh9ed821d2023-04-07 01:36:19 -07001/* Compute 2^x.
2 Copyright (C) 2012-2015 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <http://www.gnu.org/licenses/>. */
18
19#include <math.h>
20#include <math_private.h>
21#include <float.h>
22
23/* To avoid spurious underflows, use this definition to treat IBM long
24 double as approximating an IEEE-style format. */
25#if LDBL_MANT_DIG == 106
26# undef LDBL_EPSILON
27# define LDBL_EPSILON 0x1p-106L
28#endif
29
30long double
31__ieee754_exp2l (long double x)
32{
33 if (__glibc_likely (isless (x, (long double) LDBL_MAX_EXP)))
34 {
35 if (__builtin_expect (isgreaterequal (x, (long double) (LDBL_MIN_EXP
36 - LDBL_MANT_DIG
37 - 1)), 1))
38 {
39 int intx = (int) x;
40 long double fractx = x - intx;
41 if (fabsl (fractx) < LDBL_EPSILON / 4.0L)
42 return __scalbnl (1.0L + fractx, intx);
43 return __scalbnl (__ieee754_expl (M_LN2l * fractx), intx);
44 }
45 else
46 {
47 /* Underflow or exact zero. */
48 if (isinf (x))
49 return 0;
50 else
51 return LDBL_MIN * LDBL_MIN;
52 }
53 }
54 else
55 /* Infinity, NaN or overflow. */
56 return LDBL_MAX * x;
57}
58strong_alias (__ieee754_exp2l, __exp2l_finite)