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 * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
23 */
24/*
25 * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
26 * Use is subject to license terms.
27 */
28
29#pragma weak __j0f = j0f
30#pragma weak __j1f = j1f
31#pragma weak __jnf = jnf
32#pragma weak __y0f = y0f
33#pragma weak __y1f = y1f
34#pragma weak __ynf = ynf
35
36#include "libm.h"
37#include <float.h>
38
39#if defined(__i386) && !defined(__amd64)
40extern int __swapRP(int);
41#endif
42
43static const float
44	zerof	= 0.0f,
45	onef	= 1.0f;
46
47static const double C[] = {
48	0.0,
49	-0.125,
50	0.25,
51	0.375,
52	0.5,
53	1.0,
54	2.0,
55	8.0,
56	0.5641895835477562869480794515607725858441,	/* 1/sqrt(pi) */
57	0.636619772367581343075535053490057448,	/* 2/pi */
58	1.0e9,
59};
60
61#define	zero	C[0]
62#define	neighth	C[1]
63#define	quarter	C[2]
64#define	three8	C[3]
65#define	half	C[4]
66#define	one	C[5]
67#define	two	C[6]
68#define	eight   C[7]
69#define	isqrtpi	C[8]
70#define	tpi	C[9]
71#define	big	C[10]
72
73static const double Cj0y0[] = {
74	0.4861344183386052721391238447e5,	/* pr */
75	0.1377662549407112278133438945e6,
76	0.1222466364088289731869114004e6,
77	0.4107070084315176135583353374e5,
78	0.5026073801860637125889039915e4,
79	0.1783193659125479654541542419e3,
80	0.88010344055383421691677564e0,
81	0.4861344183386052721414037058e5,	/* ps */
82	0.1378196632630384670477582699e6,
83	0.1223967185341006542748936787e6,
84	0.4120150243795353639995862617e5,
85	0.5068271181053546392490184353e4,
86	0.1829817905472769960535671664e3,
87	1.0,
88	-0.1731210995701068539185611951e3,	/* qr */
89	-0.5522559165936166961235240613e3,
90	-0.5604935606637346590614529613e3,
91	-0.2200430300226009379477365011e3,
92	-0.323869355375648849771296746e2,
93	-0.14294979207907956223499258e1,
94	-0.834690374102384988158918e-2,
95	0.1107975037248683865326709645e5,	/* qs */
96	0.3544581680627082674651471873e5,
97	0.3619118937918394132179019059e5,
98	0.1439895563565398007471485822e5,
99	0.2190277023344363955930226234e4,
100	0.106695157020407986137501682e3,
101	1.0,
102};
103
104#define	pr	Cj0y0
105#define	ps	(Cj0y0+7)
106#define	qr	(Cj0y0+14)
107#define	qs	(Cj0y0+21)
108
109static const double Cj0[] = {
110	-2.500000000000003622131880894830476755537e-0001,	/* r0 */
111	1.095597547334830263234433855932375353303e-0002,
112	-1.819734750463320921799187258987098087697e-0004,
113	9.977001946806131657544212501069893930846e-0007,
114	1.0,							/* s0 */
115	1.867609810662950169966782360588199673741e-0002,
116	1.590389206181565490878430827706972074208e-0004,
117	6.520867386742583632375520147714499522721e-0007,
118	9.999999999999999942156495584397047660949e-0001,	/* r1 */
119	-2.389887722731319130476839836908143731281e-0001,
120	1.293359476138939027791270393439493640570e-0002,
121	-2.770985642343140122168852400228563364082e-0004,
122	2.905241575772067678086738389169625218912e-0006,
123	-1.636846356264052597969042009265043251279e-0008,
124	5.072306160724884775085431059052611737827e-0011,
125	-8.187060730684066824228914775146536139112e-0014,
126	5.422219326959949863954297860723723423842e-0017,
127	1.0,							/* s1 */
128	1.101122772686807702762104741932076228349e-0002,
129	6.140169310641649223411427764669143978228e-0005,
130	2.292035877515152097976946119293215705250e-0007,
131	6.356910426504644334558832036362219583789e-0010,
132	1.366626326900219555045096999553948891401e-0012,
133	2.280399586866739522891837985560481180088e-0015,
134	2.801559820648939665270492520004836611187e-0018,
135	2.073101088320349159764410261466350732968e-0021,
136};
137
138#define	r0	Cj0
139#define	s0	(Cj0+4)
140#define	r1	(Cj0+8)
141#define	s1	(Cj0+17)
142
143static const double Cy0[] = {
144	-7.380429510868722526754723020704317641941e-0002,	/* u0 */
145	1.772607102684869924301459663049874294814e-0001,
146	-1.524370666542713828604078090970799356306e-0002,
147	4.650819100693891757143771557629924591915e-0004,
148	-7.125768872339528975036316108718239946022e-0006,
149	6.411017001656104598327565004771515257146e-0008,
150	-3.694275157433032553021246812379258781665e-0010,
151	1.434364544206266624252820889648445263842e-0012,
152	-3.852064731859936455895036286874139896861e-0015,
153	7.182052899726138381739945881914874579696e-0018,
154	-9.060556574619677567323741194079797987200e-0021,
155	7.124435467408860515265552217131230511455e-0024,
156	-2.709726774636397615328813121715432044771e-0027,
157	1.0,							/* v0 */
158	4.678678931512549002587702477349214886475e-0003,
159	9.486828955529948534822800829497565178985e-0006,
160	1.001495929158861646659010844136682454906e-0008,
161	4.725338116256021660204443235685358593611e-0012,
162};
163
164#define	u0	Cy0
165#define	v0	(Cy0+13)
166
167static const double Cj1y1[] = {
168	-0.4435757816794127857114720794e7,	/* pr0 */
169	-0.9942246505077641195658377899e7,
170	-0.6603373248364939109255245434e7,
171	-0.1523529351181137383255105722e7,
172	-0.1098240554345934672737413139e6,
173	-0.1611616644324610116477412898e4,
174	-0.4435757816794127856828016962e7,	/* ps0 */
175	-0.9934124389934585658967556309e7,
176	-0.6585339479723087072826915069e7,
177	-0.1511809506634160881644546358e7,
178	-0.1072638599110382011903063867e6,
179	-0.1455009440190496182453565068e4,
180	0.3322091340985722351859704442e5,	/* qr0 */
181	0.8514516067533570196555001171e5,
182	0.6617883658127083517939992166e5,
183	0.1849426287322386679652009819e5,
184	0.1706375429020768002061283546e4,
185	0.3526513384663603218592175580e2,
186	0.7087128194102874357377502472e6,	/* qs0 */
187	0.1819458042243997298924553839e7,
188	0.1419460669603720892855755253e7,
189	0.4002944358226697511708610813e6,
190	0.3789022974577220264142952256e5,
191	0.8638367769604990967475517183e3,
192};
193
194#define	pr0	Cj1y1
195#define	ps0	(Cj1y1+6)
196#define	qr0	(Cj1y1+12)
197#define	qs0	(Cj1y1+18)
198
199static const double Cj1[] = {
200	-6.250000000000002203053200981413218949548e-0002,	/* a0 */
201	1.600998455640072901321605101981501263762e-0003,
202	-1.963888815948313758552511884390162864930e-0005,
203	8.263917341093549759781339713418201620998e-0008,
204	1.0e0,							/* b0 */
205	1.605069137643004242395356851797873766927e-0002,
206	1.149454623251299996428500249509098499383e-0004,
207	3.849701673735260970379681807910852327825e-0007,
208	4.999999999999999995517408894340485471724e-0001,
209	-6.003825028120475684835384519945468075423e-0002,
210	2.301719899263321828388344461995355419832e-0003,
211	-4.208494869238892934859525221654040304068e-0005,
212	4.377745135188837783031540029700282443388e-0007,
213	-2.854106755678624335145364226735677754179e-0009,
214	1.234002865443952024332943901323798413689e-0011,
215	-3.645498437039791058951273508838177134310e-0014,
216	7.404320596071797459925377103787837414422e-0017,
217	-1.009457448277522275262808398517024439084e-0019,
218	8.520158355824819796968771418801019930585e-0023,
219	-3.458159926081163274483854614601091361424e-0026,
220	1.0e0,							/* b1 */
221	4.923499437590484879081138588998986303306e-0003,
222	1.054389489212184156499666953501976688452e-0005,
223	1.180768373106166527048240364872043816050e-0008,
224	5.942665743476099355323245707680648588540e-0012,
225};
226
227#define	a0	Cj1
228#define	b0	(Cj1+4)
229#define	a1	(Cj1+8)
230#define	b1	(Cj1+20)
231
232static const double Cy1[] = {
233	-1.960570906462389461018983259589655961560e-0001,	/* c0 */
234	4.931824118350661953459180060007970291139e-0002,
235	-1.626975871565393656845930125424683008677e-0003,
236	1.359657517926394132692884168082224258360e-0005,
237	1.0e0,							/* d0 */
238	2.565807214838390835108224713630901653793e-0002,
239	3.374175208978404268650522752520906231508e-0004,
240	2.840368571306070719539936935220728843177e-0006,
241	1.396387402048998277638900944415752207592e-0008,
242	-1.960570906462389473336339614647555351626e-0001,	/* c1 */
243	5.336268030335074494231369159933012844735e-0002,
244	-2.684137504382748094149184541866332033280e-0003,
245	5.737671618979185736981543498580051903060e-0005,
246	-6.642696350686335339171171785557663224892e-0007,
247	4.692417922568160354012347591960362101664e-0009,
248	-2.161728635907789319335231338621412258355e-0011,
249	6.727353419738316107197644431844194668702e-0014,
250	-1.427502986803861372125234355906790573422e-0016,
251	2.020392498726806769468143219616642940371e-0019,
252	-1.761371948595104156753045457888272716340e-0022,
253	7.352828391941157905175042420249225115816e-0026,
254	1.0e0,							/* d1 */
255	5.029187436727947764916247076102283399442e-0003,
256	1.102693095808242775074856548927801750627e-0005,
257	1.268035774543174837829534603830227216291e-0008,
258	6.579416271766610825192542295821308730206e-0012,
259};
260
261#define	c0	Cy1
262#define	d0	(Cy1+4)
263#define	c1	(Cy1+9)
264#define	d1	(Cy1+21)
265
266
267/* core of j0f computation; assumes fx is finite */
268static double
269__k_j0f(float fx)
270{
271	double	x, z, s, c, ss, cc, r, t, p0, q0;
272	int	ix, i;
273
274	ix = *(int *)&fx & ~0x80000000;
275	x = fabs((double)fx);
276	if (ix > 0x41000000) {
277		/* x > 8; see comments in j0.c */
278		s = sin(x);
279		c = cos(x);
280		if (signbit(s) != signbit(c)) {
281			ss = s - c;
282			cc = -cos(x + x) / ss;
283		} else {
284			cc = s + c;
285			ss = -cos(x + x) / cc;
286		}
287		if (ix > 0x501502f9) {
288			/* x > 1.0e10 */
289			p0 = one;
290			q0 = neighth / x;
291		} else {
292			t = eight / x;
293			z = t * t;
294			p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
295			    z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
296			    (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
297			    z * (ps[4] + z * (ps[5] + z))))));
298			q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
299			    z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
300			    (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
301			    z * (qs[4] + z * (qs[5] + z))))))) * t;
302		}
303		return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
304	}
305	if (ix <= 0x3727c5ac) {
306		/* x <= 1.0e-5 */
307		if (ix <= 0x219392ef) /* x <= 1.0e-18 */
308			return (one - x);
309		return (one - x * x * quarter);
310	}
311	z = x * x;
312	if (ix <= 0x3fa3d70a) {
313		/* x <= 1.28 */
314		r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
315		s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
316		return (one + z * (r / s));
317	}
318	r = r1[8];
319	s = s1[8];
320	for (i = 7; i >= 0; i--) {
321		r = r * z + r1[i];
322		s = s * z + s1[i];
323	}
324	return (r / s);
325}
326
327float
328j0f(float fx)
329{
330	float	f;
331	int	ix;
332#if defined(__i386) && !defined(__amd64)
333	int	rp;
334#endif
335
336	ix = *(int *)&fx & ~0x80000000;
337	if (ix >= 0x7f800000) {			/* nan or inf */
338		if (ix > 0x7f800000)
339			return (fx * fx);
340		return (zerof);
341	}
342
343#if defined(__i386) && !defined(__amd64)
344	rp = __swapRP(fp_extended);
345#endif
346	f = (float)__k_j0f(fx);
347#if defined(__i386) && !defined(__amd64)
348	if (rp != fp_extended)
349		(void) __swapRP(rp);
350#endif
351	return (f);
352}
353
354/* core of y0f computation; assumes fx is finite and positive */
355static double
356__k_y0f(float fx)
357{
358	double	x, z, s, c, ss, cc, t, p0, q0, u, v;
359	int	ix, i;
360
361	ix = *(int *)&fx;
362	x = (double)fx;
363	if (ix > 0x41000000) {
364		/* x > 8; see comments in j0.c */
365		s = sin(x);
366		c = cos(x);
367		if (signbit(s) != signbit(c)) {
368			ss = s - c;
369			cc = -cos(x + x) / ss;
370		} else {
371			cc = s + c;
372			ss = -cos(x + x) / cc;
373		}
374		if (ix > 0x501502f9) {
375			/* x > 1.0e10 */
376			p0 = one;
377			q0 = neighth / x;
378		} else {
379			t = eight / x;
380			z = t * t;
381			p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
382			    z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
383			    (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
384			    z * (ps[4] + z * (ps[5] + z))))));
385			q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
386			    z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
387			    (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
388			    z * (qs[4] + z * (qs[5] + z))))))) * t;
389		}
390		return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
391	}
392	if (ix <= 0x219392ef) /* x <= 1.0e-18 */
393		return (u0[0] + tpi * log(x));
394	z = x * x;
395	u = u0[12];
396	for (i = 11; i >= 0; i--)
397		u = u * z + u0[i];
398	v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
399	return (u / v + tpi * (__k_j0f(fx) * log(x)));
400}
401
402float
403y0f(float fx)
404{
405	float	f;
406	int	ix;
407#if defined(__i386) && !defined(__amd64)
408	int	rp;
409#endif
410
411	ix = *(int *)&fx;
412	if ((ix & ~0x80000000) > 0x7f800000)	/* nan */
413		return (fx * fx);
414	if (ix <= 0) {				/* zero or negative */
415		if ((ix << 1) == 0)
416			return (-onef / zerof);
417		return (zerof / zerof);
418	}
419	if (ix == 0x7f800000)			/* +inf */
420		return (zerof);
421
422#if defined(__i386) && !defined(__amd64)
423	rp = __swapRP(fp_extended);
424#endif
425	f = (float)__k_y0f(fx);
426#if defined(__i386) && !defined(__amd64)
427	if (rp != fp_extended)
428		(void) __swapRP(rp);
429#endif
430	return (f);
431}
432
433/* core of j1f computation; assumes fx is finite */
434static double
435__k_j1f(float fx)
436{
437	double	x, z, s, c, ss, cc, r, t, p1, q1;
438	int	i, ix, sgn;
439
440	ix = *(int *)&fx;
441	sgn = (unsigned)ix >> 31;
442	ix &= ~0x80000000;
443	x = fabs((double)fx);
444	if (ix > 0x41000000) {
445		/* x > 8; see comments in j1.c */
446		s = sin(x);
447		c = cos(x);
448		if (signbit(s) != signbit(c)) {
449			cc = s - c;
450			ss = cos(x + x) / cc;
451		} else {
452			ss = -s - c;
453			cc = cos(x + x) / ss;
454		}
455		if (ix > 0x501502f9) {
456			/* x > 1.0e10 */
457			p1 = one;
458			q1 = three8 / x;
459		} else {
460			t = eight / x;
461			z = t * t;
462			p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
463			    (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
464			    (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
465			    (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
466			q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
467			    (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
468			    (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
469			    (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
470		}
471		t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
472		return ((sgn)? -t : t);
473	}
474	if (ix <= 0x3727c5ac) {
475		/* x <= 1.0e-5 */
476		if (ix <= 0x219392ef) /* x <= 1.0e-18 */
477			t = half * x;
478		else
479			t = x * (half + neighth * x * x);
480		return ((sgn)? -t : t);
481	}
482	z = x * x;
483	if (ix < 0x3fa3d70a) {
484		/* x < 1.28 */
485		r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
486		s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
487		t = x * half + x * (z * (r / s));
488	} else {
489		r = a1[11];
490		for (i = 10; i >= 0; i--)
491			r = r * z + a1[i];
492		s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
493		t = x * (r / s);
494	}
495	return ((sgn)? -t : t);
496}
497
498float
499j1f(float fx)
500{
501	float	f;
502	int	ix;
503#if defined(__i386) && !defined(__amd64)
504	int	rp;
505#endif
506
507	ix = *(int *)&fx & ~0x80000000;
508	if (ix >= 0x7f800000)			/* nan or inf */
509		return (onef / fx);
510
511#if defined(__i386) && !defined(__amd64)
512	rp = __swapRP(fp_extended);
513#endif
514	f = (float)__k_j1f(fx);
515#if defined(__i386) && !defined(__amd64)
516	if (rp != fp_extended)
517		(void) __swapRP(rp);
518#endif
519	return (f);
520}
521
522/* core of y1f computation; assumes fx is finite and positive */
523static double
524__k_y1f(float fx)
525{
526	double	x, z, s, c, ss, cc, u, v, p1, q1, t;
527	int	i, ix;
528
529	ix = *(int *)&fx;
530	x = (double)fx;
531	if (ix > 0x41000000) {
532		/* x > 8; see comments in j1.c */
533		s = sin(x);
534		c = cos(x);
535		if (signbit(s) != signbit(c)) {
536			cc = s - c;
537			ss = cos(x + x) / cc;
538		} else {
539			ss = -s - c;
540			cc = cos(x + x) / ss;
541		}
542		if (ix > 0x501502f9) {
543			/* x > 1.0e10 */
544			p1 = one;
545			q1 = three8 / x;
546		} else {
547			t = eight / x;
548			z = t * t;
549			p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
550			    (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
551			    (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
552			    (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
553			q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
554			    (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
555			    (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
556			    (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
557		}
558		return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
559	}
560	if (ix <= 0x219392ef) /* x <= 1.0e-18 */
561		return (-tpi / x);
562	z = x * x;
563	if (ix < 0x3fa3d70a) {
564		/* x < 1.28 */
565		u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
566		v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
567	} else {
568		u = c1[11];
569		for (i = 10; i >= 0; i--)
570			u = u * z + c1[i];
571		v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
572	}
573	return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
574}
575
576float
577y1f(float fx)
578{
579	float	f;
580	int	ix;
581#if defined(__i386) && !defined(__amd64)
582	int	rp;
583#endif
584
585	ix = *(int *)&fx;
586	if ((ix & ~0x80000000) > 0x7f800000)	/* nan */
587		return (fx * fx);
588	if (ix <= 0) {				/* zero or negative */
589		if ((ix << 1) == 0)
590			return (-onef / zerof);
591		return (zerof / zerof);
592	}
593	if (ix == 0x7f800000)			/* +inf */
594		return (zerof);
595
596#if defined(__i386) && !defined(__amd64)
597	rp = __swapRP(fp_extended);
598#endif
599	f = (float)__k_y1f(fx);
600#if defined(__i386) && !defined(__amd64)
601	if (rp != fp_extended)
602		(void) __swapRP(rp);
603#endif
604	return (f);
605}
606
607float
608jnf(int n, float fx)
609{
610	double	a, b, temp, x, z, w, t, q0, q1, h;
611	float	f;
612	int	i, ix, sgn, m, k;
613#if defined(__i386) && !defined(__amd64)
614	int	rp;
615#endif
616
617	if (n < 0) {
618		n = -n;
619		fx = -fx;
620	}
621	if (n == 0)
622		return (j0f(fx));
623	if (n == 1)
624		return (j1f(fx));
625
626	ix = *(int *)&fx;
627	sgn = (n & 1)? ((unsigned)ix >> 31) : 0;
628	ix &= ~0x80000000;
629	if (ix >= 0x7f800000) {		/* nan or inf */
630		if (ix > 0x7f800000)
631			return (fx * fx);
632		return ((sgn)? -zerof : zerof);
633	}
634	if ((ix << 1) == 0)
635		return ((sgn)? -zerof : zerof);
636
637#if defined(__i386) && !defined(__amd64)
638	rp = __swapRP(fp_extended);
639#endif
640	fx = fabsf(fx);
641	x = (double)fx;
642	if ((double)n <= x) {
643		/* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
644		a = __k_j0f(fx);
645		b = __k_j1f(fx);
646		for (i = 1; i < n; i++) {
647			temp = b;
648			b = b * ((double)(i + i) / x) - a;
649			a = temp;
650		}
651		f = (float)b;
652#if defined(__i386) && !defined(__amd64)
653		if (rp != fp_extended)
654			(void) __swapRP(rp);
655#endif
656		return ((sgn)? -f : f);
657	}
658	if (ix < 0x3089705f) {
659		/* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
660		if (n > 6)
661			n = 6;	/* result underflows to zero for n >= 6 */
662		b = t = half * x;
663		a = one;
664		for (i = 2; i <= n; i++) {
665			b *= t;
666			a *= (double)i;
667		}
668		b /= a;
669	} else {
670		/*
671		 * Use the backward recurrence:
672		 *
673		 * 			x      x^2	x^2
674		 *  J(n,x)/J(n-1,x) =  ---- - ------ - ------   .....
675		 *			2n    2(n+1)   2(n+2)
676		 *
677		 * Let w = 2n/x and h = 2/x.  Then the above quotient
678		 * is equal to the continued fraction:
679		 *		     1
680		 *	= -----------------------
681		 *			1
682		 *	   w - -----------------
683		 *			  1
684		 * 		w+h - ---------
685		 *			w+2h - ...
686		 *
687		 * To determine how many terms are needed, run the
688		 * recurrence
689		 *
690		 *	Q(0) = w,
691		 *	Q(1) = w(w+h) - 1,
692		 *	Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
693		 *
694		 * Then when Q(k) > 1e4, k is large enough for single
695		 * precision.
696		 */
697/* XXX NOT DONE - rework this */
698		w = (n + n) / x;
699		h = two / x;
700		q0 = w;
701		z = w + h;
702		q1 = w * z - one;
703		k = 1;
704		while (q1 < big) {
705			k++;
706			z += h;
707			temp = z * q1 - q0;
708			q0 = q1;
709			q1 = temp;
710		}
711		m = n + n;
712		t = zero;
713		for (i = (n + k) << 1; i >= m; i -= 2)
714			t = one / ((double)i / x - t);
715		a = t;
716		b = one;
717		/*
718		 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
719		 * hence, if n*(log(2n/x)) > ...
720		 *	single 8.8722839355e+01
721		 *	double 7.09782712893383973096e+02
722		 *	then recurrent value may overflow and the result is
723		 *	likely underflow to zero
724		 */
725		temp = (double)n;
726		temp *= log((two / x) * temp);
727		if (temp < 7.09782712893383973096e+02) {
728			for (i = n - 1; i > 0; i--) {
729				temp = b;
730				b = b * ((double)(i + i) / x) - a;
731				a = temp;
732			}
733		} else {
734			for (i = n - 1; i > 0; i--) {
735				temp = b;
736				b = b * ((double)(i + i) / x) - a;
737				a = temp;
738				if (b > 1.0e100) {
739					a /= b;
740					t /= b;
741					b = one;
742				}
743			}
744		}
745		b = (t * __k_j0f(fx) / b);
746	}
747	f = (float)b;
748#if defined(__i386) && !defined(__amd64)
749	if (rp != fp_extended)
750		(void) __swapRP(rp);
751#endif
752	return ((sgn)? -f : f);
753}
754
755float
756ynf(int n, float fx)
757{
758	double	a, b, temp, x;
759	float	f;
760	int	i, sign, ix;
761#if defined(__i386) && !defined(__amd64)
762	int	rp;
763#endif
764
765	sign = 0;
766	if (n < 0) {
767		n = -n;
768		if (n & 1)
769			sign = 1;
770	}
771	if (n == 0)
772		return (y0f(fx));
773	if (n == 1)
774		return ((sign)? -y1f(fx) : y1f(fx));
775
776	ix = *(int *)&fx;
777	if ((ix & ~0x80000000) > 0x7f800000)	/* nan */
778		return (fx * fx);
779	if (ix <= 0) {				/* zero or negative */
780		if ((ix << 1) == 0)
781			return (-onef / zerof);
782		return (zerof / zerof);
783	}
784	if (ix == 0x7f800000)			/* +inf */
785		return (zerof);
786
787#if defined(__i386) && !defined(__amd64)
788	rp = __swapRP(fp_extended);
789#endif
790	a = __k_y0f(fx);
791	b = __k_y1f(fx);
792	x = (double)fx;
793	for (i = 1; i < n; i++) {
794		temp = b;
795		b *= (double)(i + i) / x;
796		if (b <= -DBL_MAX)
797			break;
798		b -= a;
799		a = temp;
800	}
801	f = (float)b;
802#if defined(__i386) && !defined(__amd64)
803	if (rp != fp_extended)
804		(void) __swapRP(rp);
805#endif
806	return ((sign)? -f : f);
807}
808