common/R/besself.c

25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __j0f = j0f
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __j1f = j1f
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __jnf = jnf
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __y0f = y0f
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __y1f = y1f
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __ynf = ynf
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include <float.h>
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern int __swapRP(int);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const float
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    zerof   = 0.0f,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    onef    = 1.0f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double C[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.125,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.25,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.375,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    8.0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.5641895835477562869480794515607725858441, /* 1/sqrt(pi) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.636619772367581343075535053490057448, /* 2/pi */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0e9,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define zero    C[0]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define neighth C[1]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define quarter C[2]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define three8  C[3]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define half    C[4]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define one C[5]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define two C[6]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define eight   C[7]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define isqrtpi C[8]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define tpi C[9]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define big C[10]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double Cj0y0[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.4861344183386052721391238447e5,   /* pr */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1377662549407112278133438945e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1222466364088289731869114004e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.4107070084315176135583353374e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.5026073801860637125889039915e4,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1783193659125479654541542419e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.88010344055383421691677564e0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.4861344183386052721414037058e5,   /* ps */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1378196632630384670477582699e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1223967185341006542748936787e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.4120150243795353639995862617e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.5068271181053546392490184353e4,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1829817905472769960535671664e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1731210995701068539185611951e3,  /* qr */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.5522559165936166961235240613e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.5604935606637346590614529613e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.2200430300226009379477365011e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.323869355375648849771296746e2,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.14294979207907956223499258e1,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.834690374102384988158918e-2,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1107975037248683865326709645e5,   /* qs */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.3544581680627082674651471873e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.3619118937918394132179019059e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1439895563565398007471485822e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.2190277023344363955930226234e4,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.106695157020407986137501682e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define pr  Cj0y0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define ps  (Cj0y0+7)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define qr  (Cj0y0+14)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define qs  (Cj0y0+21)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double Cj0[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.500000000000003622131880894830476755537e-0001,   /* r0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.095597547334830263234433855932375353303e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.819734750463320921799187258987098087697e-0004,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    9.977001946806131657544212501069893930846e-0007,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0,                            /* s0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.867609810662950169966782360588199673741e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.590389206181565490878430827706972074208e-0004,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    6.520867386742583632375520147714499522721e-0007,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    9.999999999999999942156495584397047660949e-0001,    /* r1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.389887722731319130476839836908143731281e-0001,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.293359476138939027791270393439493640570e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.770985642343140122168852400228563364082e-0004,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.905241575772067678086738389169625218912e-0006,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.636846356264052597969042009265043251279e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    5.072306160724884775085431059052611737827e-0011,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -8.187060730684066824228914775146536139112e-0014,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    5.422219326959949863954297860723723423842e-0017,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0,                            /* s1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.101122772686807702762104741932076228349e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    6.140169310641649223411427764669143978228e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.292035877515152097976946119293215705250e-0007,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    6.356910426504644334558832036362219583789e-0010,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.366626326900219555045096999553948891401e-0012,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.280399586866739522891837985560481180088e-0015,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.801559820648939665270492520004836611187e-0018,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.073101088320349159764410261466350732968e-0021,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define r0  Cj0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define s0  (Cj0+4)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define r1  (Cj0+8)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define s1  (Cj0+17)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double Cy0[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -7.380429510868722526754723020704317641941e-0002,   /* u0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.772607102684869924301459663049874294814e-0001,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.524370666542713828604078090970799356306e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.650819100693891757143771557629924591915e-0004,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -7.125768872339528975036316108718239946022e-0006,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    6.411017001656104598327565004771515257146e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -3.694275157433032553021246812379258781665e-0010,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.434364544206266624252820889648445263842e-0012,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -3.852064731859936455895036286874139896861e-0015,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    7.182052899726138381739945881914874579696e-0018,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -9.060556574619677567323741194079797987200e-0021,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    7.124435467408860515265552217131230511455e-0024,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.709726774636397615328813121715432044771e-0027,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0,                            /* v0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.678678931512549002587702477349214886475e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    9.486828955529948534822800829497565178985e-0006,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.001495929158861646659010844136682454906e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.725338116256021660204443235685358593611e-0012,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define u0  Cy0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define v0  (Cy0+13)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double Cj1y1[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.4435757816794127857114720794e7,  /* pr0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.9942246505077641195658377899e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.6603373248364939109255245434e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1523529351181137383255105722e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1098240554345934672737413139e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1611616644324610116477412898e4,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.4435757816794127856828016962e7,  /* ps0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.9934124389934585658967556309e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.6585339479723087072826915069e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1511809506634160881644546358e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1072638599110382011903063867e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -0.1455009440190496182453565068e4,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.3322091340985722351859704442e5,   /* qr0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.8514516067533570196555001171e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.6617883658127083517939992166e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1849426287322386679652009819e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1706375429020768002061283546e4,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.3526513384663603218592175580e2,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.7087128194102874357377502472e6,   /* qs0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1819458042243997298924553839e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.1419460669603720892855755253e7,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.4002944358226697511708610813e6,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.3789022974577220264142952256e5,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    0.8638367769604990967475517183e3,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define pr0 Cj1y1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define ps0 (Cj1y1+6)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define qr0 (Cj1y1+12)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define qs0 (Cj1y1+18)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double Cj1[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -6.250000000000002203053200981413218949548e-0002,   /* a0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.600998455640072901321605101981501263762e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.963888815948313758552511884390162864930e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    8.263917341093549759781339713418201620998e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0e0,                          /* b0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.605069137643004242395356851797873766927e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.149454623251299996428500249509098499383e-0004,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    3.849701673735260970379681807910852327825e-0007,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.999999999999999995517408894340485471724e-0001,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -6.003825028120475684835384519945468075423e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.301719899263321828388344461995355419832e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -4.208494869238892934859525221654040304068e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.377745135188837783031540029700282443388e-0007,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.854106755678624335145364226735677754179e-0009,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.234002865443952024332943901323798413689e-0011,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -3.645498437039791058951273508838177134310e-0014,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    7.404320596071797459925377103787837414422e-0017,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.009457448277522275262808398517024439084e-0019,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    8.520158355824819796968771418801019930585e-0023,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -3.458159926081163274483854614601091361424e-0026,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0e0,                          /* b1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.923499437590484879081138588998986303306e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.054389489212184156499666953501976688452e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.180768373106166527048240364872043816050e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    5.942665743476099355323245707680648588540e-0012,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define a0  Cj1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define b0  (Cj1+4)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define a1  (Cj1+8)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define b1  (Cj1+20)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double Cy1[] = {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.960570906462389461018983259589655961560e-0001,   /* c0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.931824118350661953459180060007970291139e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.626975871565393656845930125424683008677e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.359657517926394132692884168082224258360e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0e0,                          /* d0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.565807214838390835108224713630901653793e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    3.374175208978404268650522752520906231508e-0004,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.840368571306070719539936935220728843177e-0006,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.396387402048998277638900944415752207592e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.960570906462389473336339614647555351626e-0001,   /* c1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    5.336268030335074494231369159933012844735e-0002,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.684137504382748094149184541866332033280e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    5.737671618979185736981543498580051903060e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -6.642696350686335339171171785557663224892e-0007,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    4.692417922568160354012347591960362101664e-0009,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -2.161728635907789319335231338621412258355e-0011,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    6.727353419738316107197644431844194668702e-0014,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.427502986803861372125234355906790573422e-0016,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    2.020392498726806769468143219616642940371e-0019,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    -1.761371948595104156753045457888272716340e-0022,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    7.352828391941157905175042420249225115816e-0026,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.0e0,                          /* d1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    5.029187436727947764916247076102283399442e-0003,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.102693095808242775074856548927801750627e-0005,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    1.268035774543174837829534603830227216291e-0008,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    6.579416271766610825192542295821308730206e-0012,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define c0  Cy1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define d0  (Cy1+4)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define c1  (Cy1+9)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define d1  (Cy1+21)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* core of j0f computation; assumes fx is finite */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__k_j0f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double  x, z, s, c, ss, cc, r, t, p0, q0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int ix, i;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx & ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = fabs((double)fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix > 0x41000000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x > 8; see comments in j0.c */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = sin(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        c = cos(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (signbit(s) != signbit(c)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = s - c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = -cos(x + x) / ss;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = s + c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = -cos(x + x) / cc;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix > 0x501502f9) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            /* x > 1.0e10 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p0 = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q0 = neighth / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = eight / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            z = t * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (ps[4] + z * (ps[5] + z))))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (qs[4] + z * (qs[5] + z))))))) * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0x3727c5ac) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x <= 1.0e-5 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix <= 0x219392ef) /* x <= 1.0e-18 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (one - x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (one - x * x * quarter);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    z = x * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0x3fa3d70a) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x <= 1.28 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (one + z * (r / s));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    r = r1[8];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    s = s1[8];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    for (i = 7; i >= 0; i--) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        r = r * z + r1[i];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = s * z + s1[i];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (r / s);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfloat
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisj0f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    float   f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int rp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx & ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix >= 0x7f800000) {         /* nan or inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix > 0x7f800000)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (fx * fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    rp = __swapRP(fp_extended);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    f = (float)__k_j0f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* core of y0f computation; assumes fx is finite and positive */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__k_y0f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double  x, z, s, c, ss, cc, t, p0, q0, u, v;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int ix, i;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = (double)fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix > 0x41000000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x > 8; see comments in j0.c */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = sin(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        c = cos(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (signbit(s) != signbit(c)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = s - c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = -cos(x + x) / ss;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = s + c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = -cos(x + x) / cc;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix > 0x501502f9) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            /* x > 1.0e10 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p0 = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q0 = neighth / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = eight / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            z = t * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (ps[4] + z * (ps[5] + z))))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                z * (qs[4] + z * (qs[5] + z))))))) * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0x219392ef) /* x <= 1.0e-18 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (u0[0] + tpi * log(x));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    z = x * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    u = u0[12];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    for (i = 11; i >= 0; i--)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        u = u * z + u0[i];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (u / v + tpi * (__k_j0f(fx) * log(x)));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfloat
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisy0f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    float   f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int rp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((ix & ~0x80000000) > 0x7f800000)    /* nan */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (fx * fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0) {              /* zero or negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if ((ix << 1) == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (-onef / zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof / zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix == 0x7f800000)           /* +inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    rp = __swapRP(fp_extended);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    f = (float)__k_y0f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* core of j1f computation; assumes fx is finite */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__k_j1f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double  x, z, s, c, ss, cc, r, t, p1, q1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int i, ix, sgn;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sgn = (unsigned)ix >> 31;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix &= ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = fabs((double)fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix > 0x41000000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x > 8; see comments in j1.c */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = sin(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        c = cos(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (signbit(s) != signbit(c)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = s - c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = cos(x + x) / cc;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = -s - c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = cos(x + x) / ss;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix > 0x501502f9) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            /* x > 1.0e10 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p1 = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q1 = three8 / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = eight / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            z = t * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((sgn)? -t : t);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0x3727c5ac) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x <= 1.0e-5 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix <= 0x219392ef) /* x <= 1.0e-18 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = half * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = x * (half + neighth * x * x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((sgn)? -t : t);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    z = x * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix < 0x3fa3d70a) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x < 1.28 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        t = x * half + x * (z * (r / s));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        r = a1[11];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        for (i = 10; i >= 0; i--)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            r = r * z + a1[i];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        t = x * (r / s);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return ((sgn)? -t : t);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfloat
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisj1f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    float   f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int rp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx & ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix >= 0x7f800000)           /* nan or inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (onef / fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    rp = __swapRP(fp_extended);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    f = (float)__k_j1f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* core of y1f computation; assumes fx is finite and positive */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__k_y1f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double  x, z, s, c, ss, cc, u, v, p1, q1, t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int i, ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = (double)fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix > 0x41000000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x > 8; see comments in j1.c */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = sin(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        c = cos(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (signbit(s) != signbit(c)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = s - c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = cos(x + x) / cc;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ss = -s - c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            cc = cos(x + x) / ss;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix > 0x501502f9) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            /* x > 1.0e10 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p1 = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q1 = three8 / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = eight / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            z = t * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0x219392ef) /* x <= 1.0e-18 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (-tpi / x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    z = x * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix < 0x3fa3d70a) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x < 1.28 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        u = c1[11];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        for (i = 10; i >= 0; i--)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            u = u * z + c1[i];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfloat
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisy1f(float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    float   f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int rp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((ix & ~0x80000000) > 0x7f800000)    /* nan */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (fx * fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0) {              /* zero or negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if ((ix << 1) == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (-onef / zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof / zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix == 0x7f800000)           /* +inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    rp = __swapRP(fp_extended);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    f = (float)__k_y1f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfloat
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisjnf(int n, float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double  a, b, temp, x, z, w, t, q0, q1, h;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    float   f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int i, ix, sgn, m, k;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int rp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (n < 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        n = -n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        fx = -fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (n == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (j0f(fx));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (n == 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (j1f(fx));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sgn = (n & 1)? ((unsigned)ix >> 31) : 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix &= ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix >= 0x7f800000) {     /* nan or inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix > 0x7f800000)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (fx * fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((sgn)? -zerof : zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((ix << 1) == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((sgn)? -zerof : zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    rp = __swapRP(fp_extended);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    fx = fabsf(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = (double)fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((double)n <= x) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        a = __k_j0f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b = __k_j1f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        for (i = 1; i < n; i++) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            temp = b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            b = b * ((double)(i + i) / x) - a;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            a = temp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        f = (float)b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((sgn)? -f : f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix < 0x3089705f) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (n > 6)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            n = 6;  /* result underflows to zero for n >= 6 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b = t = half * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        a = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        for (i = 2; i <= n; i++) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            b *= t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            a *= (double)i;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b /= a;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * Use the backward recurrence:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *          x      x^2  x^2
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  J(n,x)/J(n-1,x) =  ---- - ------ - ------   .....
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *          2n    2(n+1)   2(n+2)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * Let w = 2n/x and h = 2/x.  Then the above quotient
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * is equal to the continued fraction:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *           1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  = -----------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *          1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *     w - -----------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *            1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *      w+h - ---------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *          w+2h - ...
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * To determine how many terms are needed, run the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * recurrence
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  Q(0) = w,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  Q(1) = w(w+h) - 1,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * Then when Q(k) > 1e4, k is large enough for single
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * precision.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* XXX NOT DONE - rework this */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        w = (n + n) / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        h = two / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        q0 = w;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        z = w + h;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        q1 = w * z - one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        k = 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        while (q1 < big) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            k++;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            z += h;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            temp = z * q1 - q0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q0 = q1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            q1 = temp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        m = n + n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        t = zero;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        for (i = (n + k) << 1; i >= m; i -= 2)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = one / ((double)i / x - t);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        a = t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         * hence, if n*(log(2n/x)) > ...
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  single 8.8722839355e+01
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  double 7.09782712893383973096e+02
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  then recurrent value may overflow and the result is
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         *  likely underflow to zero
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis         */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        temp = (double)n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        temp *= log((two / x) * temp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (temp < 7.09782712893383973096e+02) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            for (i = n - 1; i > 0; i--) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                temp = b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                b = b * ((double)(i + i) / x) - a;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                a = temp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            for (i = n - 1; i > 0; i--) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                temp = b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                b = b * ((double)(i + i) / x) - a;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                a = temp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                if (b > 1.0e100) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                    a /= b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                    t /= b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                    b = one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b = (t * __k_j0f(fx) / b);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    f = (float)b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return ((sgn)? -f : f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfloat
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisynf(int n, float fx)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis{
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double  a, b, temp, x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    float   f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int i, sign, ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int rp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sign = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (n < 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        n = -n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (n & 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            sign = 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (n == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (y0f(fx));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (n == 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((sign)? -y1f(fx) : y1f(fx));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = *(int *)&fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((ix & ~0x80000000) > 0x7f800000)    /* nan */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (fx * fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix <= 0) {              /* zero or negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if ((ix << 1) == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (-onef / zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof / zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix == 0x7f800000)           /* +inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (zerof);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    rp = __swapRP(fp_extended);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    a = __k_y0f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    b = __k_y1f(fx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = (double)fx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    for (i = 1; i < n; i++) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        temp = b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b *= (double)(i + i) / x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (b <= -DBL_MAX)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            break;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        b -= a;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        a = temp;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    f = (float)b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__i386) && !defined(__amd64)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (rp != fp_extended)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (void) __swapRP(rp);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return ((sign)? -f : f);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}