1 |
|
|
/* $OpenBSD: s_expm1l.c,v 1.3 2016/09/12 19:47:03 guenther Exp $ */ |
2 |
|
|
|
3 |
|
|
/* |
4 |
|
|
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
5 |
|
|
* |
6 |
|
|
* Permission to use, copy, modify, and distribute this software for any |
7 |
|
|
* purpose with or without fee is hereby granted, provided that the above |
8 |
|
|
* copyright notice and this permission notice appear in all copies. |
9 |
|
|
* |
10 |
|
|
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
11 |
|
|
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
12 |
|
|
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR |
13 |
|
|
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
14 |
|
|
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN |
15 |
|
|
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF |
16 |
|
|
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. |
17 |
|
|
*/ |
18 |
|
|
|
19 |
|
|
/* expm1l.c |
20 |
|
|
* |
21 |
|
|
* Exponential function, minus 1 |
22 |
|
|
* Long double precision |
23 |
|
|
* |
24 |
|
|
* |
25 |
|
|
* SYNOPSIS: |
26 |
|
|
* |
27 |
|
|
* long double x, y, expm1l(); |
28 |
|
|
* |
29 |
|
|
* y = expm1l( x ); |
30 |
|
|
* |
31 |
|
|
* |
32 |
|
|
* |
33 |
|
|
* DESCRIPTION: |
34 |
|
|
* |
35 |
|
|
* Returns e (2.71828...) raised to the x power, minus 1. |
36 |
|
|
* |
37 |
|
|
* Range reduction is accomplished by separating the argument |
38 |
|
|
* into an integer k and fraction f such that |
39 |
|
|
* |
40 |
|
|
* x k f |
41 |
|
|
* e = 2 e. |
42 |
|
|
* |
43 |
|
|
* An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1 |
44 |
|
|
* in the basic range [-0.5 ln 2, 0.5 ln 2]. |
45 |
|
|
* |
46 |
|
|
* |
47 |
|
|
* ACCURACY: |
48 |
|
|
* |
49 |
|
|
* Relative error: |
50 |
|
|
* arithmetic domain # trials peak rms |
51 |
|
|
* IEEE -45,+MAXLOG 200,000 1.2e-19 2.5e-20 |
52 |
|
|
* |
53 |
|
|
* ERROR MESSAGES: |
54 |
|
|
* |
55 |
|
|
* message condition value returned |
56 |
|
|
* expm1l overflow x > MAXLOG MAXNUM |
57 |
|
|
* |
58 |
|
|
*/ |
59 |
|
|
|
60 |
|
|
#include <math.h> |
61 |
|
|
|
62 |
|
|
static const long double MAXLOGL = 1.1356523406294143949492E4L; |
63 |
|
|
|
64 |
|
|
/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) |
65 |
|
|
-.5 ln 2 < x < .5 ln 2 |
66 |
|
|
Theoretical peak relative error = 3.4e-22 */ |
67 |
|
|
|
68 |
|
|
static const long double |
69 |
|
|
P0 = -1.586135578666346600772998894928250240826E4L, |
70 |
|
|
P1 = 2.642771505685952966904660652518429479531E3L, |
71 |
|
|
P2 = -3.423199068835684263987132888286791620673E2L, |
72 |
|
|
P3 = 1.800826371455042224581246202420972737840E1L, |
73 |
|
|
P4 = -5.238523121205561042771939008061958820811E-1L, |
74 |
|
|
|
75 |
|
|
Q0 = -9.516813471998079611319047060563358064497E4L, |
76 |
|
|
Q1 = 3.964866271411091674556850458227710004570E4L, |
77 |
|
|
Q2 = -7.207678383830091850230366618190187434796E3L, |
78 |
|
|
Q3 = 7.206038318724600171970199625081491823079E2L, |
79 |
|
|
Q4 = -4.002027679107076077238836622982900945173E1L, |
80 |
|
|
/* Q5 = 1.000000000000000000000000000000000000000E0 */ |
81 |
|
|
|
82 |
|
|
/* C1 + C2 = ln 2 */ |
83 |
|
|
C1 = 6.93145751953125E-1L, |
84 |
|
|
C2 = 1.428606820309417232121458176568075500134E-6L, |
85 |
|
|
/* ln 2^-65 */ |
86 |
|
|
minarg = -4.5054566736396445112120088E1L; |
87 |
|
|
static const long double huge = 0x1p10000L; |
88 |
|
|
|
89 |
|
|
long double |
90 |
|
|
expm1l(long double x) |
91 |
|
|
{ |
92 |
|
|
long double px, qx, xx; |
93 |
|
|
int k; |
94 |
|
|
|
95 |
|
|
/* Overflow. */ |
96 |
|
|
if (x > MAXLOGL) |
97 |
|
|
return (huge*huge); /* overflow */ |
98 |
|
|
|
99 |
|
|
if (x == 0.0) |
100 |
|
|
return x; |
101 |
|
|
|
102 |
|
|
/* Minimum value. */ |
103 |
|
|
if (x < minarg) |
104 |
|
|
return -1.0L; |
105 |
|
|
|
106 |
|
|
xx = C1 + C2; |
107 |
|
|
|
108 |
|
|
/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ |
109 |
|
|
px = floorl (0.5 + x / xx); |
110 |
|
|
k = px; |
111 |
|
|
/* remainder times ln 2 */ |
112 |
|
|
x -= px * C1; |
113 |
|
|
x -= px * C2; |
114 |
|
|
|
115 |
|
|
/* Approximate exp(remainder ln 2). */ |
116 |
|
|
px = (((( P4 * x |
117 |
|
|
+ P3) * x |
118 |
|
|
+ P2) * x |
119 |
|
|
+ P1) * x |
120 |
|
|
+ P0) * x; |
121 |
|
|
|
122 |
|
|
qx = (((( x |
123 |
|
|
+ Q4) * x |
124 |
|
|
+ Q3) * x |
125 |
|
|
+ Q2) * x |
126 |
|
|
+ Q1) * x |
127 |
|
|
+ Q0; |
128 |
|
|
|
129 |
|
|
xx = x * x; |
130 |
|
|
qx = x + (0.5 * xx + xx * px / qx); |
131 |
|
|
|
132 |
|
|
/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). |
133 |
|
|
We have qx = exp(remainder ln 2) - 1, so |
134 |
|
|
exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */ |
135 |
|
|
px = ldexpl(1.0L, k); |
136 |
|
|
x = px * qx + (px - 1.0); |
137 |
|
|
return x; |
138 |
|
|
} |
139 |
|
|
DEF_STD(expm1l); |