Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/C/j1.c
+++ new/usr/src/lib/libm/common/C/j1.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: j1(x),y1(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 __j1 = j1
40 41 #pragma weak __y1 = y1
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-20,
52 -one = 1.0,
53 -invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
54 -tpi = 0.636619772367581343075535053490057448;
48 +#define GENERIC double
49 +
50 +static const GENERIC zero = 0.0,
51 + small = 1.0e-5,
52 + tiny = 1.0e-20,
53 + one = 1.0,
54 + invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
55 + tpi = 0.636619772367581343075535053490057448;
56 +
57 +static GENERIC pone(GENERIC);
58 +static GENERIC qone(GENERIC);
55 59
56 -static GENERIC pone(GENERIC), qone(GENERIC);
57 60 static const GENERIC r0[4] = {
58 61 -6.250000000000002203053200981413218949548e-0002,
59 62 1.600998455640072901321605101981501263762e-0003,
60 63 -1.963888815948313758552511884390162864930e-0005,
61 64 8.263917341093549759781339713418201620998e-0008,
62 65 };
66 +
63 67 static const GENERIC s0[7] = {
64 68 1.0e0,
65 69 1.605069137643004242395356851797873766927e-0002,
66 70 1.149454623251299996428500249509098499383e-0004,
67 71 3.849701673735260970379681807910852327825e-0007,
68 72 };
73 +
69 74 static const GENERIC r1[12] = {
70 75 4.999999999999999995517408894340485471724e-0001,
71 76 -6.003825028120475684835384519945468075423e-0002,
72 77 2.301719899263321828388344461995355419832e-0003,
73 78 -4.208494869238892934859525221654040304068e-0005,
74 79 4.377745135188837783031540029700282443388e-0007,
75 80 -2.854106755678624335145364226735677754179e-0009,
76 81 1.234002865443952024332943901323798413689e-0011,
77 82 -3.645498437039791058951273508838177134310e-0014,
78 83 7.404320596071797459925377103787837414422e-0017,
79 84 -1.009457448277522275262808398517024439084e-0019,
80 85 8.520158355824819796968771418801019930585e-0023,
81 86 -3.458159926081163274483854614601091361424e-0026,
82 87 };
88 +
83 89 static const GENERIC s1[5] = {
84 90 1.0e0,
85 91 4.923499437590484879081138588998986303306e-0003,
86 92 1.054389489212184156499666953501976688452e-0005,
87 93 1.180768373106166527048240364872043816050e-0008,
88 94 5.942665743476099355323245707680648588540e-0012,
89 95 };
90 96
91 97 GENERIC
92 -j1(GENERIC x) {
98 +j1(GENERIC x)
99 +{
93 100 GENERIC z, d, s, c, ss, cc, r;
94 101 int i, sgn;
95 102
96 103 if (!finite(x))
97 - return (one/x);
104 + return (one / x);
105 +
98 106 sgn = signbit(x);
99 107 x = fabs(x);
108 +
100 109 if (x > 8.00) {
101 110 s = sin(x);
102 111 c = cos(x);
103 - /*
104 - * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
105 - * where x0 = x-3pi/4
106 - * Better formula:
107 - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
108 - * = 1/sqrt(2) * (sin(x) - cos(x))
109 - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
110 - * = -1/sqrt(2) * (cos(x) + sin(x))
111 - * To avoid cancellation, use
112 - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
113 - * to compute the worse one.
114 - */
112 +
113 + /* BEGIN CSTYLED */
114 + /*
115 + * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
116 + * where x0 = x-3pi/4
117 + * Better formula:
118 + * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
119 + * = 1/sqrt(2) * (sin(x) - cos(x))
120 + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
121 + * = -1/sqrt(2) * (cos(x) + sin(x))
122 + * To avoid cancellation, use
123 + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
124 + * to compute the worse one.
125 + */
126 + /* END CSTYLED */
115 127 if (x > 8.9e307) { /* x+x may overflow */
116 - ss = -s-c;
117 - cc = s-c;
128 + ss = -s - c;
129 + cc = s - c;
118 130 } else if (signbit(s) != signbit(c)) {
119 131 cc = s - c;
120 - ss = cos(x+x)/cc;
132 + ss = cos(x + x) / cc;
121 133 } else {
122 - ss = -s-c;
123 - cc = cos(x+x)/ss;
134 + ss = -s - c;
135 + cc = cos(x + x) / ss;
124 136 }
125 - /*
126 - * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
127 - * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
128 - */
137 +
138 + /*
139 + * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
140 + * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
141 + */
129 142 if (x > 1.0e40)
130 - d = (invsqrtpi*cc)/sqrt(x);
143 + d = (invsqrtpi * cc) / sqrt(x);
131 144 else
132 - d = invsqrtpi*(pone(x)*cc-qone(x)*ss)/sqrt(x);
145 + d = invsqrtpi * (pone(x) * cc - qone(x) * ss) / sqrt(x);
133 146
134 147 if (x > X_TLOSS) {
135 - if (sgn != 0) { d = -d; x = -x; }
148 + if (sgn != 0) {
149 + d = -d;
150 + x = -x;
151 + }
152 +
136 153 return (_SVID_libm_err(x, d, 36));
137 - } else
138 - if (sgn == 0)
139 - return (d);
140 - else
141 - return (-d);
154 + } else if (sgn == 0) {
155 + return (d);
156 + } else {
157 + return (-d);
158 + }
142 159 }
160 +
143 161 if (x <= small) {
144 162 if (x <= tiny)
145 - d = 0.5*x;
163 + d = 0.5 * x;
146 164 else
147 - d = x*(0.5-x*x*0.125);
165 + d = x * (0.5 - x * x * 0.125);
166 +
148 167 if (sgn == 0)
149 168 return (d);
150 169 else
151 170 return (-d);
152 171 }
153 - z = x*x;
172 +
173 + z = x * x;
174 +
154 175 if (x < 1.28) {
155 - r = r0[3];
156 - s = s0[3];
157 - for (i = 2; i >= 0; i--) {
158 - r = r*z + r0[i];
159 - s = s*z + s0[i];
160 - }
161 - d = x*0.5+x*(z*(r/s));
176 + r = r0[3];
177 + s = s0[3];
178 +
179 + for (i = 2; i >= 0; i--) {
180 + r = r * z + r0[i];
181 + s = s * z + s0[i];
182 + }
183 +
184 + d = x * 0.5 + x * (z * (r / s));
162 185 } else {
163 - r = r1[11];
164 - for (i = 10; i >= 0; i--) r = r*z + r1[i];
165 - s = s1[0]+z*(s1[1]+z*(s1[2]+z*(s1[3]+z*s1[4])));
166 - d = x*(r/s);
186 + r = r1[11];
187 +
188 + for (i = 10; i >= 0; i--)
189 + r = r * z + r1[i];
190 +
191 + s = s1[0] + z * (s1[1] + z * (s1[2] + z * (s1[3] + z * s1[4])));
192 + d = x * (r / s);
167 193 }
194 +
168 195 if (sgn == 0)
169 196 return (d);
170 197 else
171 198 return (-d);
172 199 }
173 200
174 201 static const GENERIC u0[4] = {
175 202 -1.960570906462389461018983259589655961560e-0001,
176 203 4.931824118350661953459180060007970291139e-0002,
177 204 -1.626975871565393656845930125424683008677e-0003,
178 205 1.359657517926394132692884168082224258360e-0005,
179 206 };
207 +
180 208 static const GENERIC v0[5] = {
181 209 1.0e0,
182 210 2.565807214838390835108224713630901653793e-0002,
183 211 3.374175208978404268650522752520906231508e-0004,
184 212 2.840368571306070719539936935220728843177e-0006,
185 213 1.396387402048998277638900944415752207592e-0008,
186 214 };
215 +
187 216 static const GENERIC u1[12] = {
188 217 -1.960570906462389473336339614647555351626e-0001,
189 218 5.336268030335074494231369159933012844735e-0002,
190 219 -2.684137504382748094149184541866332033280e-0003,
191 220 5.737671618979185736981543498580051903060e-0005,
192 221 -6.642696350686335339171171785557663224892e-0007,
193 222 4.692417922568160354012347591960362101664e-0009,
194 223 -2.161728635907789319335231338621412258355e-0011,
195 224 6.727353419738316107197644431844194668702e-0014,
196 225 -1.427502986803861372125234355906790573422e-0016,
197 226 2.020392498726806769468143219616642940371e-0019,
198 227 -1.761371948595104156753045457888272716340e-0022,
199 228 7.352828391941157905175042420249225115816e-0026,
200 229 };
230 +
201 231 static const GENERIC v1[5] = {
202 232 1.0e0,
203 233 5.029187436727947764916247076102283399442e-0003,
204 234 1.102693095808242775074856548927801750627e-0005,
205 235 1.268035774543174837829534603830227216291e-0008,
206 236 6.579416271766610825192542295821308730206e-0012,
207 237 };
208 238
209 -
210 239 GENERIC
211 -y1(GENERIC x) {
240 +y1(GENERIC x)
241 +{
212 242 GENERIC z, d, s, c, ss, cc, u, v;
213 243 int i;
214 244
215 245 if (isnan(x))
216 - return (x*x); /* + -> * for Cheetah */
246 + return (x * x); /* + -> * for Cheetah */
247 +
217 248 if (x <= zero) {
218 249 if (x == zero)
219 - /* return -one/zero; */
220 - return (_SVID_libm_err(x, x, 10));
250 + /* return -one/zero; */
251 + return (_SVID_libm_err(x, x, 10));
221 252 else
222 - /* return zero/zero; */
223 - return (_SVID_libm_err(x, x, 11));
253 + /* return zero/zero; */
254 + return (_SVID_libm_err(x, x, 11));
224 255 }
256 +
225 257 if (x > 8.0) {
226 258 if (!finite(x))
227 259 return (zero);
260 +
228 261 s = sin(x);
229 262 c = cos(x);
230 - /*
231 - * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
232 - * where x0 = x-3pi/4
233 - * Better formula:
234 - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
235 - * = 1/sqrt(2) * (sin(x) - cos(x))
236 - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
237 - * = -1/sqrt(2) * (cos(x) + sin(x))
238 - * To avoid cancellation, use
239 - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
240 - * to compute the worse one.
241 - */
263 +
264 + /* BEGIN CSTYLED */
265 + /*
266 + * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x0)-q1(x)*sin(x0))
267 + * where x0 = x-3pi/4
268 + * Better formula:
269 + * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
270 + * = 1/sqrt(2) * (sin(x) - cos(x))
271 + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
272 + * = -1/sqrt(2) * (cos(x) + sin(x))
273 + * To avoid cancellation, use
274 + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
275 + * to compute the worse one.
276 + */
277 + /* END CSTYLED */
242 278 if (x > 8.9e307) { /* x+x may overflow */
243 - ss = -s-c;
244 - cc = s-c;
279 + ss = -s - c;
280 + cc = s - c;
245 281 } else if (signbit(s) != signbit(c)) {
246 282 cc = s - c;
247 - ss = cos(x+x)/cc;
283 + ss = cos(x + x) / cc;
248 284 } else {
249 - ss = -s-c;
250 - cc = cos(x+x)/ss;
285 + ss = -s - c;
286 + cc = cos(x + x) / ss;
251 287 }
252 - /*
253 - * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
254 - * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
255 - */
288 +
289 + /*
290 + * j1(x) = 1/sqrt(pi*x) * (P(1,x)*cc - Q(1,x)*ss)
291 + * y1(x) = 1/sqrt(pi*x) * (P(1,x)*ss + Q(1,x)*cc)
292 + */
256 293 if (x > 1.0e91)
257 - d = (invsqrtpi*ss)/sqrt(x);
294 + d = (invsqrtpi * ss) / sqrt(x);
258 295 else
259 - d = invsqrtpi*(pone(x)*ss+qone(x)*cc)/sqrt(x);
296 + d = invsqrtpi * (pone(x) * ss + qone(x) * cc) / sqrt(x);
260 297
261 298 if (x > X_TLOSS)
262 299 return (_SVID_libm_err(x, d, 37));
263 300 else
264 301 return (d);
265 302 }
266 - if (x <= tiny) {
267 - return (-tpi/x);
268 - }
269 - z = x*x;
303 +
304 + if (x <= tiny)
305 + return (-tpi / x);
306 +
307 + z = x * x;
308 +
270 309 if (x < 1.28) {
271 - u = u0[3]; v = v0[3]+z*v0[4];
272 - for (i = 2; i >= 0; i--) {
273 - u = u*z + u0[i];
274 - v = v*z + v0[i];
275 - }
310 + u = u0[3];
311 + v = v0[3] + z * v0[4];
312 +
313 + for (i = 2; i >= 0; i--) {
314 + u = u * z + u0[i];
315 + v = v * z + v0[i];
316 + }
276 317 } else {
277 - for (u = u1[11], i = 10; i >= 0; i--) u = u*z+u1[i];
278 - v = v1[0]+z*(v1[1]+z*(v1[2]+z*(v1[3]+z*v1[4])));
318 + for (u = u1[11], i = 10; i >= 0; i--)
319 + u = u * z + u1[i];
320 +
321 + v = v1[0] + z * (v1[1] + z * (v1[2] + z * (v1[3] + z * v1[4])));
279 322 }
280 - return (x*(u/v) + tpi*(j1(x)*log(x)-one/x));
323 +
324 + return (x * (u / v) + tpi * (j1(x) * log(x) - one / x));
281 325 }
282 326
283 327 static const GENERIC pr0[6] = {
284 - -.4435757816794127857114720794e7,
285 - -.9942246505077641195658377899e7,
286 - -.6603373248364939109255245434e7,
287 - -.1523529351181137383255105722e7,
288 - -.1098240554345934672737413139e6,
289 - -.1611616644324610116477412898e4,
328 + -.4435757816794127857114720794e7, -.9942246505077641195658377899e7,
329 + -.6603373248364939109255245434e7, -.1523529351181137383255105722e7,
330 + -.1098240554345934672737413139e6, -.1611616644324610116477412898e4,
290 331 };
332 +
291 333 static const GENERIC ps0[6] = {
292 - -.4435757816794127856828016962e7,
293 - -.9934124389934585658967556309e7,
294 - -.6585339479723087072826915069e7,
295 - -.1511809506634160881644546358e7,
296 - -.1072638599110382011903063867e6,
297 - -.1455009440190496182453565068e4,
334 + -.4435757816794127856828016962e7, -.9934124389934585658967556309e7,
335 + -.6585339479723087072826915069e7, -.1511809506634160881644546358e7,
336 + -.1072638599110382011903063867e6, -.1455009440190496182453565068e4,
298 337 };
299 -static const GENERIC huge = 1.0e10;
300 338
339 +static const GENERIC huge = 1.0e10;
301 340 static GENERIC
302 -pone(GENERIC x) {
341 +pone(GENERIC x)
342 +{
303 343 GENERIC s, r, t, z;
304 344 int i;
305 - /* assume x > 8 */
345 +
346 + /* assume x > 8 */
306 347 if (x > huge)
307 348 return (one);
308 349
309 - t = 8.0/x; z = t*t;
310 - r = pr0[5]; s = ps0[5]+z;
350 + t = 8.0 / x;
351 + z = t * t;
352 + r = pr0[5];
353 + s = ps0[5] + z;
354 +
311 355 for (i = 4; i >= 0; i--) {
312 - r = z*r + pr0[i];
313 - s = z*s + ps0[i];
356 + r = z * r + pr0[i];
357 + s = z * s + ps0[i];
314 358 }
315 - return (r/s);
316 -}
317 359
360 + return (r / s);
361 +}
318 362
319 363 static const GENERIC qr0[6] = {
320 - 0.3322091340985722351859704442e5,
321 - 0.8514516067533570196555001171e5,
322 - 0.6617883658127083517939992166e5,
323 - 0.1849426287322386679652009819e5,
324 - 0.1706375429020768002061283546e4,
325 - 0.3526513384663603218592175580e2,
364 + 0.3322091340985722351859704442e5, 0.8514516067533570196555001171e5,
365 + 0.6617883658127083517939992166e5, 0.1849426287322386679652009819e5,
366 + 0.1706375429020768002061283546e4, 0.3526513384663603218592175580e2,
326 367 };
368 +
327 369 static const GENERIC qs0[6] = {
328 - 0.7087128194102874357377502472e6,
329 - 0.1819458042243997298924553839e7,
330 - 0.1419460669603720892855755253e7,
331 - 0.4002944358226697511708610813e6,
332 - 0.3789022974577220264142952256e5,
333 - 0.8638367769604990967475517183e3,
370 + 0.7087128194102874357377502472e6, 0.1819458042243997298924553839e7,
371 + 0.1419460669603720892855755253e7, 0.4002944358226697511708610813e6,
372 + 0.3789022974577220264142952256e5, 0.8638367769604990967475517183e3,
334 373 };
335 374
336 375 static GENERIC
337 -qone(GENERIC x) {
376 +qone(GENERIC x)
377 +{
338 378 GENERIC s, r, t, z;
339 379 int i;
380 +
340 381 if (x > huge)
341 - return (0.375/x);
382 + return (0.375 / x);
383 +
384 + t = 8.0 / x;
385 + z = t * t;
386 + /* assume x > 8 */
387 + r = qr0[5];
388 + s = qs0[5] + z;
342 389
343 - t = 8.0/x; z = t*t;
344 - /* assume x > 8 */
345 - r = qr0[5]; s = qs0[5]+z;
346 390 for (i = 4; i >= 0; i--) {
347 - r = z*r + qr0[i];
348 - s = z*s + qs0[i];
391 + r = z * r + qr0[i];
392 + s = z * s + qs0[i];
349 393 }
350 - return (t*(r/s));
394 +
395 + return (t * (r / s));
351 396 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX