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 * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27 * Use is subject to license terms.
28 */
29
30/*
31 * Floating point Bessel's function of the first and second kinds
32 * of order zero: j0(x),y0(x);
33 *
34 * Special cases:
35 *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
36 *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
37 */
38
39#pragma weak __j0 = j0
40#pragma weak __y0 = y0
41
42#include "libm.h"
43#include "libm_protos.h"
44#include <math.h>
45#include <values.h>
46
47#define	GENERIC double
48static const GENERIC
49zero    = 0.0,
50small	= 1.0e-5,
51tiny	= 1.0e-18,
52one	= 1.0,
53eight   = 8.0,
54invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
55tpi	= 0.636619772367581343075535053490057448;
56
57static GENERIC pzero(GENERIC), qzero(GENERIC);
58static const GENERIC r0[4] = {	/* [1.e-5, 1.28] */
59	-2.500000000000003622131880894830476755537e-0001,
60	1.095597547334830263234433855932375353303e-0002,
61	-1.819734750463320921799187258987098087697e-0004,
62	9.977001946806131657544212501069893930846e-0007,
63};
64static const GENERIC s0[4] = {	/* [1.e-5, 1.28] */
65	1.0,
66	1.867609810662950169966782360588199673741e-0002,
67	1.590389206181565490878430827706972074208e-0004,
68	6.520867386742583632375520147714499522721e-0007,
69};
70static const GENERIC r1[9] = {	/* [1.28,8] */
71	9.999999999999999942156495584397047660949e-0001,
72	-2.389887722731319130476839836908143731281e-0001,
73	1.293359476138939027791270393439493640570e-0002,
74	-2.770985642343140122168852400228563364082e-0004,
75	2.905241575772067678086738389169625218912e-0006,
76	-1.636846356264052597969042009265043251279e-0008,
77	5.072306160724884775085431059052611737827e-0011,
78	-8.187060730684066824228914775146536139112e-0014,
79	5.422219326959949863954297860723723423842e-0017,
80};
81static const GENERIC s1[9] = {	/* [1.28,8] */
82	1.0,
83	1.101122772686807702762104741932076228349e-0002,
84	6.140169310641649223411427764669143978228e-0005,
85	2.292035877515152097976946119293215705250e-0007,
86	6.356910426504644334558832036362219583789e-0010,
87	1.366626326900219555045096999553948891401e-0012,
88	2.280399586866739522891837985560481180088e-0015,
89	2.801559820648939665270492520004836611187e-0018,
90	2.073101088320349159764410261466350732968e-0021,
91};
92
93GENERIC
94j0(GENERIC x) {
95	GENERIC z, s, c, ss, cc, r, u, v, ox;
96	int i;
97
98	if (isnan(x))
99		return (x*x);	/* + -> * for Cheetah */
100	ox = x;
101	x = fabs(x);
102	if (x > 8.0) {
103		if (!finite(x))
104			return (zero);
105		s = sin(x);
106		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	 */
119		if (x > 8.9e307) {	/* x+x may overflow */
120			ss = s-c;
121			cc = s+c;
122		} else if (signbit(s) != signbit(c)) {
123			ss = s - c;
124			cc = -cos(x+x)/ss;
125		} else {
126			cc = s + c;
127			ss = -cos(x+x)/cc;
128		}
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);
137		}
138	/* force to pass SVR4 even the result is wrong (sign) */
139		if (x > X_TLOSS)
140		    return (_SVID_libm_err(ox, z, 34));
141		else
142		    return (z);
143	}
144	if (x <= small) {
145	    if (x <= tiny)
146			return (one-x);
147	    else
148			return (one-x*x*0.25);
149	}
150	z = x*x;
151	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));
155	} 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);
161	}
162}
163
164static const GENERIC u0[13] = {
165	-7.380429510868722526754723020704317641941e-0002,
166	1.772607102684869924301459663049874294814e-0001,
167	-1.524370666542713828604078090970799356306e-0002,
168	4.650819100693891757143771557629924591915e-0004,
169	-7.125768872339528975036316108718239946022e-0006,
170	6.411017001656104598327565004771515257146e-0008,
171	-3.694275157433032553021246812379258781665e-0010,
172	1.434364544206266624252820889648445263842e-0012,
173	-3.852064731859936455895036286874139896861e-0015,
174	7.182052899726138381739945881914874579696e-0018,
175	-9.060556574619677567323741194079797987200e-0021,
176	7.124435467408860515265552217131230511455e-0024,
177	-2.709726774636397615328813121715432044771e-0027,
178};
179static const GENERIC v0[5] = {
180	1.0,
181	4.678678931512549002587702477349214886475e-0003,
182	9.486828955529948534822800829497565178985e-0006,
183	1.001495929158861646659010844136682454906e-0008,
184	4.725338116256021660204443235685358593611e-0012,
185};
186
187GENERIC
188y0(GENERIC x) {
189	GENERIC z, /* d, */ s, c, ss, cc, u, v;
190	int i;
191
192	if (isnan(x))
193		return (x*x);	/* + -> * for Cheetah */
194	if (x <= zero) {
195		if (x == zero)
196		    /* d= -one/(x-x); */
197		    return (_SVID_libm_err(x, x, 8));
198		else
199		    /* d = zero/(x-x); */
200		    return (_SVID_libm_err(x, x, 9));
201	}
202	if (x > 8.0) {
203		if (!finite(x))
204			return (zero);
205		s = sin(x);
206		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	 */
219		if (x > 8.9e307) {	/* x+x may overflow */
220			ss = s-c;
221			cc = s+c;
222		} else if (signbit(s) != signbit(c)) {
223			ss = s - c;
224			cc = -cos(x+x)/ss;
225		} else {
226			cc = s + c;
227			ss = -cos(x+x)/cc;
228		}
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	 */
233		if (x > 1.0e40)
234		    z = (invsqrtpi*ss)/sqrt(x);
235		else
236		    z =  invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x);
237		if (x > X_TLOSS)
238		    return (_SVID_libm_err(x, z, 35));
239		else
240		    return (z);
241
242	}
243	if (x <= tiny) {
244	    return (u0[0] + tpi*log(x));
245	}
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)));
250}
251
252static const GENERIC pr[7] = {	/* [8 -- inf]  pzero 6550 */
253	.4861344183386052721391238447e5,
254	.1377662549407112278133438945e6,
255	.1222466364088289731869114004e6,
256	.4107070084315176135583353374e5,
257	.5026073801860637125889039915e4,
258	.1783193659125479654541542419e3,
259	.88010344055383421691677564e0,
260};
261static const GENERIC ps[7] = {	/* [8 -- inf] pzero 6550 */
262	.4861344183386052721414037058e5,
263	.1378196632630384670477582699e6,
264	.1223967185341006542748936787e6,
265	.4120150243795353639995862617e5,
266	.5068271181053546392490184353e4,
267	.1829817905472769960535671664e3,
268	1.0,
269};
270static const GENERIC huge    = 1.0e10;
271
272static GENERIC
273pzero(GENERIC x) {
274	GENERIC s, r, t, z;
275	int i;
276	if (x > huge)
277		return (one);
278	t = eight/x; z = t*t;
279	r = pr[5]+z*pr[6];
280	s = ps[5]+z;
281	for (i = 4; i >= 0; i--) {
282	    r = r*z + pr[i];
283	    s = s*z + ps[i];
284	}
285	return (r/s);
286}
287
288static const GENERIC qr[7] = {	/* [8 -- inf]  qzero 6950 */
289	-.1731210995701068539185611951e3,
290	-.5522559165936166961235240613e3,
291	-.5604935606637346590614529613e3,
292	-.2200430300226009379477365011e3,
293	-.323869355375648849771296746e2,
294	-.14294979207907956223499258e1,
295	-.834690374102384988158918e-2,
296};
297static const GENERIC qs[7] = {	/* [8 -- inf] qzero 6950 */
298	.1107975037248683865326709645e5,
299	.3544581680627082674651471873e5,
300	.3619118937918394132179019059e5,
301	.1439895563565398007471485822e5,
302	.2190277023344363955930226234e4,
303	.106695157020407986137501682e3,
304	1.0,
305};
306
307static GENERIC
308qzero(GENERIC x) {
309	GENERIC s, r, t, z;
310	int i;
311	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;
316	for (i = 4; i >= 0; i--) {
317	    r = r*z + qr[i];
318	    s = s*z + qs[i];
319	}
320	return (t*(r/s));
321}
322