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: j1(x),y1(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 __j1 = j1
41 #pragma weak __y1 = y1
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-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);
59
60 static const GENERIC r0[4] = {
61 -6.250000000000002203053200981413218949548e-0002,
62 1.600998455640072901321605101981501263762e-0003,
63 -1.963888815948313758552511884390162864930e-0005,
64 8.263917341093549759781339713418201620998e-0008,
65 };
66
67 static const GENERIC s0[7] = {
68 1.0e0,
69 1.605069137643004242395356851797873766927e-0002,
70 1.149454623251299996428500249509098499383e-0004,
71 3.849701673735260970379681807910852327825e-0007,
72 };
73
74 static const GENERIC r1[12] = {
75 4.999999999999999995517408894340485471724e-0001,
76 -6.003825028120475684835384519945468075423e-0002,
77 2.301719899263321828388344461995355419832e-0003,
78 -4.208494869238892934859525221654040304068e-0005,
79 4.377745135188837783031540029700282443388e-0007,
80 -2.854106755678624335145364226735677754179e-0009,
81 1.234002865443952024332943901323798413689e-0011,
82 -3.645498437039791058951273508838177134310e-0014,
83 7.404320596071797459925377103787837414422e-0017,
84 -1.009457448277522275262808398517024439084e-0019,
85 8.520158355824819796968771418801019930585e-0023,
86 -3.458159926081163274483854614601091361424e-0026,
87 };
88
89 static const GENERIC s1[5] = {
90 1.0e0,
91 4.923499437590484879081138588998986303306e-0003,
92 1.054389489212184156499666953501976688452e-0005,
93 1.180768373106166527048240364872043816050e-0008,
94 5.942665743476099355323245707680648588540e-0012,
95 };
96
97 GENERIC
98 j1(GENERIC x)
99 {
100 GENERIC z, d, s, c, ss, cc, r;
101 int i, sgn;
102
103 if (!finite(x))
104 return (one / x);
105
106 sgn = signbit(x);
107 x = fabs(x);
108
109 if (x > 8.00) {
110 s = sin(x);
111 c = cos(x);
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 */
127 if (x > 8.9e307) { /* x+x may overflow */
128 ss = -s - c;
129 cc = s - c;
130 } else if (signbit(s) != signbit(c)) {
131 cc = s - c;
132 ss = cos(x + x) / cc;
133 } else {
134 ss = -s - c;
135 cc = cos(x + x) / ss;
136 }
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 */
142 if (x > 1.0e40)
143 d = (invsqrtpi * cc) / sqrt(x);
144 else
145 d = invsqrtpi * (pone(x) * cc - qone(x) * ss) / sqrt(x);
146
147 if (x > X_TLOSS) {
148 if (sgn != 0) {
149 d = -d;
150 x = -x;
151 }
152
153 return (_SVID_libm_err(x, d, 36));
154 } else if (sgn == 0) {
155 return (d);
156 } else {
157 return (-d);
158 }
159 }
160
161 if (x <= small) {
162 if (x <= tiny)
163 d = 0.5 * x;
164 else
165 d = x * (0.5 - x * x * 0.125);
166
167 if (sgn == 0)
168 return (d);
169 else
170 return (-d);
171 }
172
173 z = x * x;
174
175 if (x < 1.28) {
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));
185 } else {
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);
193 }
194
195 if (sgn == 0)
196 return (d);
197 else
198 return (-d);
199 }
200
201 static const GENERIC u0[4] = {
202 -1.960570906462389461018983259589655961560e-0001,
203 4.931824118350661953459180060007970291139e-0002,
204 -1.626975871565393656845930125424683008677e-0003,
205 1.359657517926394132692884168082224258360e-0005,
206 };
207
208 static const GENERIC v0[5] = {
209 1.0e0,
210 2.565807214838390835108224713630901653793e-0002,
211 3.374175208978404268650522752520906231508e-0004,
212 2.840368571306070719539936935220728843177e-0006,
213 1.396387402048998277638900944415752207592e-0008,
214 };
215
216 static const GENERIC u1[12] = {
217 -1.960570906462389473336339614647555351626e-0001,
218 5.336268030335074494231369159933012844735e-0002,
219 -2.684137504382748094149184541866332033280e-0003,
220 5.737671618979185736981543498580051903060e-0005,
221 -6.642696350686335339171171785557663224892e-0007,
222 4.692417922568160354012347591960362101664e-0009,
223 -2.161728635907789319335231338621412258355e-0011,
224 6.727353419738316107197644431844194668702e-0014,
225 -1.427502986803861372125234355906790573422e-0016,
226 2.020392498726806769468143219616642940371e-0019,
227 -1.761371948595104156753045457888272716340e-0022,
228 7.352828391941157905175042420249225115816e-0026,
229 };
230
231 static const GENERIC v1[5] = {
232 1.0e0,
233 5.029187436727947764916247076102283399442e-0003,
234 1.102693095808242775074856548927801750627e-0005,
235 1.268035774543174837829534603830227216291e-0008,
236 6.579416271766610825192542295821308730206e-0012,
237 };
238
239 GENERIC
240 y1(GENERIC x)
241 {
242 GENERIC z, d, s, c, ss, cc, u, v;
243 int i;
244
245 if (isnan(x))
246 return (x * x); /* + -> * for Cheetah */
247
248 if (x <= zero) {
249 if (x == zero)
250 /* return -one/zero; */
251 return (_SVID_libm_err(x, x, 10));
252 else
253 /* return zero/zero; */
254 return (_SVID_libm_err(x, x, 11));
255 }
256
257 if (x > 8.0) {
258 if (!finite(x))
259 return (zero);
260
261 s = sin(x);
262 c = cos(x);
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 */
278 if (x > 8.9e307) { /* x+x may overflow */
279 ss = -s - c;
280 cc = s - c;
281 } else if (signbit(s) != signbit(c)) {
282 cc = s - c;
283 ss = cos(x + x) / cc;
284 } else {
285 ss = -s - c;
286 cc = cos(x + x) / ss;
287 }
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 */
293 if (x > 1.0e91)
294 d = (invsqrtpi * ss) / sqrt(x);
295 else
296 d = invsqrtpi * (pone(x) * ss + qone(x) * cc) / sqrt(x);
297
298 if (x > X_TLOSS)
299 return (_SVID_libm_err(x, d, 37));
300 else
301 return (d);
302 }
303
304 if (x <= tiny)
305 return (-tpi / x);
306
307 z = x * x;
308
309 if (x < 1.28) {
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 }
317 } else {
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])));
322 }
323
324 return (x * (u / v) + tpi * (j1(x) * log(x) - one / x));
325 }
326
327 static const GENERIC pr0[6] = {
328 -.4435757816794127857114720794e7, -.9942246505077641195658377899e7,
329 -.6603373248364939109255245434e7, -.1523529351181137383255105722e7,
330 -.1098240554345934672737413139e6, -.1611616644324610116477412898e4,
331 };
332
333 static const GENERIC ps0[6] = {
334 -.4435757816794127856828016962e7, -.9934124389934585658967556309e7,
335 -.6585339479723087072826915069e7, -.1511809506634160881644546358e7,
336 -.1072638599110382011903063867e6, -.1455009440190496182453565068e4,
337 };
338
339 static const GENERIC huge = 1.0e10;
340 static GENERIC
341 pone(GENERIC x)
342 {
343 GENERIC s, r, t, z;
344 int i;
345
346 /* assume x > 8 */
347 if (x > huge)
348 return (one);
349
350 t = 8.0 / x;
351 z = t * t;
352 r = pr0[5];
353 s = ps0[5] + z;
354
355 for (i = 4; i >= 0; i--) {
356 r = z * r + pr0[i];
357 s = z * s + ps0[i];
358 }
359
360 return (r / s);
361 }
362
363 static const GENERIC qr0[6] = {
364 0.3322091340985722351859704442e5, 0.8514516067533570196555001171e5,
365 0.6617883658127083517939992166e5, 0.1849426287322386679652009819e5,
366 0.1706375429020768002061283546e4, 0.3526513384663603218592175580e2,
367 };
368
369 static const GENERIC qs0[6] = {
370 0.7087128194102874357377502472e6, 0.1819458042243997298924553839e7,
371 0.1419460669603720892855755253e7, 0.4002944358226697511708610813e6,
372 0.3789022974577220264142952256e5, 0.8638367769604990967475517183e3,
373 };
374
375 static GENERIC
376 qone(GENERIC x)
377 {
378 GENERIC s, r, t, z;
379 int i;
380
381 if (x > huge)
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;
389
390 for (i = 4; i >= 0; i--) {
391 r = z * r + qr0[i];
392 s = z * s + qs0[i];
393 }
394
395 return (t * (r / s));
396 }