overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22//
23// ExpDD.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28// long double expl( long double x );
29// long double exp2l( long double x );
30// long double expm1l( long double x );
31//
32// _ExpInnerLD() is exported for use by other functions.
33//
34
35#include "math.h"
36#include "fp_private.h"
37#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
38#include "DD.h"
39
40// Floating-point constants
41
42static const double kLn2 = 6.9314718055994529e-1; // {0x3FE62E42, 0xFEFA39EF}
43static const double kLn2b = 2.3190468138462996e-17; // {0x3C7ABC9E, 0x3B39803F}
44static const double kLn2c = 5.7077084384162121e-34; // {0x3907B57A, 0x079A1934}
45static const double kLn2d = -3.5824322106018114e-50; // {0xB5AACE93, 0xA4EBE5D1}
46static const double k1ByLn2 = 1.4426950408889634; // {0x3FF71547, 0x652B82FE}
47static const double kBig = 6.755399441055744e+15; // {0x43380000, 0x00000000}
48static const double r13 = 1.6059043836821613e-10; // {0x3DE61246, 0x13A86D09}
49
50static const Double kInfinityu = {{0x7ff00000, 0x0}};
51
52static const double coeff[] = {
53 6227020800.0, 3113510400.0, 1037836800.0, 259459200.0, 51891840.0, 8648640.0,
54 1235520.0, 154440.0, 17160.0, 1716.0, 156.0, 13.0, 1.0
55};
56
57struct ExpTableEntry {
58 double x;
59 double fhead;
60 double ftail;
61};
62
63extern uint32_t ExpTableLD[];
64
65long double expl( long double x )
66{
67 double fpenv;
68 DblDbl u;
69 double extra;
70
71 u.f = x;
72
73 if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) // x = 0.0
74 return 1.0L;
75
76 if ((u.i[0] & 0x7ff00000u) != 0x7ff00000u) { // x is not NaN, Infinity
77 FEGETENVD(fpenv);
78 FESETENVD(0.0);
79 u.f = _ExpInnerLD(u.d[0], u.d[1], 0.0, &extra, 0);
80 if ((u.i[0] & 0x7ff00000) == 0x7ff00000)
81 u.d[1] = 0.0;
82 FESETENVD(fpenv);
83 return u.f;
84 }
85
86 if (u.d[0] != u.d[0]) // NaN case
87 return x;
88
89 if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) // +Infinity
90 return x;
91 else // -Infinity
92 return 0.0L;
93}
94
95long double exp2l( long double x )
96{
97 double fpenv;
98 DblDbl u;
99 double head, tail, extra, uu, vv, ww, xx, yy, c, top, bottom;
100
101 u.f = x;
102
103 if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) // x = 0.0
104 return 1.0L;
105
106 if ((u.i[0] & 0x7ff00000u) != 0x7ff00000u) { // x is not NaN, Infinity
107 FEGETENVD(fpenv);
108 FESETENVD(0.0);
109 head = u.d[0];
110 tail = u.d[1];
111 uu = __FADD( head * kLn2c, tail * kLn2b );
112 vv = head * kLn2b;
113 ww = __FMSUB( head, kLn2b, vv );
114 xx = tail * kLn2;
115 yy = __FMSUB( tail, kLn2, xx );
116 c = ww + yy + uu;
117 top = head * kLn2;
118 ww = __FMSUB( head, kLn2, top );
119
120 uu = __FADD( vv, xx );
121 if (fabs(xx) > fabs(vv))
122 c = xx - uu + vv + c;
123 else
124 c = vv - uu + xx + c;
125
126 bottom = __FADD( uu, ww );
127 if (fabs(ww) > fabs(uu))
128 c = ww - bottom + uu + c;
129 else
130 c = uu - bottom + ww + c;
131
132 head = __FADD( top, bottom );
133 ww = (top - head) + bottom;
134 tail = __FADD( ww, c );
135 c = (ww - tail) + c;
136
137 u.f = _ExpInnerLD(head, tail, c, &extra, 3);
138 if ((u.i[0] & 0x7ff00000) == 0x7ff00000) {
139 u.d[1] = 0.0;
140 }
141 FESETENVD(fpenv);
142 return u.f;
143 }
144
145 if (u.d[0] != u.d[0]) // NaN case
146 return x;
147 if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) // +Inifnity
148 return x;
149 else // -Infinity
150 return 0.0L;
151}
152
153long double expm1l( long double x )
154{
155 double fpenv;
156 DblDbl u;
157 double extra;
158
159 u.f = x;
160
161 if (((u.i[0] & 0x7fffffffu) | u.i[1]) == 0) // x = +/-0.0
162 return x; // changed 5-23-2006 to return x instead of 0.0L.
163
164 if ((u.i[0] & 0x7ff00000u) != 0x7ff00000u) { // x is not NaN, Infinity
165 FEGETENVD(fpenv);
166 FESETENVD(0.0);
167 if ((u.i[0] & 0x7ff00000u) >= 0x2ef00000u) // |x| > 2^(-110)
168 u.f = _ExpInnerLD(u.d[0], u.d[1], 0.0, &extra, 2);
169 if ((u.i[0] & 0x7ff00000) == 0x7ff00000) // If the result is NaN or Infinity, clear out the tail.
170 u.d[1] = 0.0;
171 FESETENVD(fpenv);
172 return u.f;
173 }
174
175 if (u.d[0] != u.d[0]) // NaN case
176 return x;
177 if ((u.i[0] & 0xfff00000u) == 0x7ff00000u) // +Infinity
178 return x;
179 else // -Infinity
180 return -1.0L;
181}
182
183//***********************************************************************************
184//
185// FUNCTION: long double _ExpInnerLD(double alpha, double beta, double gamma,
186// double *extra, int exptype)
187//
188// DESCRIPTION: This routine is called internally by the following functions:
189//
190// long double Exp(long double);
191// long double Exp2(long double);
192// long double ExpMinus1(long double);
193// long double Power(long double, long double);
194// family of long double hyperbolic functions;
195//
196// 1) exptype = 0 (called by Exp()):
197// on entry: (alpha, beta)
198// on exit: (head, tail) of Exp(alpha + beta)
199//
200// 2) exptype = 1 (called by hyperbolic functions):
201// on entry: (alpha, beta)
202// on exit: (head, tail, extra) of Exp(alpha + beta)/2.0
203//
204// 3) exptype = 2 (called by ExpMinus1()):
205// on entry: (alpha, beta)
206// on exit: (head, tail) of Exp(alpha + beta) - 1.0
207//
208// 4) exptype = 3 (called by Power(), Exp2()):
209// on entry: (alpha, beta, gamma)
210// on exit: (head, tail) of Exp(alpha + beta + gamma)
211//
212// 5) exptype = 4 (called by gamma functions):
213// on entry: (alpha, beta, gamma)
214// on exit: (head, tail, extra) of Exp(alpha + beta + gamma)
215//
216// This routine assumes that the rounding direction on entry
217// is round-to-nearest, and infinity, NaN and 0 have been
218// pre-filtered so that the input is a normal/denormal nonzero
219// values.
220//
221//***********************************************************************************
222
223long double _ExpInnerLD(double alpha, double beta, double gamma, double *extra,
224 int exptype)
225{
226 int notnormal, i;
227 double tail, factor, redarg, b, blow, redarg1, ca, carry;
228 double arg, arg2a, c, d, arg2, sum, temp, suml, hi, prod, residual;
229 double sum2, suml2, t1, t2, small, high, prodlow, bump, resnew;
230 double close, top, residual2, residual3, bottom;
231 Double test, accndx, en, rscale, scale, power;
232 DblDbl result;
233 struct ExpTableEntry *TblPtr = (struct ExpTableEntry *)ExpTableLD + 24;
234
235 test.f = alpha;
236 scale.f = 0.0;
237 en.f = alpha*k1ByLn2 + kBig;
238 factor = en.f - kBig;
239 redarg = __FNMSUB( kLn2, factor, alpha ); // alpha - kLn2*factor;
240 notnormal = 0;
241 if ((test.i[0] & 0x7fffffffu) > 0x40862000u) { // |alpha| > 708
242 if (alpha > 0) {
243 if (alpha <= 1026.0 * kLn2) {
244 en.i[1] -= 4;
245 notnormal = (1023 + 4) * 1048576;
246 }
247 else { // result overflow
248 if ((exptype == 1) || (exptype == 4))
249 *extra = 0.0;
250 result.d[0] = kInfinityu.f;
251 result.d[1] = 0.0;
252 return result.f;
253 }
254 }
255 else {
256 if ((alpha >= -800.0) && (exptype != 2)) {
257 en.i[1] += 256;
258 notnormal = (1023 - 256)*1048576;
259 }
260 else { // result underflow or (alpha <= -708, and exptype = 2 or 4)
261 if ((exptype == 1) || (exptype == 4))
262 *extra = 0.0;
263 result.d[0] = (exptype == 2) ? - 1.0 : 0.0;
264 result.d[1] = 0.0;
265 return result.f;
266 }
267 }
268 }
269 b = kLn2b*factor;
270 blow = __FMSUB( kLn2b, factor, b ); // kLn2b*factor - b;
271 redarg1 = __FNMSUB( kLn2b, factor, beta ); // beta - b;
272 accndx.f = redarg*64.0 + kBig;
273 ca = kLn2c*factor;
274 if (fabs(b) > fabs(beta))
275 carry = beta - (b + redarg1) - blow;
276 else
277 carry = (beta - redarg1) - b - blow;
278 redarg -= TblPtr[(int32_t)accndx.i[1]].x;
279 arg = __FADD( redarg, redarg1 );
280 arg2a = (redarg - arg + redarg1);
281 c = __FSUB( carry, ca );
282 d = __FNMSUB( kLn2c, factor, ca ) - (ca - (carry - c)) - kLn2d*factor; // (ca - kLn2c*factor) - (ca - (carry - c)) - kLn2d*factor;
283 scale.i[0] = (en.i[1] + 1023) << 20;
284 arg2 = arg2a + c + d;
285 if (exptype >= 3)
286 arg2 += gamma; // extpye = 3 or 4
287
288 sum = coeff[12];
289 sum = __FMADD( sum, arg, coeff[11] ); // coeff[11] + sum*arg;
290 sum = __FMADD( sum, arg, coeff[10] ); // coeff[10] + sum*arg;
291 sum = __FMADD( sum, arg, coeff[9] ); // coeff[ 9] + sum*arg;
292 sum = __FMADD( sum, arg, coeff[8] ); // coeff[ 8] + sum*arg;
293 temp = __FMADD( sum, arg, coeff[7] ); // coeff[7] + sum*arg;
294 suml = __FMADD( sum, arg, coeff[7] - temp ); // coeff[7] - temp + sum*arg;
295 sum = temp;
296 for (i = 6; i > 0; i--) {
297 hi = __FMADD( sum, arg, coeff[i] ); // coeff[i] + sum*arg;
298 suml = __FMADD( suml, arg, __FMADD( sum, arg, coeff[i] - hi ) ); // coeff[i] - hi + sum*arg + suml*arg;
299 sum = hi;
300 }
301 prod = sum*r13;
302 residual = __FNMSUB( coeff[0], prod, sum ) + suml; // (sum - coeff[0]*prod) + suml;
303 suml = residual*r13;
304 sum2 = __FMADD( prod, arg, 1.0); // __FADD( 1.0, prod*arg );
305 suml2 = __FMADD( suml, arg, __FMADD( prod, arg, 1.0 - sum2 ) ); // 1.0 - sum2 + prod*arg + suml*arg;
306 sum = __FMADD( suml2, arg, sum2*arg ); // sum2*arg + (suml2*arg);
307 suml = __FMADD( suml2, arg, __FMSUB( sum2, arg, sum ) ); // (sum2*arg - sum) + (suml2*arg);
308
309 t1 = TblPtr[(int32_t)accndx.i[1]].fhead;
310 t2 = TblPtr[(int32_t)accndx.i[1]].ftail;
311 small = t2 + sum*(t2 + arg2*t1);
312 prod = __FMUL( t1, sum );
313 high = __FADD( prod, t1 );
314 prodlow = __FMSUB( t1, sum, prod );
315 residual = (t1 - high) + prod;
316
317 if (exptype == 2) {
318 rscale.f = 0.0;
319 rscale.i[0] = 2046*1048576 - scale.i[0];
320 temp = __FSUB( high, rscale.f );
321 if (temp > 0)
322 bump = high - temp - rscale.f;
323 else
324 bump = high - (temp + rscale.f);
325 resnew = __FADD( residual, bump );
326 if (fabs(residual) > fabs(bump))
327 small = ((residual - resnew) + bump) + small;
328 else
329 small = ((bump - resnew) + residual) + small;
330 high = temp;
331 residual = resnew;
332 }
333
334 tail = small + arg2*t1 + prodlow + suml*t1;
335 close = __FADD( tail, residual );
336 top = __FADD( high, close );
337 residual2 = (residual - close) + tail;
338 residual3 = (high - top) + close;
339 bottom = __FADD( residual3, residual2 );
340
341 if ((exptype == 1) || (exptype == 4)) {
342 if (exptype == 1) scale.f *= 0.5;
343 if (fabs(residual2) > fabs(residual3))
344 *extra = (residual2 - bottom + residual3)*scale.f;
345 else
346 *extra = (residual3 - bottom + residual2)*scale.f;
347 }
348
349 if (notnormal != 0) {
350 power.f = 0.0;
351 power.i[0] = notnormal;
352 top *= power.f;
353 bottom *= power.f;
354 if ((exptype == 1) || (exptype == 4))
355 *extra *= power.f;
356 }
357
358 result.d[0] = top*scale.f;
359 result.d[1] = bottom*scale.f;
360 return result.f;
361}
362#endif