#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: e_jnf.c,v 1.7 1999/07/02 15:37:40 simonb Exp $");
#endif
#include "math.h"
#include "math_private.h"
static const float
two = 2.0000000000e+00,
one = 1.0000000000e+00;
static const float zero = 0.0000000000e+00;
#define __ieee754_j0f j0f
#define __ieee754_j1f j1f
#define __ieee754_y0f y0f
#define __ieee754_y1f y1f
#define __ieee754_logf logf
float jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
float a, b, temp, di;
float z, w;
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
if(ix>0x7f800000) return x+x;
if(n<0){
n = -n;
x = -x;
hx ^= 0x80000000;
}
if(n==0) return(__ieee754_j0f(x));
if(n==1) return(__ieee754_j1f(x));
sgn = (n&1)&(hx>>31);
x = fabsf(x);
if(ix==0||ix>=0x7f800000)
b = zero;
else if((float)n<=x) {
a = __ieee754_j0f(x);
b = __ieee754_j1f(x);
for(i=1;i<n;i++){
temp = b;
b = b*((float)(i+i)/x) - a;
a = temp;
}
} else {
if(ix<0x30800000) {
if(n>33)
b = zero;
else {
temp = x*(float)0.5; b = temp;
for (a=one,i=2;i<=n;i++) {
a *= (float)i;
b *= temp;
}
b = b/a;
}
} else {
float t,v;
float q0,q1,h,tmp; int32_t k,m;
w = (n+n)/(float)x; h = (float)2.0/(float)x;
q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1;
while(q1<(float)1.0e9) {
k += 1; z += h;
tmp = z*q1 - q0;
q0 = q1;
q1 = tmp;
}
m = n+n;
for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
a = t;
b = one;
tmp = n;
v = two/x;
tmp = tmp*__ieee754_logf(fabsf(v*tmp));
if(tmp<(float)8.8721679688e+01) {
for(i=n-1,di=(float)(i+i);i>0;i--){
temp = b;
b *= di;
b = b/x - a;
a = temp;
di -= two;
}
} else {
for(i=n-1,di=(float)(i+i);i>0;i--){
temp = b;
b *= di;
b = b/x - a;
a = temp;
di -= two;
if(b>(float)1e10) {
a /= b;
t /= b;
b = one;
}
}
}
b = (t*__ieee754_j0f(x)/b);
}
}
if(sgn==1) return -b; else return b;
}
float ynf(int n, float x)
{
int32_t i,hx,ix,ib;
int32_t sign;
float a, b, temp;
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
if(ix>0x7f800000) return x+x;
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
sign = 1;
if(n<0){
n = -n;
sign = 1 - ((n&1)<<1);
}
if(n==0) return(__ieee754_y0f(x));
if(n==1) return(sign*__ieee754_y1f(x));
if(ix==0x7f800000) return zero;
a = __ieee754_y0f(x);
b = __ieee754_y1f(x);
GET_FLOAT_WORD(ib,b);
for(i=1;i<n&&(unsigned long)ib!=0xff800000ul;i++){
temp = b;
b = ((float)(i+i)/x)*b - a;
GET_FLOAT_WORD(ib,b);
a = temp;
}
if(sign>0) return b; else return -b;
}