125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis/*
2225c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2325c28e8Piotr Jasiukajtis */
2425c28e8Piotr Jasiukajtis/*
2525c28e8Piotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
2625c28e8Piotr Jasiukajtis * Use is subject to license terms.
2725c28e8Piotr Jasiukajtis */
2825c28e8Piotr Jasiukajtis
29ddc0e0bRichard Lowe#pragma weak __j0f = j0f
30ddc0e0bRichard Lowe#pragma weak __j1f = j1f
31ddc0e0bRichard Lowe#pragma weak __jnf = jnf
32ddc0e0bRichard Lowe#pragma weak __y0f = y0f
33ddc0e0bRichard Lowe#pragma weak __y1f = y1f
34ddc0e0bRichard Lowe#pragma weak __ynf = ynf
3525c28e8Piotr Jasiukajtis
3625c28e8Piotr Jasiukajtis#include "libm.h"
3725c28e8Piotr Jasiukajtis#include <float.h>
3825c28e8Piotr Jasiukajtis
3925c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
4025c28e8Piotr Jasiukajtisextern int __swapRP(int);
4125c28e8Piotr Jasiukajtis#endif
4225c28e8Piotr Jasiukajtis
4325c28e8Piotr Jasiukajtisstatic const float
4425c28e8Piotr Jasiukajtis	zerof	= 0.0f,
4525c28e8Piotr Jasiukajtis	onef	= 1.0f;
4625c28e8Piotr Jasiukajtis
4725c28e8Piotr Jasiukajtisstatic const double C[] = {
4825c28e8Piotr Jasiukajtis	0.0,
4925c28e8Piotr Jasiukajtis	-0.125,
5025c28e8Piotr Jasiukajtis	0.25,
5125c28e8Piotr Jasiukajtis	0.375,
5225c28e8Piotr Jasiukajtis	0.5,
5325c28e8Piotr Jasiukajtis	1.0,
5425c28e8Piotr Jasiukajtis	2.0,
5525c28e8Piotr Jasiukajtis	8.0,
5625c28e8Piotr Jasiukajtis	0.5641895835477562869480794515607725858441,	/* 1/sqrt(pi) */
5725c28e8Piotr Jasiukajtis	0.636619772367581343075535053490057448,	/* 2/pi */
5825c28e8Piotr Jasiukajtis	1.0e9,
5925c28e8Piotr Jasiukajtis};
6025c28e8Piotr Jasiukajtis
6125c28e8Piotr Jasiukajtis#define	zero	C[0]
6225c28e8Piotr Jasiukajtis#define	neighth	C[1]
6325c28e8Piotr Jasiukajtis#define	quarter	C[2]
6425c28e8Piotr Jasiukajtis#define	three8	C[3]
6525c28e8Piotr Jasiukajtis#define	half	C[4]
6625c28e8Piotr Jasiukajtis#define	one	C[5]
6725c28e8Piotr Jasiukajtis#define	two	C[6]
6825c28e8Piotr Jasiukajtis#define	eight   C[7]
6925c28e8Piotr Jasiukajtis#define	isqrtpi	C[8]
7025c28e8Piotr Jasiukajtis#define	tpi	C[9]
7125c28e8Piotr Jasiukajtis#define	big	C[10]
7225c28e8Piotr Jasiukajtis
7325c28e8Piotr Jasiukajtisstatic const double Cj0y0[] = {
7425c28e8Piotr Jasiukajtis	0.4861344183386052721391238447e5,	/* pr */
7525c28e8Piotr Jasiukajtis	0.1377662549407112278133438945e6,
7625c28e8Piotr Jasiukajtis	0.1222466364088289731869114004e6,
7725c28e8Piotr Jasiukajtis	0.4107070084315176135583353374e5,
7825c28e8Piotr Jasiukajtis	0.5026073801860637125889039915e4,
7925c28e8Piotr Jasiukajtis	0.1783193659125479654541542419e3,
8025c28e8Piotr Jasiukajtis	0.88010344055383421691677564e0,
8125c28e8Piotr Jasiukajtis	0.4861344183386052721414037058e5,	/* ps */
8225c28e8Piotr Jasiukajtis	0.1378196632630384670477582699e6,
8325c28e8Piotr Jasiukajtis	0.1223967185341006542748936787e6,
8425c28e8Piotr Jasiukajtis	0.4120150243795353639995862617e5,
8525c28e8Piotr Jasiukajtis	0.5068271181053546392490184353e4,
8625c28e8Piotr Jasiukajtis	0.1829817905472769960535671664e3,
8725c28e8Piotr Jasiukajtis	1.0,
8825c28e8Piotr Jasiukajtis	-0.1731210995701068539185611951e3,	/* qr */
8925c28e8Piotr Jasiukajtis	-0.5522559165936166961235240613e3,
9025c28e8Piotr Jasiukajtis	-0.5604935606637346590614529613e3,
9125c28e8Piotr Jasiukajtis	-0.2200430300226009379477365011e3,
9225c28e8Piotr Jasiukajtis	-0.323869355375648849771296746e2,
9325c28e8Piotr Jasiukajtis	-0.14294979207907956223499258e1,
9425c28e8Piotr Jasiukajtis	-0.834690374102384988158918e-2,
9525c28e8Piotr Jasiukajtis	0.1107975037248683865326709645e5,	/* qs */
9625c28e8Piotr Jasiukajtis	0.3544581680627082674651471873e5,
9725c28e8Piotr Jasiukajtis	0.3619118937918394132179019059e5,
9825c28e8Piotr Jasiukajtis	0.1439895563565398007471485822e5,
9925c28e8Piotr Jasiukajtis	0.2190277023344363955930226234e4,
10025c28e8Piotr Jasiukajtis	0.106695157020407986137501682e3,
10125c28e8Piotr Jasiukajtis	1.0,
10225c28e8Piotr Jasiukajtis};
10325c28e8Piotr Jasiukajtis
10425c28e8Piotr Jasiukajtis#define	pr	Cj0y0
10525c28e8Piotr Jasiukajtis#define	ps	(Cj0y0+7)
10625c28e8Piotr Jasiukajtis#define	qr	(Cj0y0+14)
10725c28e8Piotr Jasiukajtis#define	qs	(Cj0y0+21)
10825c28e8Piotr Jasiukajtis
10925c28e8Piotr Jasiukajtisstatic const double Cj0[] = {
11025c28e8Piotr Jasiukajtis	-2.500000000000003622131880894830476755537e-0001,	/* r0 */
11125c28e8Piotr Jasiukajtis	1.095597547334830263234433855932375353303e-0002,
11225c28e8Piotr Jasiukajtis	-1.819734750463320921799187258987098087697e-0004,
11325c28e8Piotr Jasiukajtis	9.977001946806131657544212501069893930846e-0007,
11425c28e8Piotr Jasiukajtis	1.0,							/* s0 */
11525c28e8Piotr Jasiukajtis	1.867609810662950169966782360588199673741e-0002,
11625c28e8Piotr Jasiukajtis	1.590389206181565490878430827706972074208e-0004,
11725c28e8Piotr Jasiukajtis	6.520867386742583632375520147714499522721e-0007,
11825c28e8Piotr Jasiukajtis	9.999999999999999942156495584397047660949e-0001,	/* r1 */
11925c28e8Piotr Jasiukajtis	-2.389887722731319130476839836908143731281e-0001,
12025c28e8Piotr Jasiukajtis	1.293359476138939027791270393439493640570e-0002,
12125c28e8Piotr Jasiukajtis	-2.770985642343140122168852400228563364082e-0004,
12225c28e8Piotr Jasiukajtis	2.905241575772067678086738389169625218912e-0006,
12325c28e8Piotr Jasiukajtis	-1.636846356264052597969042009265043251279e-0008,
12425c28e8Piotr Jasiukajtis	5.072306160724884775085431059052611737827e-0011,
12525c28e8Piotr Jasiukajtis	-8.187060730684066824228914775146536139112e-0014,
12625c28e8Piotr Jasiukajtis	5.422219326959949863954297860723723423842e-0017,
12725c28e8Piotr Jasiukajtis	1.0,							/* s1 */
12825c28e8Piotr Jasiukajtis	1.101122772686807702762104741932076228349e-0002,
12925c28e8Piotr Jasiukajtis	6.140169310641649223411427764669143978228e-0005,
13025c28e8Piotr Jasiukajtis	2.292035877515152097976946119293215705250e-0007,
13125c28e8Piotr Jasiukajtis	6.356910426504644334558832036362219583789e-0010,
13225c28e8Piotr Jasiukajtis	1.366626326900219555045096999553948891401e-0012,
13325c28e8Piotr Jasiukajtis	2.280399586866739522891837985560481180088e-0015,
13425c28e8Piotr Jasiukajtis	2.801559820648939665270492520004836611187e-0018,
13525c28e8Piotr Jasiukajtis	2.073101088320349159764410261466350732968e-0021,
13625c28e8Piotr Jasiukajtis};
13725c28e8Piotr Jasiukajtis
13825c28e8Piotr Jasiukajtis#define	r0	Cj0
13925c28e8Piotr Jasiukajtis#define	s0	(Cj0+4)
14025c28e8Piotr Jasiukajtis#define	r1	(Cj0+8)
14125c28e8Piotr Jasiukajtis#define	s1	(Cj0+17)
14225c28e8Piotr Jasiukajtis
14325c28e8Piotr Jasiukajtisstatic const double Cy0[] = {
14425c28e8Piotr Jasiukajtis	-7.380429510868722526754723020704317641941e-0002,	/* u0 */
14525c28e8Piotr Jasiukajtis	1.772607102684869924301459663049874294814e-0001,
14625c28e8Piotr Jasiukajtis	-1.524370666542713828604078090970799356306e-0002,
14725c28e8Piotr Jasiukajtis	4.650819100693891757143771557629924591915e-0004,
14825c28e8Piotr Jasiukajtis	-7.125768872339528975036316108718239946022e-0006,
14925c28e8Piotr Jasiukajtis	6.411017001656104598327565004771515257146e-0008,
15025c28e8Piotr Jasiukajtis	-3.694275157433032553021246812379258781665e-0010,
15125c28e8Piotr Jasiukajtis	1.434364544206266624252820889648445263842e-0012,
15225c28e8Piotr Jasiukajtis	-3.852064731859936455895036286874139896861e-0015,
15325c28e8Piotr Jasiukajtis	7.182052899726138381739945881914874579696e-0018,
15425c28e8Piotr Jasiukajtis	-9.060556574619677567323741194079797987200e-0021,
15525c28e8Piotr Jasiukajtis	7.124435467408860515265552217131230511455e-0024,
15625c28e8Piotr Jasiukajtis	-2.709726774636397615328813121715432044771e-0027,
15725c28e8Piotr Jasiukajtis	1.0,							/* v0 */
15825c28e8Piotr Jasiukajtis	4.678678931512549002587702477349214886475e-0003,
15925c28e8Piotr Jasiukajtis	9.486828955529948534822800829497565178985e-0006,
16025c28e8Piotr Jasiukajtis	1.001495929158861646659010844136682454906e-0008,
16125c28e8Piotr Jasiukajtis	4.725338116256021660204443235685358593611e-0012,
16225c28e8Piotr Jasiukajtis};
16325c28e8Piotr Jasiukajtis
16425c28e8Piotr Jasiukajtis#define	u0	Cy0
16525c28e8Piotr Jasiukajtis#define	v0	(Cy0+13)
16625c28e8Piotr Jasiukajtis
16725c28e8Piotr Jasiukajtisstatic const double Cj1y1[] = {
16825c28e8Piotr Jasiukajtis	-0.4435757816794127857114720794e7,	/* pr0 */
16925c28e8Piotr Jasiukajtis	-0.9942246505077641195658377899e7,
17025c28e8Piotr Jasiukajtis	-0.6603373248364939109255245434e7,
17125c28e8Piotr Jasiukajtis	-0.1523529351181137383255105722e7,
17225c28e8Piotr Jasiukajtis	-0.1098240554345934672737413139e6,
17325c28e8Piotr Jasiukajtis	-0.1611616644324610116477412898e4,
17425c28e8Piotr Jasiukajtis	-0.4435757816794127856828016962e7,	/* ps0 */
17525c28e8Piotr Jasiukajtis	-0.9934124389934585658967556309e7,
17625c28e8Piotr Jasiukajtis	-0.6585339479723087072826915069e7,
17725c28e8Piotr Jasiukajtis	-0.1511809506634160881644546358e7,
17825c28e8Piotr Jasiukajtis	-0.1072638599110382011903063867e6,
17925c28e8Piotr Jasiukajtis	-0.1455009440190496182453565068e4,
18025c28e8Piotr Jasiukajtis	0.3322091340985722351859704442e5,	/* qr0 */
18125c28e8Piotr Jasiukajtis	0.8514516067533570196555001171e5,
18225c28e8Piotr Jasiukajtis	0.6617883658127083517939992166e5,
18325c28e8Piotr Jasiukajtis	0.1849426287322386679652009819e5,
18425c28e8Piotr Jasiukajtis	0.1706375429020768002061283546e4,
18525c28e8Piotr Jasiukajtis	0.3526513384663603218592175580e2,
18625c28e8Piotr Jasiukajtis	0.7087128194102874357377502472e6,	/* qs0 */
18725c28e8Piotr Jasiukajtis	0.1819458042243997298924553839e7,
18825c28e8Piotr Jasiukajtis	0.1419460669603720892855755253e7,
18925c28e8Piotr Jasiukajtis	0.4002944358226697511708610813e6,
19025c28e8Piotr Jasiukajtis	0.3789022974577220264142952256e5,
19125c28e8Piotr Jasiukajtis	0.8638367769604990967475517183e3,
19225c28e8Piotr Jasiukajtis};
19325c28e8Piotr Jasiukajtis
19425c28e8Piotr Jasiukajtis#define	pr0	Cj1y1
19525c28e8Piotr Jasiukajtis#define	ps0	(Cj1y1+6)
19625c28e8Piotr Jasiukajtis#define	qr0	(Cj1y1+12)
19725c28e8Piotr Jasiukajtis#define	qs0	(Cj1y1+18)
19825c28e8Piotr Jasiukajtis
19925c28e8Piotr Jasiukajtisstatic const double Cj1[] = {
20025c28e8Piotr Jasiukajtis	-6.250000000000002203053200981413218949548e-0002,	/* a0 */
20125c28e8Piotr Jasiukajtis	1.600998455640072901321605101981501263762e-0003,
20225c28e8Piotr Jasiukajtis	-1.963888815948313758552511884390162864930e-0005,
20325c28e8Piotr Jasiukajtis	8.263917341093549759781339713418201620998e-0008,
20425c28e8Piotr Jasiukajtis	1.0e0,							/* b0 */
20525c28e8Piotr Jasiukajtis	1.605069137643004242395356851797873766927e-0002,
20625c28e8Piotr Jasiukajtis	1.149454623251299996428500249509098499383e-0004,
20725c28e8Piotr Jasiukajtis	3.849701673735260970379681807910852327825e-0007,
20825c28e8Piotr Jasiukajtis	4.999999999999999995517408894340485471724e-0001,
20925c28e8Piotr Jasiukajtis	-6.003825028120475684835384519945468075423e-0002,
21025c28e8Piotr Jasiukajtis	2.301719899263321828388344461995355419832e-0003,
21125c28e8Piotr Jasiukajtis	-4.208494869238892934859525221654040304068e-0005,
21225c28e8Piotr Jasiukajtis	4.377745135188837783031540029700282443388e-0007,
21325c28e8Piotr Jasiukajtis	-2.854106755678624335145364226735677754179e-0009,
21425c28e8Piotr Jasiukajtis	1.234002865443952024332943901323798413689e-0011,
21525c28e8Piotr Jasiukajtis	-3.645498437039791058951273508838177134310e-0014,
21625c28e8Piotr Jasiukajtis	7.404320596071797459925377103787837414422e-0017,
21725c28e8Piotr Jasiukajtis	-1.009457448277522275262808398517024439084e-0019,
21825c28e8Piotr Jasiukajtis	8.520158355824819796968771418801019930585e-0023,
21925c28e8Piotr Jasiukajtis	-3.458159926081163274483854614601091361424e-0026,
22025c28e8Piotr Jasiukajtis	1.0e0,							/* b1 */
22125c28e8Piotr Jasiukajtis	4.923499437590484879081138588998986303306e-0003,
22225c28e8Piotr Jasiukajtis	1.054389489212184156499666953501976688452e-0005,
22325c28e8Piotr Jasiukajtis	1.180768373106166527048240364872043816050e-0008,
22425c28e8Piotr Jasiukajtis	5.942665743476099355323245707680648588540e-0012,
22525c28e8Piotr Jasiukajtis};
22625c28e8Piotr Jasiukajtis
22725c28e8Piotr Jasiukajtis#define	a0	Cj1
22825c28e8Piotr Jasiukajtis#define	b0	(Cj1+4)
22925c28e8Piotr Jasiukajtis#define	a1	(Cj1+8)
23025c28e8Piotr Jasiukajtis#define	b1	(Cj1+20)
23125c28e8Piotr Jasiukajtis
23225c28e8Piotr Jasiukajtisstatic const double Cy1[] = {
23325c28e8Piotr Jasiukajtis	-1.960570906462389461018983259589655961560e-0001,	/* c0 */
23425c28e8Piotr Jasiukajtis	4.931824118350661953459180060007970291139e-0002,
23525c28e8Piotr Jasiukajtis	-1.626975871565393656845930125424683008677e-0003,
23625c28e8Piotr Jasiukajtis	1.359657517926394132692884168082224258360e-0005,
23725c28e8Piotr Jasiukajtis	1.0e0,							/* d0 */
23825c28e8Piotr Jasiukajtis	2.565807214838390835108224713630901653793e-0002,
23925c28e8Piotr Jasiukajtis	3.374175208978404268650522752520906231508e-0004,
24025c28e8Piotr Jasiukajtis	2.840368571306070719539936935220728843177e-0006,
24125c28e8Piotr Jasiukajtis	1.396387402048998277638900944415752207592e-0008,
24225c28e8Piotr Jasiukajtis	-1.960570906462389473336339614647555351626e-0001,	/* c1 */
24325c28e8Piotr Jasiukajtis	5.336268030335074494231369159933012844735e-0002,
24425c28e8Piotr Jasiukajtis	-2.684137504382748094149184541866332033280e-0003,
24525c28e8Piotr Jasiukajtis	5.737671618979185736981543498580051903060e-0005,
24625c28e8Piotr Jasiukajtis	-6.642696350686335339171171785557663224892e-0007,
24725c28e8Piotr Jasiukajtis	4.692417922568160354012347591960362101664e-0009,
24825c28e8Piotr Jasiukajtis	-2.161728635907789319335231338621412258355e-0011,
24925c28e8Piotr Jasiukajtis	6.727353419738316107197644431844194668702e-0014,
25025c28e8Piotr Jasiukajtis	-1.427502986803861372125234355906790573422e-0016,
25125c28e8Piotr Jasiukajtis	2.020392498726806769468143219616642940371e-0019,
25225c28e8Piotr Jasiukajtis	-1.761371948595104156753045457888272716340e-0022,
25325c28e8Piotr Jasiukajtis	7.352828391941157905175042420249225115816e-0026,
25425c28e8Piotr Jasiukajtis	1.0e0,							/* d1 */
25525c28e8Piotr Jasiukajtis	5.029187436727947764916247076102283399442e-0003,
25625c28e8Piotr Jasiukajtis	1.102693095808242775074856548927801750627e-0005,
25725c28e8Piotr Jasiukajtis	1.268035774543174837829534603830227216291e-0008,
25825c28e8Piotr Jasiukajtis	6.579416271766610825192542295821308730206e-0012,
25925c28e8Piotr Jasiukajtis};
26025c28e8Piotr Jasiukajtis
26125c28e8Piotr Jasiukajtis#define	c0	Cy1
26225c28e8Piotr Jasiukajtis#define	d0	(Cy1+4)
26325c28e8Piotr Jasiukajtis#define	c1	(Cy1+9)
26425c28e8Piotr Jasiukajtis#define	d1	(Cy1+21)
26525c28e8Piotr Jasiukajtis
26625c28e8Piotr Jasiukajtis
26725c28e8Piotr Jasiukajtis/* core of j0f computation; assumes fx is finite */
26825c28e8Piotr Jasiukajtisstatic double
26925c28e8Piotr Jasiukajtis__k_j0f(float fx)
27025c28e8Piotr Jasiukajtis{
27125c28e8Piotr Jasiukajtis	double	x, z, s, c, ss, cc, r, t, p0, q0;
27225c28e8Piotr Jasiukajtis	int	ix, i;
27325c28e8Piotr Jasiukajtis
27425c28e8Piotr Jasiukajtis	ix = *(int *)&fx & ~0x80000000;
27525c28e8Piotr Jasiukajtis	x = fabs((double)fx);
27625c28e8Piotr Jasiukajtis	if (ix > 0x41000000) {
27725c28e8Piotr Jasiukajtis		/* x > 8; see comments in j0.c */
27825c28e8Piotr Jasiukajtis		s = sin(x);
27925c28e8Piotr Jasiukajtis		c = cos(x);
28025c28e8Piotr Jasiukajtis		if (signbit(s) != signbit(c)) {
28125c28e8Piotr Jasiukajtis			ss = s - c;
28225c28e8Piotr Jasiukajtis			cc = -cos(x + x) / ss;
28325c28e8Piotr Jasiukajtis		} else {
28425c28e8Piotr Jasiukajtis			cc = s + c;
28525c28e8Piotr Jasiukajtis			ss = -cos(x + x) / cc;
28625c28e8Piotr Jasiukajtis		}
28725c28e8Piotr Jasiukajtis		if (ix > 0x501502f9) {
28825c28e8Piotr Jasiukajtis			/* x > 1.0e10 */
28925c28e8Piotr Jasiukajtis			p0 = one;
29025c28e8Piotr Jasiukajtis			q0 = neighth / x;
29125c28e8Piotr Jasiukajtis		} else {
29225c28e8Piotr Jasiukajtis			t = eight / x;
29325c28e8Piotr Jasiukajtis			z = t * t;
29425c28e8Piotr Jasiukajtis			p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
29525c28e8Piotr Jasiukajtis			    z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
29625c28e8Piotr Jasiukajtis			    (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
29725c28e8Piotr Jasiukajtis			    z * (ps[4] + z * (ps[5] + z))))));
29825c28e8Piotr Jasiukajtis			q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
29925c28e8Piotr Jasiukajtis			    z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
30025c28e8Piotr Jasiukajtis			    (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
30125c28e8Piotr Jasiukajtis			    z * (qs[4] + z * (qs[5] + z))))))) * t;
30225c28e8Piotr Jasiukajtis		}
30325c28e8Piotr Jasiukajtis		return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
30425c28e8Piotr Jasiukajtis	}
30525c28e8Piotr Jasiukajtis	if (ix <= 0x3727c5ac) {
30625c28e8Piotr Jasiukajtis		/* x <= 1.0e-5 */
30725c28e8Piotr Jasiukajtis		if (ix <= 0x219392ef) /* x <= 1.0e-18 */
30825c28e8Piotr Jasiukajtis			return (one - x);
30925c28e8Piotr Jasiukajtis		return (one - x * x * quarter);
31025c28e8Piotr Jasiukajtis	}
31125c28e8Piotr Jasiukajtis	z = x * x;
31225c28e8Piotr Jasiukajtis	if (ix <= 0x3fa3d70a) {
31325c28e8Piotr Jasiukajtis		/* x <= 1.28 */
31425c28e8Piotr Jasiukajtis		r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
31525c28e8Piotr Jasiukajtis		s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
31625c28e8Piotr Jasiukajtis		return (one + z * (r / s));
31725c28e8Piotr Jasiukajtis	}
31825c28e8Piotr Jasiukajtis	r = r1[8];
31925c28e8Piotr Jasiukajtis	s = s1[8];
32025c28e8Piotr Jasiukajtis	for (i = 7; i >= 0; i--) {
32125c28e8Piotr Jasiukajtis		r = r * z + r1[i];
32225c28e8Piotr Jasiukajtis		s = s * z + s1[i];
32325c28e8Piotr Jasiukajtis	}
32425c28e8Piotr Jasiukajtis	return (r / s);
32525c28e8Piotr Jasiukajtis}
32625c28e8Piotr Jasiukajtis
32725c28e8Piotr Jasiukajtisfloat
32825c28e8Piotr Jasiukajtisj0f(float fx)
32925c28e8Piotr Jasiukajtis{
33025c28e8Piotr Jasiukajtis	float	f;
33125c28e8Piotr Jasiukajtis	int	ix;
33225c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
33325c28e8Piotr Jasiukajtis	int	rp;
33425c28e8Piotr Jasiukajtis#endif
33525c28e8Piotr Jasiukajtis
33625c28e8Piotr Jasiukajtis	ix = *(int *)&fx & ~0x80000000;
33725c28e8Piotr Jasiukajtis	if (ix >= 0x7f800000) {			/* nan or inf */
33825c28e8Piotr Jasiukajtis		if (ix > 0x7f800000)
33925c28e8Piotr Jasiukajtis			return (fx * fx);
34025c28e8Piotr Jasiukajtis		return (zerof);
34125c28e8Piotr Jasiukajtis	}
34225c28e8Piotr Jasiukajtis
34325c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
34425c28e8Piotr Jasiukajtis	rp = __swapRP(fp_extended);
34525c28e8Piotr Jasiukajtis#endif
34625c28e8Piotr Jasiukajtis	f = (float)__k_j0f(fx);
34725c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
34825c28e8Piotr Jasiukajtis	if (rp != fp_extended)
34925c28e8Piotr Jasiukajtis		(void) __swapRP(rp);
35025c28e8Piotr Jasiukajtis#endif
35125c28e8Piotr Jasiukajtis	return (f);
35225c28e8Piotr Jasiukajtis}
35325c28e8Piotr Jasiukajtis
35425c28e8Piotr Jasiukajtis/* core of y0f computation; assumes fx is finite and positive */
35525c28e8Piotr Jasiukajtisstatic double
35625c28e8Piotr Jasiukajtis__k_y0f(float fx)
35725c28e8Piotr Jasiukajtis{
35825c28e8Piotr Jasiukajtis	double	x, z, s, c, ss, cc, t, p0, q0, u, v;
35925c28e8Piotr Jasiukajtis	int	ix, i;
36025c28e8Piotr Jasiukajtis
36125c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
36225c28e8Piotr Jasiukajtis	x = (double)fx;
36325c28e8Piotr Jasiukajtis	if (ix > 0x41000000) {
36425c28e8Piotr Jasiukajtis		/* x > 8; see comments in j0.c */
36525c28e8Piotr Jasiukajtis		s = sin(x);
36625c28e8Piotr Jasiukajtis		c = cos(x);
36725c28e8Piotr Jasiukajtis		if (signbit(s) != signbit(c)) {
36825c28e8Piotr Jasiukajtis			ss = s - c;
36925c28e8Piotr Jasiukajtis			cc = -cos(x + x) / ss;
37025c28e8Piotr Jasiukajtis		} else {
37125c28e8Piotr Jasiukajtis			cc = s + c;
37225c28e8Piotr Jasiukajtis			ss = -cos(x + x) / cc;
37325c28e8Piotr Jasiukajtis		}
37425c28e8Piotr Jasiukajtis		if (ix > 0x501502f9) {
37525c28e8Piotr Jasiukajtis			/* x > 1.0e10 */
37625c28e8Piotr Jasiukajtis			p0 = one;
37725c28e8Piotr Jasiukajtis			q0 = neighth / x;
37825c28e8Piotr Jasiukajtis		} else {
37925c28e8Piotr Jasiukajtis			t = eight / x;
38025c28e8Piotr Jasiukajtis			z = t * t;
38125c28e8Piotr Jasiukajtis			p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
38225c28e8Piotr Jasiukajtis			    z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
38325c28e8Piotr Jasiukajtis			    (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
38425c28e8Piotr Jasiukajtis			    z * (ps[4] + z * (ps[5] + z))))));
38525c28e8Piotr Jasiukajtis			q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
38625c28e8Piotr Jasiukajtis			    z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
38725c28e8Piotr Jasiukajtis			    (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
38825c28e8Piotr Jasiukajtis			    z * (qs[4] + z * (qs[5] + z))))))) * t;
38925c28e8Piotr Jasiukajtis		}
39025c28e8Piotr Jasiukajtis		return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
39125c28e8Piotr Jasiukajtis	}
39225c28e8Piotr Jasiukajtis	if (ix <= 0x219392ef) /* x <= 1.0e-18 */
39325c28e8Piotr Jasiukajtis		return (u0[0] + tpi * log(x));
39425c28e8Piotr Jasiukajtis	z = x * x;
39525c28e8Piotr Jasiukajtis	u = u0[12];
39625c28e8Piotr Jasiukajtis	for (i = 11; i >= 0; i--)
39725c28e8Piotr Jasiukajtis		u = u * z + u0[i];
39825c28e8Piotr Jasiukajtis	v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
39925c28e8Piotr Jasiukajtis	return (u / v + tpi * (__k_j0f(fx) * log(x)));
40025c28e8Piotr Jasiukajtis}
40125c28e8Piotr Jasiukajtis
40225c28e8Piotr Jasiukajtisfloat
40325c28e8Piotr Jasiukajtisy0f(float fx)
40425c28e8Piotr Jasiukajtis{
40525c28e8Piotr Jasiukajtis	float	f;
40625c28e8Piotr Jasiukajtis	int	ix;
40725c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
40825c28e8Piotr Jasiukajtis	int	rp;
40925c28e8Piotr Jasiukajtis#endif
41025c28e8Piotr Jasiukajtis
41125c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
41225c28e8Piotr Jasiukajtis	if ((ix & ~0x80000000) > 0x7f800000)	/* nan */
41325c28e8Piotr Jasiukajtis		return (fx * fx);
41425c28e8Piotr Jasiukajtis	if (ix <= 0) {				/* zero or negative */
41525c28e8Piotr Jasiukajtis		if ((ix << 1) == 0)
41625c28e8Piotr Jasiukajtis			return (-onef / zerof);
41725c28e8Piotr Jasiukajtis		return (zerof / zerof);
41825c28e8Piotr Jasiukajtis	}
41925c28e8Piotr Jasiukajtis	if (ix == 0x7f800000)			/* +inf */
42025c28e8Piotr Jasiukajtis		return (zerof);
42125c28e8Piotr Jasiukajtis
42225c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
42325c28e8Piotr Jasiukajtis	rp = __swapRP(fp_extended);
42425c28e8Piotr Jasiukajtis#endif
42525c28e8Piotr Jasiukajtis	f = (float)__k_y0f(fx);
42625c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
42725c28e8Piotr Jasiukajtis	if (rp != fp_extended)
42825c28e8Piotr Jasiukajtis		(void) __swapRP(rp);
42925c28e8Piotr Jasiukajtis#endif
43025c28e8Piotr Jasiukajtis	return (f);
43125c28e8Piotr Jasiukajtis}
43225c28e8Piotr Jasiukajtis
43325c28e8Piotr Jasiukajtis/* core of j1f computation; assumes fx is finite */
43425c28e8Piotr Jasiukajtisstatic double
43525c28e8Piotr Jasiukajtis__k_j1f(float fx)
43625c28e8Piotr Jasiukajtis{
43725c28e8Piotr Jasiukajtis	double	x, z, s, c, ss, cc, r, t, p1, q1;
43825c28e8Piotr Jasiukajtis	int	i, ix, sgn;
43925c28e8Piotr Jasiukajtis
44025c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
44125c28e8Piotr Jasiukajtis	sgn = (unsigned)ix >> 31;
44225c28e8Piotr Jasiukajtis	ix &= ~0x80000000;
44325c28e8Piotr Jasiukajtis	x = fabs((double)fx);
44425c28e8Piotr Jasiukajtis	if (ix > 0x41000000) {
44525c28e8Piotr Jasiukajtis		/* x > 8; see comments in j1.c */
44625c28e8Piotr Jasiukajtis		s = sin(x);
44725c28e8Piotr Jasiukajtis		c = cos(x);
44825c28e8Piotr Jasiukajtis		if (signbit(s) != signbit(c)) {
44925c28e8Piotr Jasiukajtis			cc = s - c;
45025c28e8Piotr Jasiukajtis			ss = cos(x + x) / cc;
45125c28e8Piotr Jasiukajtis		} else {
45225c28e8Piotr Jasiukajtis			ss = -s - c;
45325c28e8Piotr Jasiukajtis			cc = cos(x + x) / ss;
45425c28e8Piotr Jasiukajtis		}
45525c28e8Piotr Jasiukajtis		if (ix > 0x501502f9) {
45625c28e8Piotr Jasiukajtis			/* x > 1.0e10 */
45725c28e8Piotr Jasiukajtis			p1 = one;
45825c28e8Piotr Jasiukajtis			q1 = three8 / x;
45925c28e8Piotr Jasiukajtis		} else {
46025c28e8Piotr Jasiukajtis			t = eight / x;
46125c28e8Piotr Jasiukajtis			z = t * t;
46225c28e8Piotr Jasiukajtis			p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
46325c28e8Piotr Jasiukajtis			    (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
46425c28e8Piotr Jasiukajtis			    (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
46525c28e8Piotr Jasiukajtis			    (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
46625c28e8Piotr Jasiukajtis			q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
46725c28e8Piotr Jasiukajtis			    (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
46825c28e8Piotr Jasiukajtis			    (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
46925c28e8Piotr Jasiukajtis			    (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
47025c28e8Piotr Jasiukajtis		}
47125c28e8Piotr Jasiukajtis		t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
47225c28e8Piotr Jasiukajtis		return ((sgn)? -t : t);
47325c28e8Piotr Jasiukajtis	}
47425c28e8Piotr Jasiukajtis	if (ix <= 0x3727c5ac) {
47525c28e8Piotr Jasiukajtis		/* x <= 1.0e-5 */
47625c28e8Piotr Jasiukajtis		if (ix <= 0x219392ef) /* x <= 1.0e-18 */
47725c28e8Piotr Jasiukajtis			t = half * x;
47825c28e8Piotr Jasiukajtis		else
47925c28e8Piotr Jasiukajtis			t = x * (half + neighth * x * x);
48025c28e8Piotr Jasiukajtis		return ((sgn)? -t : t);
48125c28e8Piotr Jasiukajtis	}
48225c28e8Piotr Jasiukajtis	z = x * x;
48325c28e8Piotr Jasiukajtis	if (ix < 0x3fa3d70a) {
48425c28e8Piotr Jasiukajtis		/* x < 1.28 */
48525c28e8Piotr Jasiukajtis		r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
48625c28e8Piotr Jasiukajtis		s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
48725c28e8Piotr Jasiukajtis		t = x * half + x * (z * (r / s));
48825c28e8Piotr Jasiukajtis	} else {
48925c28e8Piotr Jasiukajtis		r = a1[11];
49025c28e8Piotr Jasiukajtis		for (i = 10; i >= 0; i--)
49125c28e8Piotr Jasiukajtis			r = r * z + a1[i];
49225c28e8Piotr Jasiukajtis		s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
49325c28e8Piotr Jasiukajtis		t = x * (r / s);
49425c28e8Piotr Jasiukajtis	}
49525c28e8Piotr Jasiukajtis	return ((sgn)? -t : t);
49625c28e8Piotr Jasiukajtis}
49725c28e8Piotr Jasiukajtis
49825c28e8Piotr Jasiukajtisfloat
49925c28e8Piotr Jasiukajtisj1f(float fx)
50025c28e8Piotr Jasiukajtis{
50125c28e8Piotr Jasiukajtis	float	f;
50225c28e8Piotr Jasiukajtis	int	ix;
50325c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
50425c28e8Piotr Jasiukajtis	int	rp;
50525c28e8Piotr Jasiukajtis#endif
50625c28e8Piotr Jasiukajtis
50725c28e8Piotr Jasiukajtis	ix = *(int *)&fx & ~0x80000000;
50825c28e8Piotr Jasiukajtis	if (ix >= 0x7f800000)			/* nan or inf */
50925c28e8Piotr Jasiukajtis		return (onef / fx);
51025c28e8Piotr Jasiukajtis
51125c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
51225c28e8Piotr Jasiukajtis	rp = __swapRP(fp_extended);
51325c28e8Piotr Jasiukajtis#endif
51425c28e8Piotr Jasiukajtis	f = (float)__k_j1f(fx);
51525c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
51625c28e8Piotr Jasiukajtis	if (rp != fp_extended)
51725c28e8Piotr Jasiukajtis		(void) __swapRP(rp);
51825c28e8Piotr Jasiukajtis#endif
51925c28e8Piotr Jasiukajtis	return (f);
52025c28e8Piotr Jasiukajtis}
52125c28e8Piotr Jasiukajtis
52225c28e8Piotr Jasiukajtis/* core of y1f computation; assumes fx is finite and positive */
52325c28e8Piotr Jasiukajtisstatic double
52425c28e8Piotr Jasiukajtis__k_y1f(float fx)
52525c28e8Piotr Jasiukajtis{
52625c28e8Piotr Jasiukajtis	double	x, z, s, c, ss, cc, u, v, p1, q1, t;
52725c28e8Piotr Jasiukajtis	int	i, ix;
52825c28e8Piotr Jasiukajtis
52925c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
53025c28e8Piotr Jasiukajtis	x = (double)fx;
53125c28e8Piotr Jasiukajtis	if (ix > 0x41000000) {
53225c28e8Piotr Jasiukajtis		/* x > 8; see comments in j1.c */
53325c28e8Piotr Jasiukajtis		s = sin(x);
53425c28e8Piotr Jasiukajtis		c = cos(x);
53525c28e8Piotr Jasiukajtis		if (signbit(s) != signbit(c)) {
53625c28e8Piotr Jasiukajtis			cc = s - c;
53725c28e8Piotr Jasiukajtis			ss = cos(x + x) / cc;
53825c28e8Piotr Jasiukajtis		} else {
53925c28e8Piotr Jasiukajtis			ss = -s - c;
54025c28e8Piotr Jasiukajtis			cc = cos(x + x) / ss;
54125c28e8Piotr Jasiukajtis		}
54225c28e8Piotr Jasiukajtis		if (ix > 0x501502f9) {
54325c28e8Piotr Jasiukajtis			/* x > 1.0e10 */
54425c28e8Piotr Jasiukajtis			p1 = one;
54525c28e8Piotr Jasiukajtis			q1 = three8 / x;
54625c28e8Piotr Jasiukajtis		} else {
54725c28e8Piotr Jasiukajtis			t = eight / x;
54825c28e8Piotr Jasiukajtis			z = t * t;
54925c28e8Piotr Jasiukajtis			p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
55025c28e8Piotr Jasiukajtis			    (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
55125c28e8Piotr Jasiukajtis			    (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
55225c28e8Piotr Jasiukajtis			    (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
55325c28e8Piotr Jasiukajtis			q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
55425c28e8Piotr Jasiukajtis			    (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
55525c28e8Piotr Jasiukajtis			    (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
55625c28e8Piotr Jasiukajtis			    (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
55725c28e8Piotr Jasiukajtis		}
55825c28e8Piotr Jasiukajtis		return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
55925c28e8Piotr Jasiukajtis	}
56025c28e8Piotr Jasiukajtis	if (ix <= 0x219392ef) /* x <= 1.0e-18 */
56125c28e8Piotr Jasiukajtis		return (-tpi / x);
56225c28e8Piotr Jasiukajtis	z = x * x;
56325c28e8Piotr Jasiukajtis	if (ix < 0x3fa3d70a) {
56425c28e8Piotr Jasiukajtis		/* x < 1.28 */
56525c28e8Piotr Jasiukajtis		u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
56625c28e8Piotr Jasiukajtis		v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
56725c28e8Piotr Jasiukajtis	} else {
56825c28e8Piotr Jasiukajtis		u = c1[11];
56925c28e8Piotr Jasiukajtis		for (i = 10; i >= 0; i--)
57025c28e8Piotr Jasiukajtis			u = u * z + c1[i];
57125c28e8Piotr Jasiukajtis		v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
57225c28e8Piotr Jasiukajtis	}
57325c28e8Piotr Jasiukajtis	return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
57425c28e8Piotr Jasiukajtis}
57525c28e8Piotr Jasiukajtis
57625c28e8Piotr Jasiukajtisfloat
57725c28e8Piotr Jasiukajtisy1f(float fx)
57825c28e8Piotr Jasiukajtis{
57925c28e8Piotr Jasiukajtis	float	f;
58025c28e8Piotr Jasiukajtis	int	ix;
58125c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
58225c28e8Piotr Jasiukajtis	int	rp;
58325c28e8Piotr Jasiukajtis#endif
58425c28e8Piotr Jasiukajtis
58525c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
58625c28e8Piotr Jasiukajtis	if ((ix & ~0x80000000) > 0x7f800000)	/* nan */
58725c28e8Piotr Jasiukajtis		return (fx * fx);
58825c28e8Piotr Jasiukajtis	if (ix <= 0) {				/* zero or negative */
58925c28e8Piotr Jasiukajtis		if ((ix << 1) == 0)
59025c28e8Piotr Jasiukajtis			return (-onef / zerof);
59125c28e8Piotr Jasiukajtis		return (zerof / zerof);
59225c28e8Piotr Jasiukajtis	}
59325c28e8Piotr Jasiukajtis	if (ix == 0x7f800000)			/* +inf */
59425c28e8Piotr Jasiukajtis		return (zerof);
59525c28e8Piotr Jasiukajtis
59625c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
59725c28e8Piotr Jasiukajtis	rp = __swapRP(fp_extended);
59825c28e8Piotr Jasiukajtis#endif
59925c28e8Piotr Jasiukajtis	f = (float)__k_y1f(fx);
60025c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
60125c28e8Piotr Jasiukajtis	if (rp != fp_extended)
60225c28e8Piotr Jasiukajtis		(void) __swapRP(rp);
60325c28e8Piotr Jasiukajtis#endif
60425c28e8Piotr Jasiukajtis	return (f);
60525c28e8Piotr Jasiukajtis}
60625c28e8Piotr Jasiukajtis
60725c28e8Piotr Jasiukajtisfloat
60825c28e8Piotr Jasiukajtisjnf(int n, float fx)
60925c28e8Piotr Jasiukajtis{
61025c28e8Piotr Jasiukajtis	double	a, b, temp, x, z, w, t, q0, q1, h;
61125c28e8Piotr Jasiukajtis	float	f;
61225c28e8Piotr Jasiukajtis	int	i, ix, sgn, m, k;
61325c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
61425c28e8Piotr Jasiukajtis	int	rp;
61525c28e8Piotr Jasiukajtis#endif
61625c28e8Piotr Jasiukajtis
61725c28e8Piotr Jasiukajtis	if (n < 0) {
61825c28e8Piotr Jasiukajtis		n = -n;
61925c28e8Piotr Jasiukajtis		fx = -fx;
62025c28e8Piotr Jasiukajtis	}
62125c28e8Piotr Jasiukajtis	if (n == 0)
62225c28e8Piotr Jasiukajtis		return (j0f(fx));
62325c28e8Piotr Jasiukajtis	if (n == 1)
62425c28e8Piotr Jasiukajtis		return (j1f(fx));
62525c28e8Piotr Jasiukajtis
62625c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
62725c28e8Piotr Jasiukajtis	sgn = (n & 1)? ((unsigned)ix >> 31) : 0;
62825c28e8Piotr Jasiukajtis	ix &= ~0x80000000;
62925c28e8Piotr Jasiukajtis	if (ix >= 0x7f800000) {		/* nan or inf */
63025c28e8Piotr Jasiukajtis		if (ix > 0x7f800000)
63125c28e8Piotr Jasiukajtis			return (fx * fx);
63225c28e8Piotr Jasiukajtis		return ((sgn)? -zerof : zerof);
63325c28e8Piotr Jasiukajtis	}
63425c28e8Piotr Jasiukajtis	if ((ix << 1) == 0)
63525c28e8Piotr Jasiukajtis		return ((sgn)? -zerof : zerof);
63625c28e8Piotr Jasiukajtis
63725c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
63825c28e8Piotr Jasiukajtis	rp = __swapRP(fp_extended);
63925c28e8Piotr Jasiukajtis#endif
64025c28e8Piotr Jasiukajtis	fx = fabsf(fx);
64125c28e8Piotr Jasiukajtis	x = (double)fx;
64225c28e8Piotr Jasiukajtis	if ((double)n <= x) {
64325c28e8Piotr Jasiukajtis		/* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
64425c28e8Piotr Jasiukajtis		a = __k_j0f(fx);
64525c28e8Piotr Jasiukajtis		b = __k_j1f(fx);
64625c28e8Piotr Jasiukajtis		for (i = 1; i < n; i++) {
64725c28e8Piotr Jasiukajtis			temp = b;
64825c28e8Piotr Jasiukajtis			b = b * ((double)(i + i) / x) - a;
64925c28e8Piotr Jasiukajtis			a = temp;
65025c28e8Piotr Jasiukajtis		}
65125c28e8Piotr Jasiukajtis		f = (float)b;
65225c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
65325c28e8Piotr Jasiukajtis		if (rp != fp_extended)
65425c28e8Piotr Jasiukajtis			(void) __swapRP(rp);
65525c28e8Piotr Jasiukajtis#endif
65625c28e8Piotr Jasiukajtis		return ((sgn)? -f : f);
65725c28e8Piotr Jasiukajtis	}
65825c28e8Piotr Jasiukajtis	if (ix < 0x3089705f) {
65925c28e8Piotr Jasiukajtis		/* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
66025c28e8Piotr Jasiukajtis		if (n > 6)
66125c28e8Piotr Jasiukajtis			n = 6;	/* result underflows to zero for n >= 6 */
66225c28e8Piotr Jasiukajtis		b = t = half * x;
66325c28e8Piotr Jasiukajtis		a = one;
66425c28e8Piotr Jasiukajtis		for (i = 2; i <= n; i++) {
66525c28e8Piotr Jasiukajtis			b *= t;
66625c28e8Piotr Jasiukajtis			a *= (double)i;
66725c28e8Piotr Jasiukajtis		}
66825c28e8Piotr Jasiukajtis		b /= a;
66925c28e8Piotr Jasiukajtis	} else {
67025c28e8Piotr Jasiukajtis		/*
67125c28e8Piotr Jasiukajtis		 * Use the backward recurrence:
67225c28e8Piotr Jasiukajtis		 *
67325c28e8Piotr Jasiukajtis		 * 			x      x^2	x^2
67425c28e8Piotr Jasiukajtis		 *  J(n,x)/J(n-1,x) =  ---- - ------ - ------   .....
67525c28e8Piotr Jasiukajtis		 *			2n    2(n+1)   2(n+2)
67625c28e8Piotr Jasiukajtis		 *
67725c28e8Piotr Jasiukajtis		 * Let w = 2n/x and h = 2/x.  Then the above quotient
67825c28e8Piotr Jasiukajtis		 * is equal to the continued fraction:
67925c28e8Piotr Jasiukajtis		 *		     1
68025c28e8Piotr Jasiukajtis		 *	= -----------------------
68125c28e8Piotr Jasiukajtis		 *			1
68225c28e8Piotr Jasiukajtis		 *	   w - -----------------
68325c28e8Piotr Jasiukajtis		 *			  1
68425c28e8Piotr Jasiukajtis		 * 		w+h - ---------
68525c28e8Piotr Jasiukajtis		 *			w+2h - ...
68625c28e8Piotr Jasiukajtis		 *
68725c28e8Piotr Jasiukajtis		 * To determine how many terms are needed, run the
68825c28e8Piotr Jasiukajtis		 * recurrence
68925c28e8Piotr Jasiukajtis		 *
69025c28e8Piotr Jasiukajtis		 *	Q(0) = w,
69125c28e8Piotr Jasiukajtis		 *	Q(1) = w(w+h) - 1,
69225c28e8Piotr Jasiukajtis		 *	Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
69325c28e8Piotr Jasiukajtis		 *
69425c28e8Piotr Jasiukajtis		 * Then when Q(k) > 1e4, k is large enough for single
69525c28e8Piotr Jasiukajtis		 * precision.
69625c28e8Piotr Jasiukajtis		 */
69725c28e8Piotr Jasiukajtis/* XXX NOT DONE - rework this */
69825c28e8Piotr Jasiukajtis		w = (n + n) / x;
69925c28e8Piotr Jasiukajtis		h = two / x;
70025c28e8Piotr Jasiukajtis		q0 = w;
70125c28e8Piotr Jasiukajtis		z = w + h;
70225c28e8Piotr Jasiukajtis		q1 = w * z - one;
70325c28e8Piotr Jasiukajtis		k = 1;
70425c28e8Piotr Jasiukajtis		while (q1 < big) {
70525c28e8Piotr Jasiukajtis			k++;
70625c28e8Piotr Jasiukajtis			z += h;
70725c28e8Piotr Jasiukajtis			temp = z * q1 - q0;
70825c28e8Piotr Jasiukajtis			q0 = q1;
70925c28e8Piotr Jasiukajtis			q1 = temp;
71025c28e8Piotr Jasiukajtis		}
71125c28e8Piotr Jasiukajtis		m = n + n;
71225c28e8Piotr Jasiukajtis		t = zero;
71325c28e8Piotr Jasiukajtis		for (i = (n + k) << 1; i >= m; i -= 2)
71425c28e8Piotr Jasiukajtis			t = one / ((double)i / x - t);
71525c28e8Piotr Jasiukajtis		a = t;
71625c28e8Piotr Jasiukajtis		b = one;
71725c28e8Piotr Jasiukajtis		/*
71825c28e8Piotr Jasiukajtis		 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
71925c28e8Piotr Jasiukajtis		 * hence, if n*(log(2n/x)) > ...
72025c28e8Piotr Jasiukajtis		 *	single 8.8722839355e+01
72125c28e8Piotr Jasiukajtis		 *	double 7.09782712893383973096e+02
72225c28e8Piotr Jasiukajtis		 *	then recurrent value may overflow and the result is
72325c28e8Piotr Jasiukajtis		 *	likely underflow to zero
72425c28e8Piotr Jasiukajtis		 */
72525c28e8Piotr Jasiukajtis		temp = (double)n;
72625c28e8Piotr Jasiukajtis		temp *= log((two / x) * temp);
72725c28e8Piotr Jasiukajtis		if (temp < 7.09782712893383973096e+02) {
72825c28e8Piotr Jasiukajtis			for (i = n - 1; i > 0; i--) {
72925c28e8Piotr Jasiukajtis				temp = b;
73025c28e8Piotr Jasiukajtis				b = b * ((double)(i + i) / x) - a;
73125c28e8Piotr Jasiukajtis				a = temp;
73225c28e8Piotr Jasiukajtis			}
73325c28e8Piotr Jasiukajtis		} else {
73425c28e8Piotr Jasiukajtis			for (i = n - 1; i > 0; i--) {
73525c28e8Piotr Jasiukajtis				temp = b;
73625c28e8Piotr Jasiukajtis				b = b * ((double)(i + i) / x) - a;
73725c28e8Piotr Jasiukajtis				a = temp;
73825c28e8Piotr Jasiukajtis				if (b > 1.0e100) {
73925c28e8Piotr Jasiukajtis					a /= b;
74025c28e8Piotr Jasiukajtis					t /= b;
74125c28e8Piotr Jasiukajtis					b = one;
74225c28e8Piotr Jasiukajtis				}
74325c28e8Piotr Jasiukajtis			}
74425c28e8Piotr Jasiukajtis		}
74525c28e8Piotr Jasiukajtis		b = (t * __k_j0f(fx) / b);
74625c28e8Piotr Jasiukajtis	}
74725c28e8Piotr Jasiukajtis	f = (float)b;
74825c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
74925c28e8Piotr Jasiukajtis	if (rp != fp_extended)
75025c28e8Piotr Jasiukajtis		(void) __swapRP(rp);
75125c28e8Piotr Jasiukajtis#endif
75225c28e8Piotr Jasiukajtis	return ((sgn)? -f : f);
75325c28e8Piotr Jasiukajtis}
75425c28e8Piotr Jasiukajtis
75525c28e8Piotr Jasiukajtisfloat
75625c28e8Piotr Jasiukajtisynf(int n, float fx)
75725c28e8Piotr Jasiukajtis{
75825c28e8Piotr Jasiukajtis	double	a, b, temp, x;
75925c28e8Piotr Jasiukajtis	float	f;
76025c28e8Piotr Jasiukajtis	int	i, sign, ix;
76125c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
76225c28e8Piotr Jasiukajtis	int	rp;
76325c28e8Piotr Jasiukajtis#endif
76425c28e8Piotr Jasiukajtis
76525c28e8Piotr Jasiukajtis	sign = 0;
76625c28e8Piotr Jasiukajtis	if (n < 0) {
76725c28e8Piotr Jasiukajtis		n = -n;
76825c28e8Piotr Jasiukajtis		if (n & 1)
76925c28e8Piotr Jasiukajtis			sign = 1;
77025c28e8Piotr Jasiukajtis	}
77125c28e8Piotr Jasiukajtis	if (n == 0)
77225c28e8Piotr Jasiukajtis		return (y0f(fx));
77325c28e8Piotr Jasiukajtis	if (n == 1)
77425c28e8Piotr Jasiukajtis		return ((sign)? -y1f(fx) : y1f(fx));
77525c28e8Piotr Jasiukajtis
77625c28e8Piotr Jasiukajtis	ix = *(int *)&fx;
77725c28e8Piotr Jasiukajtis	if ((ix & ~0x80000000) > 0x7f800000)	/* nan */
77825c28e8Piotr Jasiukajtis		return (fx * fx);
77925c28e8Piotr Jasiukajtis	if (ix <= 0) {				/* zero or negative */
78025c28e8Piotr Jasiukajtis		if ((ix << 1) == 0)
78125c28e8Piotr Jasiukajtis			return (-onef / zerof);
78225c28e8Piotr Jasiukajtis		return (zerof / zerof);
78325c28e8Piotr Jasiukajtis	}
78425c28e8Piotr Jasiukajtis	if (ix == 0x7f800000)			/* +inf */
78525c28e8Piotr Jasiukajtis		return (zerof);
78625c28e8Piotr Jasiukajtis
78725c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
78825c28e8Piotr Jasiukajtis	rp = __swapRP(fp_extended);
78925c28e8Piotr Jasiukajtis#endif
79025c28e8Piotr Jasiukajtis	a = __k_y0f(fx);
79125c28e8Piotr Jasiukajtis	b = __k_y1f(fx);
79225c28e8Piotr Jasiukajtis	x = (double)fx;
79325c28e8Piotr Jasiukajtis	for (i = 1; i < n; i++) {
79425c28e8Piotr Jasiukajtis		temp = b;
79525c28e8Piotr Jasiukajtis		b *= (double)(i + i) / x;
79625c28e8Piotr Jasiukajtis		if (b <= -DBL_MAX)
79725c28e8Piotr Jasiukajtis			break;
79825c28e8Piotr Jasiukajtis		b -= a;
79925c28e8Piotr Jasiukajtis		a = temp;
80025c28e8Piotr Jasiukajtis	}
80125c28e8Piotr Jasiukajtis	f = (float)b;
80225c28e8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
80325c28e8Piotr Jasiukajtis	if (rp != fp_extended)
80425c28e8Piotr Jasiukajtis		(void) __swapRP(rp);
80525c28e8Piotr Jasiukajtis#endif
80625c28e8Piotr Jasiukajtis	return ((sign)? -f : f);
80725c28e8Piotr Jasiukajtis}
808