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