1/* Single-precision floating point 2^x.
2 Copyright (C) 1997-2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
5
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <http://www.gnu.org/licenses/>. */
19
20/* The basic design here is from
21 Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical
22 Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft.,
23 17 (1), March 1991, pp. 26-45.
24 It has been slightly modified to compute 2^x instead of e^x, and for
25 single-precision.
26 */
27#ifndef _GNU_SOURCE
28# define _GNU_SOURCE
29#endif
30#include <stdlib.h>
31#include <float.h>
32#include <ieee754.h>
33#include <math.h>
34#include <fenv.h>
35#include <inttypes.h>
36#include <math_private.h>
37
38#include "t_exp2f.h"
39
40static const float TWOM100 = 7.88860905e-31;
41static const float TWO127 = 1.7014118346e+38;
42
43float
44__ieee754_exp2f (float x)
45{
46 static const float himark = (float) FLT_MAX_EXP;
47 static const float lomark = (float) (FLT_MIN_EXP - FLT_MANT_DIG - 1);
48
49 /* Check for usual case. */
50 if (isless (x, himark) && isgreaterequal (x, lomark))
51 {
52 static const float THREEp14 = 49152.0;
53 int tval, unsafe;
54 float rx, x22, result;
55 union ieee754_float ex2_u, scale_u;
56
57 if (fabsf (x) < FLT_EPSILON / 4.0f)
58 return 1.0f + x;
59
60 {
61 SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
62
63 /* 1. Argument reduction.
64 Choose integers ex, -128 <= t < 128, and some real
65 -1/512 <= x1 <= 1/512 so that
66 x = ex + t/512 + x1.
67
68 First, calculate rx = ex + t/256. */
69 rx = x + THREEp14;
70 rx -= THREEp14;
71 x -= rx; /* Compute x=x1. */
72 /* Compute tval = (ex*256 + t)+128.
73 Now, t = (tval mod 256)-128 and ex=tval/256 [that's mod, NOT %;
74 and /-round-to-nearest not the usual c integer /]. */
75 tval = (int) (rx * 256.0f + 128.0f);
76
77 /* 2. Adjust for accurate table entry.
78 Find e so that
79 x = ex + t/256 + e + x2
80 where -7e-4 < e < 7e-4, and
81 (float)(2^(t/256+e))
82 is accurate to one part in 2^-64. */
83
84 /* 'tval & 255' is the same as 'tval%256' except that it's always
85 positive.
86 Compute x = x2. */
87 x -= __exp2f_deltatable[tval & 255];
88
89 /* 3. Compute ex2 = 2^(t/255+e+ex). */
90 ex2_u.f = __exp2f_atable[tval & 255];
91 tval >>= 8;
92 /* x2 is an integer multiple of 2^-30; avoid intermediate
93 underflow from the calculation of x22 * x. */
94 unsafe = abs(tval) >= -FLT_MIN_EXP - 32;
95 ex2_u.ieee.exponent += tval >> unsafe;
96 scale_u.f = 1.0;
97 scale_u.ieee.exponent += tval - (tval >> unsafe);
98
99 /* 4. Approximate 2^x2 - 1, using a second-degree polynomial,
100 with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14]
101 less than 1.3e-10. */
102
103 x22 = (.24022656679f * x + .69314736128f) * ex2_u.f;
104 }
105
106 /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
107 result = x22 * x + ex2_u.f;
108
109 if (!unsafe)
110 return result;
111 else
112 {
113 result *= scale_u.f;
114 math_check_force_underflow_nonneg (result);
115 return result;
116 }
117 }
118 /* Exceptional cases: */
119 else if (isless (x, himark))
120 {
121 if (isinf (x))
122 /* e^-inf == 0, with no error. */
123 return 0;
124 else
125 /* Underflow */
126 return TWOM100 * TWOM100;
127 }
128 else
129 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */
130 return TWO127*x;
131}
132strong_alias (__ieee754_exp2f, __exp2f_finite)
133