1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
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 * Floating point Bessel's function of the first and second kinds
33 * of order zero: j0(x),y0(x);
34 *
35 * Special cases:
36 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
37 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
38 */
39
40 #pragma weak __j0 = j0
41 #pragma weak __y0 = y0
42
43 #include "libm.h"
44 #include "libm_protos.h"
45 #include <math.h>
46 #include <values.h>
47
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);
60
61 static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
62 -2.500000000000003622131880894830476755537e-0001,
63 1.095597547334830263234433855932375353303e-0002,
64 -1.819734750463320921799187258987098087697e-0004,
65 9.977001946806131657544212501069893930846e-0007,
66 };
67
68 static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
69 1.0,
70 1.867609810662950169966782360588199673741e-0002,
71 1.590389206181565490878430827706972074208e-0004,
72 6.520867386742583632375520147714499522721e-0007,
73 };
74
75 static const GENERIC r1[9] = { /* [1.28,8] */
76 9.999999999999999942156495584397047660949e-0001,
77 -2.389887722731319130476839836908143731281e-0001,
78 1.293359476138939027791270393439493640570e-0002,
79 -2.770985642343140122168852400228563364082e-0004,
80 2.905241575772067678086738389169625218912e-0006,
81 -1.636846356264052597969042009265043251279e-0008,
82 5.072306160724884775085431059052611737827e-0011,
83 -8.187060730684066824228914775146536139112e-0014,
84 5.422219326959949863954297860723723423842e-0017,
85 };
86
87 static const GENERIC s1[9] = { /* [1.28,8] */
88 1.0,
89 1.101122772686807702762104741932076228349e-0002,
90 6.140169310641649223411427764669143978228e-0005,
91 2.292035877515152097976946119293215705250e-0007,
92 6.356910426504644334558832036362219583789e-0010,
93 1.366626326900219555045096999553948891401e-0012,
94 2.280399586866739522891837985560481180088e-0015,
95 2.801559820648939665270492520004836611187e-0018,
96 2.073101088320349159764410261466350732968e-0021,
97 };
98
99 GENERIC
100 j0(GENERIC x)
101 {
102 GENERIC z, s, c, ss, cc, r, u, v, ox;
103 int i;
104
105 if (isnan(x))
106 return (x * x); /* + -> * for Cheetah */
107
108 ox = x;
109 x = fabs(x);
110
111 if (x > 8.0) {
112 if (!finite(x))
113 return (zero);
114
115 s = sin(x);
116 c = cos(x);
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 */
132 if (x > 8.9e307) { /* x+x may overflow */
133 ss = s - c;
134 cc = s + c;
135 } else if (signbit(s) != signbit(c)) {
136 ss = s - c;
137 cc = -cos(x + x) / ss;
138 } else {
139 cc = s + c;
140 ss = -cos(x + x) / cc;
141 }
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);
153 }
154
155 /* force to pass SVR4 even the result is wrong (sign) */
156 if (x > X_TLOSS)
157 return (_SVID_libm_err(ox, z, 34));
158 else
159 return (z);
160 }
161
162 if (x <= small) {
163 if (x <= tiny)
164 return (one - x);
165 else
166 return (one - x * x * 0.25);
167 }
168
169 z = x * x;
170
171 if (x <= 1.28) {
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));
175 } else {
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);
182 }
183 }
184
185 static const GENERIC u0[13] = {
186 -7.380429510868722526754723020704317641941e-0002,
187 1.772607102684869924301459663049874294814e-0001,
188 -1.524370666542713828604078090970799356306e-0002,
189 4.650819100693891757143771557629924591915e-0004,
190 -7.125768872339528975036316108718239946022e-0006,
191 6.411017001656104598327565004771515257146e-0008,
192 -3.694275157433032553021246812379258781665e-0010,
193 1.434364544206266624252820889648445263842e-0012,
194 -3.852064731859936455895036286874139896861e-0015,
195 7.182052899726138381739945881914874579696e-0018,
196 -9.060556574619677567323741194079797987200e-0021,
197 7.124435467408860515265552217131230511455e-0024,
198 -2.709726774636397615328813121715432044771e-0027,
199 };
200
201 static const GENERIC v0[5] = {
202 1.0,
203 4.678678931512549002587702477349214886475e-0003,
204 9.486828955529948534822800829497565178985e-0006,
205 1.001495929158861646659010844136682454906e-0008,
206 4.725338116256021660204443235685358593611e-0012,
207 };
208
209 GENERIC
210 y0(GENERIC x)
211 {
212 GENERIC z, /* d, */ s, c, ss, cc, u, v;
213 int i;
214
215 if (isnan(x))
216 return (x * x); /* + -> * for Cheetah */
217
218 if (x <= zero) {
219 if (x == zero)
220 /* d= -one/(x-x); */
221 return (_SVID_libm_err(x, x, 8));
222 else
223 /* d = zero/(x-x); */
224 return (_SVID_libm_err(x, x, 9));
225 }
226
227 if (x > 8.0) {
228 if (!finite(x))
229 return (zero);
230
231 s = sin(x);
232 c = cos(x);
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 */
248 if (x > 8.9e307) { /* x+x may overflow */
249 ss = s - c;
250 cc = s + c;
251 } else if (signbit(s) != signbit(c)) {
252 ss = s - c;
253 cc = -cos(x + x) / ss;
254 } else {
255 cc = s + c;
256 ss = -cos(x + x) / cc;
257 }
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 */
263 if (x > 1.0e40)
264 z = (invsqrtpi * ss) / sqrt(x);
265 else
266 z = invsqrtpi * (pzero(x) * ss + qzero(x) * cc) /
267 sqrt(x);
268
269 if (x > X_TLOSS)
270 return (_SVID_libm_err(x, z, 35));
271 else
272 return (z);
273 }
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)));
285 }
286
287 static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
288 .4861344183386052721391238447e5, .1377662549407112278133438945e6,
289 .1222466364088289731869114004e6, .4107070084315176135583353374e5,
290 .5026073801860637125889039915e4, .1783193659125479654541542419e3,
291 .88010344055383421691677564e0,
292 };
293
294 static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
295 .4861344183386052721414037058e5, .1378196632630384670477582699e6,
296 .1223967185341006542748936787e6, .4120150243795353639995862617e5,
297 .5068271181053546392490184353e4, .1829817905472769960535671664e3,
298 1.0,
299 };
300
301 static const GENERIC huge = 1.0e10;
302 static GENERIC
303 pzero(GENERIC x)
304 {
305 GENERIC s, r, t, z;
306 int i;
307
308 if (x > huge)
309 return (one);
310
311 t = eight / x;
312 z = t * t;
313 r = pr[5] + z * pr[6];
314 s = ps[5] + z;
315
316 for (i = 4; i >= 0; i--) {
317 r = r * z + pr[i];
318 s = s * z + ps[i];
319 }
320
321 return (r / s);
322 }
323
324 static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
325 -.1731210995701068539185611951e3, -.5522559165936166961235240613e3,
326 -.5604935606637346590614529613e3, -.2200430300226009379477365011e3,
327 -.323869355375648849771296746e2, -.14294979207907956223499258e1,
328 -.834690374102384988158918e-2,
329 };
330
331 static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
332 .1107975037248683865326709645e5, .3544581680627082674651471873e5,
333 .3619118937918394132179019059e5, .1439895563565398007471485822e5,
334 .2190277023344363955930226234e4, .106695157020407986137501682e3,
335 1.0,
336 };
337
338 static GENERIC
339 qzero(GENERIC x)
340 {
341 GENERIC s, r, t, z;
342 int i;
343
344 if (x > huge)
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
352 for (i = 4; i >= 0; i--) {
353 r = r * z + qr[i];
354 s = s * z + qs[i];
355 }
356
357 return (t * (r / s));
358 }