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_synonyms.h"
44 #include "libm_protos.h"
45 #include <math.h>
46 #include <values.h>
47
48 #define GENERIC double
49 static const GENERIC
50 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), qzero(GENERIC);
59 static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
60 -2.500000000000003622131880894830476755537e-0001,
61 1.095597547334830263234433855932375353303e-0002,
62 -1.819734750463320921799187258987098087697e-0004,
63 9.977001946806131657544212501069893930846e-0007,
64 };
65 static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
66 1.0,
67 1.867609810662950169966782360588199673741e-0002,
68 1.590389206181565490878430827706972074208e-0004,
69 6.520867386742583632375520147714499522721e-0007,
70 };
71 static const GENERIC r1[9] = { /* [1.28,8] */
72 9.999999999999999942156495584397047660949e-0001,
73 -2.389887722731319130476839836908143731281e-0001,
74 1.293359476138939027791270393439493640570e-0002,
75 -2.770985642343140122168852400228563364082e-0004,
76 2.905241575772067678086738389169625218912e-0006,
77 -1.636846356264052597969042009265043251279e-0008,
78 5.072306160724884775085431059052611737827e-0011,
79 -8.187060730684066824228914775146536139112e-0014,
80 5.422219326959949863954297860723723423842e-0017,
81 };
82 static const GENERIC s1[9] = { /* [1.28,8] */
83 1.0,
84 1.101122772686807702762104741932076228349e-0002,
85 6.140169310641649223411427764669143978228e-0005,
86 2.292035877515152097976946119293215705250e-0007,
87 6.356910426504644334558832036362219583789e-0010,
88 1.366626326900219555045096999553948891401e-0012,
89 2.280399586866739522891837985560481180088e-0015,
90 2.801559820648939665270492520004836611187e-0018,
91 2.073101088320349159764410261466350732968e-0021,
92 };
93
94 GENERIC
95 j0(GENERIC x) {
96 GENERIC z, s,c,ss,cc,r,u,v,ox;
97 int i;
98
99 if(isnan(x)) return x*x; /* + -> * for Cheetah */
100 ox= x;
101 x = fabs(x);
102 if(x > 8.0){
103 if(!finite(x)) return zero;
104 s = sin(x);
105 c = cos(x);
106 /* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
107 * where x0 = x-pi/4
108 * Better formula:
109 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
110 * = 1/sqrt(2) * (cos(x) + sin(x))
111 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
112 * = 1/sqrt(2) * (sin(x) - cos(x))
113 * To avoid cancellation, use
114 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
115 * to compute the worse one.
116 */
117 if(x>8.9e307) { /* x+x may overflow */
118 ss = s-c;
119 cc = s+c;
120 } else if(signbit(s)!=signbit(c)) {
121 ss = s - c;
122 cc = -cos(x+x)/ss;
123 } else {
124 cc = s + c;
125 ss = -cos(x+x)/cc;
126 }
127 /*
128 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
129 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
130 */
131 if(x>1.0e40) z= (invsqrtpi*cc)/sqrt(x);
132 else {
133 u = pzero(x); v = qzero(x);
134 z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
135 }
136 /* force to pass SVR4 even the result is wrong (sign) */
137 if (x > X_TLOSS)
138 return _SVID_libm_err(ox,z,34);
139 else
140 return z;
141 }
142 if(x<=small) {
143 if(x<=tiny) return one-x;
144 else return one-x*x*0.25;
145 }
146 z = x*x;
147 if(x<=1.28) {
148 r = r0[0]+z*(r0[1]+z*(r0[2]+z*r0[3]));
149 s = s0[0]+z*(s0[1]+z*(s0[2]+z*s0[3]));
150 return one + z*(r/s);
151 } else {
152 for(r=r1[8],s=s1[8],i=7;i>=0;i--) {
153 r = r*z + r1[i];
154 s = s*z + s1[i];
155 }
156 return(r/s);
157 }
158 }
159
160 static const GENERIC u0[13] = {
161 -7.380429510868722526754723020704317641941e-0002,
162 1.772607102684869924301459663049874294814e-0001,
163 -1.524370666542713828604078090970799356306e-0002,
164 4.650819100693891757143771557629924591915e-0004,
165 -7.125768872339528975036316108718239946022e-0006,
166 6.411017001656104598327565004771515257146e-0008,
167 -3.694275157433032553021246812379258781665e-0010,
168 1.434364544206266624252820889648445263842e-0012,
169 -3.852064731859936455895036286874139896861e-0015,
170 7.182052899726138381739945881914874579696e-0018,
171 -9.060556574619677567323741194079797987200e-0021,
172 7.124435467408860515265552217131230511455e-0024,
173 -2.709726774636397615328813121715432044771e-0027,
174 };
175 static const GENERIC v0[5] = {
176 1.0,
177 4.678678931512549002587702477349214886475e-0003,
178 9.486828955529948534822800829497565178985e-0006,
179 1.001495929158861646659010844136682454906e-0008,
180 4.725338116256021660204443235685358593611e-0012,
181 };
182
183 GENERIC
184 y0(GENERIC x) {
185 GENERIC z, /* d, */ s,c,ss,cc,u,v;
186 int i;
187
188 if(isnan(x)) return x*x; /* + -> * for Cheetah */
189 if(x <= zero){
190 if(x==zero)
191 /* d= -one/(x-x); */
192 return _SVID_libm_err(x,x,8);
193 else
194 /* d = zero/(x-x); */
195 return _SVID_libm_err(x,x,9);
196 }
197 if(x > 8.0){
198 if(!finite(x)) return zero;
199 s = sin(x);
200 c = cos(x);
201 /* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
202 * where x0 = x-pi/4
203 * Better formula:
204 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
205 * = 1/sqrt(2) * (cos(x) + sin(x))
206 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
207 * = 1/sqrt(2) * (sin(x) - cos(x))
208 * To avoid cancellation, use
209 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
210 * to compute the worse one.
211 */
212 if(x>8.9e307) { /* x+x may overflow */
213 ss = s-c;
214 cc = s+c;
215 } else if(signbit(s)!=signbit(c)) {
216 ss = s - c;
217 cc = -cos(x+x)/ss;
218 } else {
219 cc = s + c;
220 ss = -cos(x+x)/cc;
221 }
222 /*
223 * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
224 * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
225 */
226 if(x>1.0e40)
227 z = (invsqrtpi*ss)/sqrt(x);
228 else
229 z = invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x);
230 if (x > X_TLOSS)
231 return _SVID_libm_err(x,z,35);
232 else
233 return z;
234
235 }
236 if(x<=tiny) {
237 return(u0[0] + tpi*log(x));
238 }
239 z = x*x;
240 for(u=u0[12],i=11;i>=0;i--) u = u*z + u0[i];
241 v = v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])));
242 return(u/v + tpi*(j0(x)*log(x)));
243 }
244
245 static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
246 .4861344183386052721391238447e5,
247 .1377662549407112278133438945e6,
248 .1222466364088289731869114004e6,
249 .4107070084315176135583353374e5,
250 .5026073801860637125889039915e4,
251 .1783193659125479654541542419e3,
252 .88010344055383421691677564e0,
253 };
254 static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
255 .4861344183386052721414037058e5,
256 .1378196632630384670477582699e6,
257 .1223967185341006542748936787e6,
258 .4120150243795353639995862617e5,
259 .5068271181053546392490184353e4,
260 .1829817905472769960535671664e3,
261 1.0,
262 };
263 static const GENERIC huge = 1.0e10;
264
265 static GENERIC
266 pzero(GENERIC x) {
267 GENERIC s,r,t,z;
268 int i;
269 if(x>huge) return one;
270 t = eight/x; z = t*t;
271 r = pr[5]+z*pr[6];
272 s = ps[5]+z;
273 for(i=4;i>=0;i--) {
274 r = r*z + pr[i];
275 s = s*z + ps[i];
276 }
277 return r/s;
278 }
279
280 static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
281 -.1731210995701068539185611951e3,
282 -.5522559165936166961235240613e3,
283 -.5604935606637346590614529613e3,
284 -.2200430300226009379477365011e3,
285 -.323869355375648849771296746e2,
286 -.14294979207907956223499258e1,
287 -.834690374102384988158918e-2,
288 };
289 static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
290 .1107975037248683865326709645e5,
291 .3544581680627082674651471873e5,
292 .3619118937918394132179019059e5,
293 .1439895563565398007471485822e5,
294 .2190277023344363955930226234e4,
295 .106695157020407986137501682e3,
296 1.0,
297 };
298
299 static GENERIC
300 qzero(GENERIC x) {
301 GENERIC s,r,t,z;
302 int i;
303 if(x>huge) return -0.125/x;
304 t = eight/x; z = t*t;
305 r = qr[5]+z*qr[6];
306 s = qs[5]+z;
307 for(i=4;i>=0;i--) {
308 r = r*z + qr[i];
309 s = s*z + qs[i];
310 }
311 return t*(r/s);
312 }