1 |
|
|
/* $OpenBSD: s_exp2f.c,v 1.2 2012/12/05 23:20:04 deraadt Exp $ */ |
2 |
|
|
/*- |
3 |
|
|
* Copyright (c) 2005 David Schultz <das@FreeBSD.ORG> |
4 |
|
|
* All rights reserved. |
5 |
|
|
* |
6 |
|
|
* Redistribution and use in source and binary forms, with or without |
7 |
|
|
* modification, are permitted provided that the following conditions |
8 |
|
|
* are met: |
9 |
|
|
* 1. Redistributions of source code must retain the above copyright |
10 |
|
|
* notice, this list of conditions and the following disclaimer. |
11 |
|
|
* 2. Redistributions in binary form must reproduce the above copyright |
12 |
|
|
* notice, this list of conditions and the following disclaimer in the |
13 |
|
|
* documentation and/or other materials provided with the distribution. |
14 |
|
|
* |
15 |
|
|
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND |
16 |
|
|
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
17 |
|
|
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
18 |
|
|
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
19 |
|
|
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
20 |
|
|
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
21 |
|
|
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
22 |
|
|
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
23 |
|
|
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
24 |
|
|
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
25 |
|
|
* SUCH DAMAGE. |
26 |
|
|
*/ |
27 |
|
|
|
28 |
|
|
#include <float.h> |
29 |
|
|
|
30 |
|
|
#include "math.h" |
31 |
|
|
#include "math_private.h" |
32 |
|
|
|
33 |
|
|
#define TBLBITS 4 |
34 |
|
|
#define TBLSIZE (1 << TBLBITS) |
35 |
|
|
|
36 |
|
|
static const float |
37 |
|
|
huge = 0x1p100f, |
38 |
|
|
redux = 0x1.8p23f / TBLSIZE, |
39 |
|
|
P1 = 0x1.62e430p-1f, |
40 |
|
|
P2 = 0x1.ebfbe0p-3f, |
41 |
|
|
P3 = 0x1.c6b348p-5f, |
42 |
|
|
P4 = 0x1.3b2c9cp-7f; |
43 |
|
|
|
44 |
|
|
static volatile float twom100 = 0x1p-100f; |
45 |
|
|
|
46 |
|
|
static const double exp2ft[TBLSIZE] = { |
47 |
|
|
0x1.6a09e667f3bcdp-1, |
48 |
|
|
0x1.7a11473eb0187p-1, |
49 |
|
|
0x1.8ace5422aa0dbp-1, |
50 |
|
|
0x1.9c49182a3f090p-1, |
51 |
|
|
0x1.ae89f995ad3adp-1, |
52 |
|
|
0x1.c199bdd85529cp-1, |
53 |
|
|
0x1.d5818dcfba487p-1, |
54 |
|
|
0x1.ea4afa2a490dap-1, |
55 |
|
|
0x1.0000000000000p+0, |
56 |
|
|
0x1.0b5586cf9890fp+0, |
57 |
|
|
0x1.172b83c7d517bp+0, |
58 |
|
|
0x1.2387a6e756238p+0, |
59 |
|
|
0x1.306fe0a31b715p+0, |
60 |
|
|
0x1.3dea64c123422p+0, |
61 |
|
|
0x1.4bfdad5362a27p+0, |
62 |
|
|
0x1.5ab07dd485429p+0, |
63 |
|
|
}; |
64 |
|
|
|
65 |
|
|
/* |
66 |
|
|
* exp2f(x): compute the base 2 exponential of x |
67 |
|
|
* |
68 |
|
|
* Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927. |
69 |
|
|
* |
70 |
|
|
* Method: (equally-spaced tables) |
71 |
|
|
* |
72 |
|
|
* Reduce x: |
73 |
|
|
* x = 2**k + y, for integer k and |y| <= 1/2. |
74 |
|
|
* Thus we have exp2f(x) = 2**k * exp2(y). |
75 |
|
|
* |
76 |
|
|
* Reduce y: |
77 |
|
|
* y = i/TBLSIZE + z for integer i near y * TBLSIZE. |
78 |
|
|
* Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z), |
79 |
|
|
* with |z| <= 2**-(TBLSIZE+1). |
80 |
|
|
* |
81 |
|
|
* We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a |
82 |
|
|
* degree-4 minimax polynomial with maximum error under 1.4 * 2**-33. |
83 |
|
|
* Using double precision for everything except the reduction makes |
84 |
|
|
* roundoff error insignificant and simplifies the scaling step. |
85 |
|
|
* |
86 |
|
|
* This method is due to Tang, but I do not use his suggested parameters: |
87 |
|
|
* |
88 |
|
|
* Tang, P. Table-driven Implementation of the Exponential Function |
89 |
|
|
* in IEEE Floating-Point Arithmetic. TOMS 15(2), 144-157 (1989). |
90 |
|
|
*/ |
91 |
|
|
float |
92 |
|
|
exp2f(float x) |
93 |
|
|
{ |
94 |
|
|
double tv, twopk, u, z; |
95 |
|
|
float t; |
96 |
|
|
uint32_t hx, ix, i0; |
97 |
|
|
int32_t k; |
98 |
|
|
|
99 |
|
|
/* Filter out exceptional cases. */ |
100 |
|
|
GET_FLOAT_WORD(hx, x); |
101 |
|
|
ix = hx & 0x7fffffff; /* high word of |x| */ |
102 |
|
|
if(ix >= 0x43000000) { /* |x| >= 128 */ |
103 |
|
|
if(ix >= 0x7f800000) { |
104 |
|
|
if ((ix & 0x7fffff) != 0 || (hx & 0x80000000) == 0) |
105 |
|
|
return (x + x); /* x is NaN or +Inf */ |
106 |
|
|
else |
107 |
|
|
return (0.0); /* x is -Inf */ |
108 |
|
|
} |
109 |
|
|
if(x >= 0x1.0p7f) |
110 |
|
|
return (huge * huge); /* overflow */ |
111 |
|
|
if(x <= -0x1.2cp7f) |
112 |
|
|
return (twom100 * twom100); /* underflow */ |
113 |
|
|
} else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */ |
114 |
|
|
return (1.0f + x); |
115 |
|
|
} |
116 |
|
|
|
117 |
|
|
/* Reduce x, computing z, i0, and k. */ |
118 |
|
|
STRICT_ASSIGN(float, t, x + redux); |
119 |
|
|
GET_FLOAT_WORD(i0, t); |
120 |
|
|
i0 += TBLSIZE / 2; |
121 |
|
|
k = (i0 >> TBLBITS) << 20; |
122 |
|
|
i0 &= TBLSIZE - 1; |
123 |
|
|
t -= redux; |
124 |
|
|
z = x - t; |
125 |
|
|
INSERT_WORDS(twopk, 0x3ff00000 + k, 0); |
126 |
|
|
|
127 |
|
|
/* Compute r = exp2(y) = exp2ft[i0] * p(z). */ |
128 |
|
|
tv = exp2ft[i0]; |
129 |
|
|
u = tv * z; |
130 |
|
|
tv = tv + u * (P1 + z * P2) + u * (z * z) * (P3 + z * P4); |
131 |
|
|
|
132 |
|
|
/* Scale by 2**(k>>20). */ |
133 |
|
|
return (tv * twopk); |
134 |
|
|
} |