5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 /* INDENT OFF */
31 /*
32 * long double __k_cexpl(long double x, int *n);
33 * Returns the exponential of x in the form of 2**n * y, y=__k_cexpl(x,&n).
34 *
35 * 1. Argument Reduction: given the input x, find r and integer k
36 * and j such that
37 * x = (32k+j)*ln2 + r, |r| <= (1/64)*ln2 .
38 *
39 * 2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
40 * Note:
41 * a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
42 * b. 2^(j/32) is represented as
43 * exp2_32_hi[j]+exp2_32_lo[j]
44 * where
45 * exp2_32_hi[j] = 2^(j/32) rounded
46 * exp2_32_lo[j] = 2^(j/32) - exp2_32_hi[j].
47 *
48 * Special cases:
49 * expl(INF) is INF, expl(NaN) is NaN;
50 * expl(-INF)= 0;
51 * for finite argument, only expl(0)=1 is exact.
52 *
53 * Accuracy:
54 * according to an error analysis, the error is always less than
55 * an ulp (unit in the last place).
56 *
57 * Misc. info.
58 * When |x| is really big, say |x| > 1000000, the accuracy
59 * is not important because the ultimate result will over or under
60 * flow. So we will simply replace n = 1000000 and r = 0.0. For
61 * moderate size x, according to an error analysis, the error is
62 * always less than 1 ulp (unit in the last place).
63 *
64 * Constants:
65 * Only decimal values are given. We assume that the compiler will convert
66 * from decimal to binary accurately enough to produce the correct
67 * hexadecimal values.
68 */
69 /* INDENT ON */
70
71 #include "libm.h" /* __k_cexpl */
72 #include "complex_wrapper.h" /* HI_XWORD */
73
74 /* INDENT OFF */
75 /* ln2/32 = 0.0216608493924982909192885037955680177523593791987579766912713 */
76 #if defined(__x86)
77 static const long double
78 /* 43 significant bits, 21 trailing zeros */
79 ln2_32hi = 2.166084939249657281834515742957592010498046875e-2L,
80 ln2_32lo = 1.7181009433463659920976473789104487579766912713e-15L;
81 static const long double exp2_32_hi[] = { /* exp2_32[j] = 2^(j/32) */
82 1.0000000000000000000000000e+00L,
83 1.0218971486541166782081522e+00L,
84 1.0442737824274138402382006e+00L,
85 1.0671404006768236181297224e+00L,
86 1.0905077326652576591003302e+00L,
87 1.1143867425958925362894369e+00L,
88 1.1387886347566916536971221e+00L,
89 1.1637248587775775137938619e+00L,
90 1.1892071150027210666875674e+00L,
91 1.2152473599804688780476325e+00L,
92 1.2418578120734840485256747e+00L,
93 1.2690509571917332224885722e+00L,
94 1.2968395546510096659215822e+00L,
95 1.3252366431597412945939118e+00L,
96 1.3542555469368927282668852e+00L,
97 1.3839098819638319548151403e+00L,
98 1.4142135623730950487637881e+00L,
99 1.4451808069770466200253470e+00L,
100 1.4768261459394993113155431e+00L,
101 1.5091644275934227397133885e+00L,
102 1.5422108254079408235859630e+00L,
103 1.5759808451078864864006862e+00L,
104 1.6104903319492543080837174e+00L,
105 1.6457554781539648445110730e+00L,
106 1.6817928305074290860378350e+00L,
107 1.7186192981224779156032914e+00L,
108 1.7562521603732994831094730e+00L,
109 1.7947090750031071864148413e+00L,
110 1.8340080864093424633989166e+00L,
111 1.8741676341102999013002103e+00L,
112 1.9152065613971472938202589e+00L,
113 1.9571441241754002689657438e+00L,
114 };
115 static const long double exp2_32_lo[] = {
116 0.0000000000000000000000000e+00L,
117 2.6327965667180882569382524e-20L,
118 8.3765863521895191129661899e-20L,
119 3.9798705777454504249209575e-20L,
120 1.0668046596651558640993042e-19L,
121 1.9376009847285360448117114e-20L,
122 6.7081819456112953751277576e-21L,
123 1.9711680502629186462729727e-20L,
124 2.9932584438449523689104569e-20L,
125 6.8887754153039109411061914e-20L,
126 6.8002718741225378942847820e-20L,
127 6.5846917376975403439742349e-20L,
128 1.2171958727511372194876001e-20L,
129 3.5625253228704087115438260e-20L,
130 3.1129551559077560956309179e-20L,
131 5.7519192396164779846216492e-20L,
132 3.7900651177865141593101239e-20L,
133 1.1659262405698741798080115e-20L,
134 7.1364385105284695967172478e-20L,
135 5.2631003710812203588788949e-20L,
136 2.6328853788732632868460580e-20L,
137 5.4583950085438242788190141e-20L,
138 9.5803254376938269960718656e-20L,
139 7.6837733983874245823512279e-21L,
140 2.4415965910835093824202087e-20L,
141 2.6052966871016580981769728e-20L,
142 2.6876456344632553875309579e-21L,
143 1.2861930155613700201703279e-20L,
144 8.8166633394037485606572294e-20L,
145 2.9788615389580190940837037e-20L,
146 5.2352341619805098677422139e-20L,
147 5.2578463064010463732242363e-20L,
148 };
149 #else /* sparc */
150 static const long double
151 /* 0x3FF962E4 2FEFA39E F35793C7 00000000 */
152 ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
153 ln2_32lo = 5.209643502595475652782654157501186731779e-27L;
154 static const long double exp2_32_hi[] = { /* exp2_32[j] = 2^(j/32) */
155 1.000000000000000000000000000000000000000e+0000L,
156 1.021897148654116678234480134783299439782e+0000L,
157 1.044273782427413840321966478739929008785e+0000L,
158 1.067140400676823618169521120992809162607e+0000L,
159 1.090507732665257659207010655760707978993e+0000L,
160 1.114386742595892536308812956919603067800e+0000L,
161 1.138788634756691653703830283841511254720e+0000L,
162 1.163724858777577513813573599092185312343e+0000L,
163 1.189207115002721066717499970560475915293e+0000L,
164 1.215247359980468878116520251338798457624e+0000L,
165 1.241857812073484048593677468726595605511e+0000L,
166 1.269050957191733222554419081032338004715e+0000L,
167 1.296839554651009665933754117792451159835e+0000L,
168 1.325236643159741294629537095498721674113e+0000L,
169 1.354255546936892728298014740140702804343e+0000L,
170 1.383909881963831954872659527265192818002e+0000L,
171 1.414213562373095048801688724209698078570e+0000L,
172 1.445180806977046620037006241471670905678e+0000L,
173 1.476826145939499311386907480374049923924e+0000L,
205 -5.182484853064646457536893018566956189817e-0035L,
206 +9.422242548621832065692116736394064879758e-0035L,
207 -3.967500825398862309167306130216418281103e-0035L,
208 +7.143528991563300614523273615092767243521e-0035L,
209 +1.159871252867985124246517834100444327747e-0035L,
210 +4.696933478358115495309739213201874466685e-0035L,
211 -3.386513175995004710799241984999819165197e-0035L,
212 -8.587318774298247068868655935103874453522e-0035L,
213 -9.605951548749350503185499362246069088835e-0035L,
214 +9.609733932128012784507558697141785813655e-0035L,
215 +6.378397921440028439244761449780848545957e-0035L,
216 +7.792430785695864249456461125169277701177e-0035L,
217 +7.361337767588456524131930836633932195088e-0035L,
218 -6.472995147913347230035214575612170525266e-0035L,
219 +8.587474417953698694278798062295229624207e-0035L,
220 +2.371815422825174835691651228302690977951e-0035L,
221 -3.026891682096118773004597373421900314256e-0037L,
222 };
223 #endif
224
225 static const long double
226 one = 1.0L,
227 two = 2.0L,
228 ln2_64 = 1.083042469624914545964425189778400898568e-2L,
229 invln2_32 = 4.616624130844682903551758979206054839765e+1L;
230
231 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
232 static const long double
233 t1 = 1.666666666666666666666666666660876387437e-1L,
234 t2 = -2.777777777777777777777707812093173478756e-3L,
235 t3 = 6.613756613756613482074280932874221202424e-5L,
236 t4 = -1.653439153392139954169609822742235851120e-6L,
237 t5 = 4.175314851769539751387852116610973796053e-8L;
238 /* INDENT ON */
239
240 long double
241 __k_cexpl(long double x, int *n) {
242 int hx, ix, j, k;
243 long double t, r;
244
245 *n = 0;
246 hx = HI_XWORD(x);
247 ix = hx & 0x7fffffff;
248 if (hx >= 0x7fff0000)
249 return (x + x); /* NaN of +inf */
250 if (((unsigned) hx) >= 0xffff0000)
251 return (-one / x); /* NaN or -inf */
252 if (ix < 0x3fc30000)
253 return (one + x); /* |x|<2^-60 */
254 if (hx > 0) {
255 if (hx > 0x401086a0) { /* x > 200000 */
256 *n = 200000;
257 return (one);
258 }
259 k = (int) (invln2_32 * (x + ln2_64));
260 } else {
261 if (ix > 0x401086a0) { /* x < -200000 */
262 *n = -200000;
263 return (one);
264 }
265 k = (int) (invln2_32 * (x - ln2_64));
266 }
267 j = k & 0x1f;
268 *n = k >> 5;
269 t = (long double) k;
270 x = (x - t * ln2_32hi) - t * ln2_32lo;
271 t = x * x;
272 r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
273 x = exp2_32_hi[j] - ((exp2_32_hi[j] * (x + x)) / r - exp2_32_lo[j]);
274 k >>= 5;
275 if (k > 240) {
276 XFSCALE(x, 240);
277 *n -= 240;
278 } else if (k > 0) {
279 XFSCALE(x, k);
280 *n = 0;
281 }
282 return (x);
283 }
|
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31
32 /*
33 * long double __k_cexpl(long double x, int *n);
34 * Returns the exponential of x in the form of 2**n * y, y=__k_cexpl(x,&n).
35 *
36 * 1. Argument Reduction: given the input x, find r and integer k
37 * and j such that
38 * x = (32k+j)*ln2 + r, |r| <= (1/64)*ln2 .
39 *
40 * 2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
41 * Note:
42 * a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
43 * b. 2^(j/32) is represented as
44 * exp2_32_hi[j]+exp2_32_lo[j]
45 * where
46 * exp2_32_hi[j] = 2^(j/32) rounded
47 * exp2_32_lo[j] = 2^(j/32) - exp2_32_hi[j].
48 *
49 * Special cases:
50 * expl(INF) is INF, expl(NaN) is NaN;
51 * expl(-INF)= 0;
52 * for finite argument, only expl(0)=1 is exact.
53 *
54 * Accuracy:
55 * according to an error analysis, the error is always less than
56 * an ulp (unit in the last place).
57 *
58 * Misc. info.
59 * When |x| is really big, say |x| > 1000000, the accuracy
60 * is not important because the ultimate result will over or under
61 * flow. So we will simply replace n = 1000000 and r = 0.0. For
62 * moderate size x, according to an error analysis, the error is
63 * always less than 1 ulp (unit in the last place).
64 *
65 * Constants:
66 * Only decimal values are given. We assume that the compiler will convert
67 * from decimal to binary accurately enough to produce the correct
68 * hexadecimal values.
69 */
70
71 #include "libm.h" /* __k_cexpl */
72 #include "complex_wrapper.h" /* HI_XWORD */
73
74 /*
75 * ln2/32 = 0.0216608493924982909192885037955680177523593791987579766912713
76 */
77 #if defined(__x86)
78 static const long double
79 /* 43 significant bits, 21 trailing zeros */
80 ln2_32hi = 2.166084939249657281834515742957592010498046875e-2L,
81 ln2_32lo = 1.7181009433463659920976473789104487579766912713e-15L;
82
83 static const long double exp2_32_hi[] = { /* exp2_32[j] = 2^(j/32) */
84 1.0000000000000000000000000e+00L, 1.0218971486541166782081522e+00L,
85 1.0442737824274138402382006e+00L, 1.0671404006768236181297224e+00L,
86 1.0905077326652576591003302e+00L, 1.1143867425958925362894369e+00L,
87 1.1387886347566916536971221e+00L, 1.1637248587775775137938619e+00L,
88 1.1892071150027210666875674e+00L, 1.2152473599804688780476325e+00L,
89 1.2418578120734840485256747e+00L, 1.2690509571917332224885722e+00L,
90 1.2968395546510096659215822e+00L, 1.3252366431597412945939118e+00L,
91 1.3542555469368927282668852e+00L, 1.3839098819638319548151403e+00L,
92 1.4142135623730950487637881e+00L, 1.4451808069770466200253470e+00L,
93 1.4768261459394993113155431e+00L, 1.5091644275934227397133885e+00L,
94 1.5422108254079408235859630e+00L, 1.5759808451078864864006862e+00L,
95 1.6104903319492543080837174e+00L, 1.6457554781539648445110730e+00L,
96 1.6817928305074290860378350e+00L, 1.7186192981224779156032914e+00L,
97 1.7562521603732994831094730e+00L, 1.7947090750031071864148413e+00L,
98 1.8340080864093424633989166e+00L, 1.8741676341102999013002103e+00L,
99 1.9152065613971472938202589e+00L, 1.9571441241754002689657438e+00L,
100 };
101
102 static const long double exp2_32_lo[] = {
103 0.0000000000000000000000000e+00L, 2.6327965667180882569382524e-20L,
104 8.3765863521895191129661899e-20L, 3.9798705777454504249209575e-20L,
105 1.0668046596651558640993042e-19L, 1.9376009847285360448117114e-20L,
106 6.7081819456112953751277576e-21L, 1.9711680502629186462729727e-20L,
107 2.9932584438449523689104569e-20L, 6.8887754153039109411061914e-20L,
108 6.8002718741225378942847820e-20L, 6.5846917376975403439742349e-20L,
109 1.2171958727511372194876001e-20L, 3.5625253228704087115438260e-20L,
110 3.1129551559077560956309179e-20L, 5.7519192396164779846216492e-20L,
111 3.7900651177865141593101239e-20L, 1.1659262405698741798080115e-20L,
112 7.1364385105284695967172478e-20L, 5.2631003710812203588788949e-20L,
113 2.6328853788732632868460580e-20L, 5.4583950085438242788190141e-20L,
114 9.5803254376938269960718656e-20L, 7.6837733983874245823512279e-21L,
115 2.4415965910835093824202087e-20L, 2.6052966871016580981769728e-20L,
116 2.6876456344632553875309579e-21L, 1.2861930155613700201703279e-20L,
117 8.8166633394037485606572294e-20L, 2.9788615389580190940837037e-20L,
118 5.2352341619805098677422139e-20L, 5.2578463064010463732242363e-20L,
119 };
120 #else /* sparc */
121 static const long double
122 /* 0x3FF962E4 2FEFA39E F35793C7 00000000 */
123 ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
124 ln2_32lo = 5.209643502595475652782654157501186731779e-27L;
125 static const long double exp2_32_hi[] = { /* exp2_32[j] = 2^(j/32) */
126 1.000000000000000000000000000000000000000e+0000L,
127 1.021897148654116678234480134783299439782e+0000L,
128 1.044273782427413840321966478739929008785e+0000L,
129 1.067140400676823618169521120992809162607e+0000L,
130 1.090507732665257659207010655760707978993e+0000L,
131 1.114386742595892536308812956919603067800e+0000L,
132 1.138788634756691653703830283841511254720e+0000L,
133 1.163724858777577513813573599092185312343e+0000L,
134 1.189207115002721066717499970560475915293e+0000L,
135 1.215247359980468878116520251338798457624e+0000L,
136 1.241857812073484048593677468726595605511e+0000L,
137 1.269050957191733222554419081032338004715e+0000L,
138 1.296839554651009665933754117792451159835e+0000L,
139 1.325236643159741294629537095498721674113e+0000L,
140 1.354255546936892728298014740140702804343e+0000L,
141 1.383909881963831954872659527265192818002e+0000L,
142 1.414213562373095048801688724209698078570e+0000L,
143 1.445180806977046620037006241471670905678e+0000L,
144 1.476826145939499311386907480374049923924e+0000L,
176 -5.182484853064646457536893018566956189817e-0035L,
177 +9.422242548621832065692116736394064879758e-0035L,
178 -3.967500825398862309167306130216418281103e-0035L,
179 +7.143528991563300614523273615092767243521e-0035L,
180 +1.159871252867985124246517834100444327747e-0035L,
181 +4.696933478358115495309739213201874466685e-0035L,
182 -3.386513175995004710799241984999819165197e-0035L,
183 -8.587318774298247068868655935103874453522e-0035L,
184 -9.605951548749350503185499362246069088835e-0035L,
185 +9.609733932128012784507558697141785813655e-0035L,
186 +6.378397921440028439244761449780848545957e-0035L,
187 +7.792430785695864249456461125169277701177e-0035L,
188 +7.361337767588456524131930836633932195088e-0035L,
189 -6.472995147913347230035214575612170525266e-0035L,
190 +8.587474417953698694278798062295229624207e-0035L,
191 +2.371815422825174835691651228302690977951e-0035L,
192 -3.026891682096118773004597373421900314256e-0037L,
193 };
194 #endif
195
196 static const long double one = 1.0L,
197 two = 2.0L,
198 ln2_64 = 1.083042469624914545964425189778400898568e-2L,
199 invln2_32 = 4.616624130844682903551758979206054839765e+1L;
200
201 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
202 static const long double t1 = 1.666666666666666666666666666660876387437e-1L,
203 t2 = -2.777777777777777777777707812093173478756e-3L,
204 t3 = 6.613756613756613482074280932874221202424e-5L,
205 t4 = -1.653439153392139954169609822742235851120e-6L,
206 t5 = 4.175314851769539751387852116610973796053e-8L;
207
208 long double
209 __k_cexpl(long double x, int *n)
210 {
211 int hx, ix, j, k;
212 long double t, r;
213
214 *n = 0;
215 hx = HI_XWORD(x);
216 ix = hx & 0x7fffffff;
217
218 if (hx >= 0x7fff0000)
219 return (x + x); /* NaN of +inf */
220
221 if (((unsigned)hx) >= 0xffff0000)
222 return (-one / x); /* NaN or -inf */
223
224 if (ix < 0x3fc30000)
225 return (one + x); /* |x|<2^-60 */
226
227 if (hx > 0) {
228 if (hx > 0x401086a0) { /* x > 200000 */
229 *n = 200000;
230 return (one);
231 }
232
233 k = (int)(invln2_32 * (x + ln2_64));
234 } else {
235 if (ix > 0x401086a0) { /* x < -200000 */
236 *n = -200000;
237 return (one);
238 }
239
240 k = (int)(invln2_32 * (x - ln2_64));
241 }
242
243 j = k & 0x1f;
244 *n = k >> 5;
245 t = (long double)k;
246 x = (x - t * ln2_32hi) - t * ln2_32lo;
247 t = x * x;
248 r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
249 x = exp2_32_hi[j] - ((exp2_32_hi[j] * (x + x)) / r - exp2_32_lo[j]);
250 k >>= 5;
251
252 if (k > 240) {
253 XFSCALE(x, 240);
254 *n -= 240;
255 } else if (k > 0) {
256 XFSCALE(x, k);
257 *n = 0;
258 }
259
260 return (x);
261 }
|