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