Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/C/j0.c
+++ new/usr/src/lib/libm/common/C/j0.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
|
↓ open down ↓ |
14 lines elided |
↑ open up ↑ |
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 +
25 26 /*
26 27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 28 * Use is subject to license terms.
28 29 */
29 30
30 31 /*
31 32 * Floating point Bessel's function of the first and second kinds
32 33 * of order zero: j0(x),y0(x);
33 34 *
34 35 * Special cases:
35 36 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
36 37 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
|
↓ open down ↓ |
2 lines elided |
↑ open up ↑ |
37 38 */
38 39
39 40 #pragma weak __j0 = j0
40 41 #pragma weak __y0 = y0
41 42
42 43 #include "libm.h"
43 44 #include "libm_protos.h"
44 45 #include <math.h>
45 46 #include <values.h>
46 47
47 -#define GENERIC double
48 -static const GENERIC
49 -zero = 0.0,
50 -small = 1.0e-5,
51 -tiny = 1.0e-18,
52 -one = 1.0,
53 -eight = 8.0,
54 -invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
55 -tpi = 0.636619772367581343075535053490057448;
48 +#define GENERIC double
49 +
50 +static const GENERIC zero = 0.0,
51 + small = 1.0e-5,
52 + tiny = 1.0e-18,
53 + one = 1.0,
54 + eight = 8.0,
55 + invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
56 + tpi = 0.636619772367581343075535053490057448;
57 +
58 +static GENERIC pzero(GENERIC);
59 +static GENERIC qzero(GENERIC);
56 60
57 -static GENERIC pzero(GENERIC), qzero(GENERIC);
58 -static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
61 +static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
59 62 -2.500000000000003622131880894830476755537e-0001,
60 63 1.095597547334830263234433855932375353303e-0002,
61 64 -1.819734750463320921799187258987098087697e-0004,
62 65 9.977001946806131657544212501069893930846e-0007,
63 66 };
64 -static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
67 +
68 +static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
65 69 1.0,
66 70 1.867609810662950169966782360588199673741e-0002,
67 71 1.590389206181565490878430827706972074208e-0004,
68 72 6.520867386742583632375520147714499522721e-0007,
69 73 };
70 -static const GENERIC r1[9] = { /* [1.28,8] */
74 +
75 +static const GENERIC r1[9] = { /* [1.28,8] */
71 76 9.999999999999999942156495584397047660949e-0001,
72 77 -2.389887722731319130476839836908143731281e-0001,
73 78 1.293359476138939027791270393439493640570e-0002,
74 79 -2.770985642343140122168852400228563364082e-0004,
75 80 2.905241575772067678086738389169625218912e-0006,
76 81 -1.636846356264052597969042009265043251279e-0008,
77 82 5.072306160724884775085431059052611737827e-0011,
78 83 -8.187060730684066824228914775146536139112e-0014,
79 84 5.422219326959949863954297860723723423842e-0017,
80 85 };
81 -static const GENERIC s1[9] = { /* [1.28,8] */
86 +
87 +static const GENERIC s1[9] = { /* [1.28,8] */
82 88 1.0,
83 89 1.101122772686807702762104741932076228349e-0002,
84 90 6.140169310641649223411427764669143978228e-0005,
85 91 2.292035877515152097976946119293215705250e-0007,
86 92 6.356910426504644334558832036362219583789e-0010,
87 93 1.366626326900219555045096999553948891401e-0012,
88 94 2.280399586866739522891837985560481180088e-0015,
89 95 2.801559820648939665270492520004836611187e-0018,
90 96 2.073101088320349159764410261466350732968e-0021,
91 97 };
92 98
93 99 GENERIC
94 -j0(GENERIC x) {
100 +j0(GENERIC x)
101 +{
95 102 GENERIC z, s, c, ss, cc, r, u, v, ox;
96 103 int i;
97 104
98 105 if (isnan(x))
99 - return (x*x); /* + -> * for Cheetah */
106 + return (x * x); /* + -> * for Cheetah */
107 +
100 108 ox = x;
101 109 x = fabs(x);
110 +
102 111 if (x > 8.0) {
103 112 if (!finite(x))
104 113 return (zero);
114 +
105 115 s = sin(x);
106 116 c = cos(x);
107 - /*
108 - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
109 - * where x0 = x-pi/4
110 - * Better formula:
111 - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
112 - * = 1/sqrt(2) * (cos(x) + sin(x))
113 - * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
114 - * = 1/sqrt(2) * (sin(x) - cos(x))
115 - * To avoid cancellation, use
116 - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
117 - * to compute the worse one.
118 - */
117 +
118 + /* BEGIN CSTYLED */
119 + /*
120 + * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
121 + * where x0 = x-pi/4
122 + * Better formula:
123 + * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
124 + * = 1/sqrt(2) * (cos(x) + sin(x))
125 + * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
126 + * = 1/sqrt(2) * (sin(x) - cos(x))
127 + * To avoid cancellation, use
128 + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
129 + * to compute the worse one.
130 + */
131 + /* END CSTYLED */
119 132 if (x > 8.9e307) { /* x+x may overflow */
120 - ss = s-c;
121 - cc = s+c;
133 + ss = s - c;
134 + cc = s + c;
122 135 } else if (signbit(s) != signbit(c)) {
123 136 ss = s - c;
124 - cc = -cos(x+x)/ss;
137 + cc = -cos(x + x) / ss;
125 138 } else {
126 139 cc = s + c;
127 - ss = -cos(x+x)/cc;
140 + ss = -cos(x + x) / cc;
128 141 }
129 - /*
130 - * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
131 - * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
132 - */
133 - if (x > 1.0e40) z = (invsqrtpi*cc)/sqrt(x);
134 - else {
135 - u = pzero(x); v = qzero(x);
136 - z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
142 +
143 + /*
144 + * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
145 + * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
146 + */
147 + if (x > 1.0e40) {
148 + z = (invsqrtpi * cc) / sqrt(x);
149 + } else {
150 + u = pzero(x);
151 + v = qzero(x);
152 + z = invsqrtpi * (u * cc - v * ss) / sqrt(x);
137 153 }
138 - /* force to pass SVR4 even the result is wrong (sign) */
154 +
155 + /* force to pass SVR4 even the result is wrong (sign) */
139 156 if (x > X_TLOSS)
140 - return (_SVID_libm_err(ox, z, 34));
157 + return (_SVID_libm_err(ox, z, 34));
141 158 else
142 - return (z);
159 + return (z);
143 160 }
161 +
144 162 if (x <= small) {
145 - if (x <= tiny)
146 - return (one-x);
147 - else
148 - return (one-x*x*0.25);
163 + if (x <= tiny)
164 + return (one - x);
165 + else
166 + return (one - x * x * 0.25);
149 167 }
150 - z = x*x;
168 +
169 + z = x * x;
170 +
151 171 if (x <= 1.28) {
152 - r = r0[0]+z*(r0[1]+z*(r0[2]+z*r0[3]));
153 - s = s0[0]+z*(s0[1]+z*(s0[2]+z*s0[3]));
154 - return (one + z*(r/s));
172 + r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
173 + s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
174 + return (one + z * (r / s));
155 175 } else {
156 - for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
157 - r = r*z + r1[i];
158 - s = s*z + s1[i];
159 - }
160 - return (r/s);
176 + for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
177 + r = r * z + r1[i];
178 + s = s * z + s1[i];
179 + }
180 +
181 + return (r / s);
161 182 }
162 183 }
163 184
164 185 static const GENERIC u0[13] = {
165 186 -7.380429510868722526754723020704317641941e-0002,
166 187 1.772607102684869924301459663049874294814e-0001,
167 188 -1.524370666542713828604078090970799356306e-0002,
168 189 4.650819100693891757143771557629924591915e-0004,
169 190 -7.125768872339528975036316108718239946022e-0006,
170 191 6.411017001656104598327565004771515257146e-0008,
171 192 -3.694275157433032553021246812379258781665e-0010,
172 193 1.434364544206266624252820889648445263842e-0012,
173 194 -3.852064731859936455895036286874139896861e-0015,
174 195 7.182052899726138381739945881914874579696e-0018,
175 196 -9.060556574619677567323741194079797987200e-0021,
176 197 7.124435467408860515265552217131230511455e-0024,
177 198 -2.709726774636397615328813121715432044771e-0027,
178 199 };
200 +
179 201 static const GENERIC v0[5] = {
180 202 1.0,
181 203 4.678678931512549002587702477349214886475e-0003,
182 204 9.486828955529948534822800829497565178985e-0006,
183 205 1.001495929158861646659010844136682454906e-0008,
184 206 4.725338116256021660204443235685358593611e-0012,
185 207 };
186 208
187 209 GENERIC
188 -y0(GENERIC x) {
210 +y0(GENERIC x)
211 +{
189 212 GENERIC z, /* d, */ s, c, ss, cc, u, v;
190 213 int i;
191 214
192 215 if (isnan(x))
193 - return (x*x); /* + -> * for Cheetah */
216 + return (x * x); /* + -> * for Cheetah */
217 +
194 218 if (x <= zero) {
195 219 if (x == zero)
196 - /* d= -one/(x-x); */
197 - return (_SVID_libm_err(x, x, 8));
220 + /* d= -one/(x-x); */
221 + return (_SVID_libm_err(x, x, 8));
198 222 else
199 - /* d = zero/(x-x); */
200 - return (_SVID_libm_err(x, x, 9));
223 + /* d = zero/(x-x); */
224 + return (_SVID_libm_err(x, x, 9));
201 225 }
226 +
202 227 if (x > 8.0) {
203 228 if (!finite(x))
204 229 return (zero);
230 +
205 231 s = sin(x);
206 232 c = cos(x);
207 - /*
208 - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
209 - * where x0 = x-pi/4
210 - * Better formula:
211 - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
212 - * = 1/sqrt(2) * (cos(x) + sin(x))
213 - * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
214 - * = 1/sqrt(2) * (sin(x) - cos(x))
215 - * To avoid cancellation, use
216 - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
217 - * to compute the worse one.
218 - */
233 +
234 + /* BEGIN CSTYLED */
235 + /*
236 + * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
237 + * where x0 = x-pi/4
238 + * Better formula:
239 + * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
240 + * = 1/sqrt(2) * (cos(x) + sin(x))
241 + * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
242 + * = 1/sqrt(2) * (sin(x) - cos(x))
243 + * To avoid cancellation, use
244 + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
245 + * to compute the worse one.
246 + */
247 + /* END CSTYLED */
219 248 if (x > 8.9e307) { /* x+x may overflow */
220 - ss = s-c;
221 - cc = s+c;
249 + ss = s - c;
250 + cc = s + c;
222 251 } else if (signbit(s) != signbit(c)) {
223 252 ss = s - c;
224 - cc = -cos(x+x)/ss;
253 + cc = -cos(x + x) / ss;
225 254 } else {
226 255 cc = s + c;
227 - ss = -cos(x+x)/cc;
256 + ss = -cos(x + x) / cc;
228 257 }
229 - /*
230 - * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
231 - * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
232 - */
258 +
259 + /*
260 + * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
261 + * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
262 + */
233 263 if (x > 1.0e40)
234 - z = (invsqrtpi*ss)/sqrt(x);
264 + z = (invsqrtpi * ss) / sqrt(x);
235 265 else
236 - z = invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x);
266 + z = invsqrtpi * (pzero(x) * ss + qzero(x) * cc) /
267 + sqrt(x);
268 +
237 269 if (x > X_TLOSS)
238 - return (_SVID_libm_err(x, z, 35));
270 + return (_SVID_libm_err(x, z, 35));
239 271 else
240 - return (z);
241 -
242 - }
243 - if (x <= tiny) {
244 - return (u0[0] + tpi*log(x));
272 + return (z);
245 273 }
246 - z = x*x;
247 - for (u = u0[12], i = 11; i >= 0; i--) u = u*z + u0[i];
248 - v = v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])));
249 - return (u/v + tpi*(j0(x)*log(x)));
274 +
275 + if (x <= tiny)
276 + return (u0[0] + tpi * log(x));
277 +
278 + z = x * x;
279 +
280 + for (u = u0[12], i = 11; i >= 0; i--)
281 + u = u * z + u0[i];
282 +
283 + v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
284 + return (u / v + tpi * (j0(x) * log(x)));
250 285 }
251 286
252 -static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
253 - .4861344183386052721391238447e5,
254 - .1377662549407112278133438945e6,
255 - .1222466364088289731869114004e6,
256 - .4107070084315176135583353374e5,
257 - .5026073801860637125889039915e4,
258 - .1783193659125479654541542419e3,
287 +static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
288 + .4861344183386052721391238447e5, .1377662549407112278133438945e6,
289 + .1222466364088289731869114004e6, .4107070084315176135583353374e5,
290 + .5026073801860637125889039915e4, .1783193659125479654541542419e3,
259 291 .88010344055383421691677564e0,
260 292 };
261 -static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
262 - .4861344183386052721414037058e5,
263 - .1378196632630384670477582699e6,
264 - .1223967185341006542748936787e6,
265 - .4120150243795353639995862617e5,
266 - .5068271181053546392490184353e4,
267 - .1829817905472769960535671664e3,
293 +
294 +static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
295 + .4861344183386052721414037058e5, .1378196632630384670477582699e6,
296 + .1223967185341006542748936787e6, .4120150243795353639995862617e5,
297 + .5068271181053546392490184353e4, .1829817905472769960535671664e3,
268 298 1.0,
269 299 };
270 -static const GENERIC huge = 1.0e10;
271 300
301 +static const GENERIC huge = 1.0e10;
272 302 static GENERIC
273 -pzero(GENERIC x) {
303 +pzero(GENERIC x)
304 +{
274 305 GENERIC s, r, t, z;
275 306 int i;
307 +
276 308 if (x > huge)
277 309 return (one);
278 - t = eight/x; z = t*t;
279 - r = pr[5]+z*pr[6];
280 - s = ps[5]+z;
310 +
311 + t = eight / x;
312 + z = t * t;
313 + r = pr[5] + z * pr[6];
314 + s = ps[5] + z;
315 +
281 316 for (i = 4; i >= 0; i--) {
282 - r = r*z + pr[i];
283 - s = s*z + ps[i];
317 + r = r * z + pr[i];
318 + s = s * z + ps[i];
284 319 }
285 - return (r/s);
320 +
321 + return (r / s);
286 322 }
287 323
288 -static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
289 - -.1731210995701068539185611951e3,
290 - -.5522559165936166961235240613e3,
291 - -.5604935606637346590614529613e3,
292 - -.2200430300226009379477365011e3,
293 - -.323869355375648849771296746e2,
294 - -.14294979207907956223499258e1,
324 +static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
325 + -.1731210995701068539185611951e3, -.5522559165936166961235240613e3,
326 + -.5604935606637346590614529613e3, -.2200430300226009379477365011e3,
327 + -.323869355375648849771296746e2, -.14294979207907956223499258e1,
295 328 -.834690374102384988158918e-2,
296 329 };
297 -static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
298 - .1107975037248683865326709645e5,
299 - .3544581680627082674651471873e5,
300 - .3619118937918394132179019059e5,
301 - .1439895563565398007471485822e5,
302 - .2190277023344363955930226234e4,
303 - .106695157020407986137501682e3,
330 +
331 +static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
332 + .1107975037248683865326709645e5, .3544581680627082674651471873e5,
333 + .3619118937918394132179019059e5, .1439895563565398007471485822e5,
334 + .2190277023344363955930226234e4, .106695157020407986137501682e3,
304 335 1.0,
305 336 };
306 337
307 338 static GENERIC
308 -qzero(GENERIC x) {
339 +qzero(GENERIC x)
340 +{
309 341 GENERIC s, r, t, z;
310 342 int i;
343 +
311 344 if (x > huge)
312 - return (-0.125/x);
313 - t = eight/x; z = t*t;
314 - r = qr[5]+z*qr[6];
315 - s = qs[5]+z;
345 + return (-0.125 / x);
346 +
347 + t = eight / x;
348 + z = t * t;
349 + r = qr[5] + z * qr[6];
350 + s = qs[5] + z;
351 +
316 352 for (i = 4; i >= 0; i--) {
317 - r = r*z + qr[i];
318 - s = s*z + qs[i];
353 + r = r * z + qr[i];
354 + s = s * z + qs[i];
319 355 }
320 - return (t*(r/s));
356 +
357 + return (t * (r / s));
321 358 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX