e_jnf.c revision 2116
12116Sjkh/* e_jnf.c -- float version of e_jn.c. 22116Sjkh * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 32116Sjkh */ 42116Sjkh 52116Sjkh/* 62116Sjkh * ==================================================== 72116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 82116Sjkh * 92116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 102116Sjkh * Permission to use, copy, modify, and distribute this 112116Sjkh * software is freely granted, provided that this notice 122116Sjkh * is preserved. 132116Sjkh * ==================================================== 142116Sjkh */ 152116Sjkh 162116Sjkh#ifndef lint 172116Sjkhstatic char rcsid[] = "$Id: e_jnf.c,v 1.2 1994/08/18 23:05:39 jtc Exp $"; 182116Sjkh#endif 192116Sjkh 202116Sjkh#include "math.h" 212116Sjkh#include "math_private.h" 222116Sjkh 232116Sjkh#ifdef __STDC__ 242116Sjkhstatic const float 252116Sjkh#else 262116Sjkhstatic float 272116Sjkh#endif 282116Sjkhinvsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ 292116Sjkhtwo = 2.0000000000e+00, /* 0x40000000 */ 302116Sjkhone = 1.0000000000e+00; /* 0x3F800000 */ 312116Sjkh 322116Sjkh#ifdef __STDC__ 332116Sjkhstatic const float zero = 0.0000000000e+00; 342116Sjkh#else 352116Sjkhstatic float zero = 0.0000000000e+00; 362116Sjkh#endif 372116Sjkh 382116Sjkh#ifdef __STDC__ 392116Sjkh float __ieee754_jnf(int n, float x) 402116Sjkh#else 412116Sjkh float __ieee754_jnf(n,x) 422116Sjkh int n; float x; 432116Sjkh#endif 442116Sjkh{ 452116Sjkh int32_t i,hx,ix, sgn; 462116Sjkh float a, b, temp, di; 472116Sjkh float z, w; 482116Sjkh 492116Sjkh /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) 502116Sjkh * Thus, J(-n,x) = J(n,-x) 512116Sjkh */ 522116Sjkh GET_FLOAT_WORD(hx,x); 532116Sjkh ix = 0x7fffffff&hx; 542116Sjkh /* if J(n,NaN) is NaN */ 552116Sjkh if(ix>0x7f800000) return x+x; 562116Sjkh if(n<0){ 572116Sjkh n = -n; 582116Sjkh x = -x; 592116Sjkh hx ^= 0x80000000; 602116Sjkh } 612116Sjkh if(n==0) return(__ieee754_j0f(x)); 622116Sjkh if(n==1) return(__ieee754_j1f(x)); 632116Sjkh sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ 642116Sjkh x = fabsf(x); 652116Sjkh if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */ 662116Sjkh b = zero; 672116Sjkh else if((float)n<=x) { 682116Sjkh /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ 692116Sjkh a = __ieee754_j0f(x); 702116Sjkh b = __ieee754_j1f(x); 712116Sjkh for(i=1;i<n;i++){ 722116Sjkh temp = b; 732116Sjkh b = b*((float)(i+i)/x) - a; /* avoid underflow */ 742116Sjkh a = temp; 752116Sjkh } 762116Sjkh } else { 772116Sjkh if(ix<0x30800000) { /* x < 2**-29 */ 782116Sjkh /* x is tiny, return the first Taylor expansion of J(n,x) 792116Sjkh * J(n,x) = 1/n!*(x/2)^n - ... 802116Sjkh */ 812116Sjkh if(n>33) /* underflow */ 822116Sjkh b = zero; 832116Sjkh else { 842116Sjkh temp = x*(float)0.5; b = temp; 852116Sjkh for (a=one,i=2;i<=n;i++) { 862116Sjkh a *= (float)i; /* a = n! */ 872116Sjkh b *= temp; /* b = (x/2)^n */ 882116Sjkh } 892116Sjkh b = b/a; 902116Sjkh } 912116Sjkh } else { 922116Sjkh /* use backward recurrence */ 932116Sjkh /* x x^2 x^2 942116Sjkh * J(n,x)/J(n-1,x) = ---- ------ ------ ..... 952116Sjkh * 2n - 2(n+1) - 2(n+2) 962116Sjkh * 972116Sjkh * 1 1 1 982116Sjkh * (for large x) = ---- ------ ------ ..... 992116Sjkh * 2n 2(n+1) 2(n+2) 1002116Sjkh * -- - ------ - ------ - 1012116Sjkh * x x x 1022116Sjkh * 1032116Sjkh * Let w = 2n/x and h=2/x, then the above quotient 1042116Sjkh * is equal to the continued fraction: 1052116Sjkh * 1 1062116Sjkh * = ----------------------- 1072116Sjkh * 1 1082116Sjkh * w - ----------------- 1092116Sjkh * 1 1102116Sjkh * w+h - --------- 1112116Sjkh * w+2h - ... 1122116Sjkh * 1132116Sjkh * To determine how many terms needed, let 1142116Sjkh * Q(0) = w, Q(1) = w(w+h) - 1, 1152116Sjkh * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), 1162116Sjkh * When Q(k) > 1e4 good for single 1172116Sjkh * When Q(k) > 1e9 good for double 1182116Sjkh * When Q(k) > 1e17 good for quadruple 1192116Sjkh */ 1202116Sjkh /* determine k */ 1212116Sjkh float t,v; 1222116Sjkh float q0,q1,h,tmp; int32_t k,m; 1232116Sjkh w = (n+n)/(float)x; h = (float)2.0/(float)x; 1242116Sjkh q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1; 1252116Sjkh while(q1<(float)1.0e9) { 1262116Sjkh k += 1; z += h; 1272116Sjkh tmp = z*q1 - q0; 1282116Sjkh q0 = q1; 1292116Sjkh q1 = tmp; 1302116Sjkh } 1312116Sjkh m = n+n; 1322116Sjkh for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); 1332116Sjkh a = t; 1342116Sjkh b = one; 1352116Sjkh /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) 1362116Sjkh * Hence, if n*(log(2n/x)) > ... 1372116Sjkh * single 8.8722839355e+01 1382116Sjkh * double 7.09782712893383973096e+02 1392116Sjkh * long double 1.1356523406294143949491931077970765006170e+04 1402116Sjkh * then recurrent value may overflow and the result is 1412116Sjkh * likely underflow to zero 1422116Sjkh */ 1432116Sjkh tmp = n; 1442116Sjkh v = two/x; 1452116Sjkh tmp = tmp*__ieee754_logf(fabsf(v*tmp)); 1462116Sjkh if(tmp<(float)8.8721679688e+01) { 1472116Sjkh for(i=n-1,di=(float)(i+i);i>0;i--){ 1482116Sjkh temp = b; 1492116Sjkh b *= di; 1502116Sjkh b = b/x - a; 1512116Sjkh a = temp; 1522116Sjkh di -= two; 1532116Sjkh } 1542116Sjkh } else { 1552116Sjkh for(i=n-1,di=(float)(i+i);i>0;i--){ 1562116Sjkh temp = b; 1572116Sjkh b *= di; 1582116Sjkh b = b/x - a; 1592116Sjkh a = temp; 1602116Sjkh di -= two; 1612116Sjkh /* scale b to avoid spurious overflow */ 1622116Sjkh if(b>(float)1e10) { 1632116Sjkh a /= b; 1642116Sjkh t /= b; 1652116Sjkh b = one; 1662116Sjkh } 1672116Sjkh } 1682116Sjkh } 1692116Sjkh b = (t*__ieee754_j0f(x)/b); 1702116Sjkh } 1712116Sjkh } 1722116Sjkh if(sgn==1) return -b; else return b; 1732116Sjkh} 1742116Sjkh 1752116Sjkh#ifdef __STDC__ 1762116Sjkh float __ieee754_ynf(int n, float x) 1772116Sjkh#else 1782116Sjkh float __ieee754_ynf(n,x) 1792116Sjkh int n; float x; 1802116Sjkh#endif 1812116Sjkh{ 1822116Sjkh int32_t i,hx,ix,ib; 1832116Sjkh int32_t sign; 1842116Sjkh float a, b, temp; 1852116Sjkh 1862116Sjkh GET_FLOAT_WORD(hx,x); 1872116Sjkh ix = 0x7fffffff&hx; 1882116Sjkh /* if Y(n,NaN) is NaN */ 1892116Sjkh if(ix>0x7f800000) return x+x; 1902116Sjkh if(ix==0) return -one/zero; 1912116Sjkh if(hx<0) return zero/zero; 1922116Sjkh sign = 1; 1932116Sjkh if(n<0){ 1942116Sjkh n = -n; 1952116Sjkh sign = 1 - ((n&1)<<2); 1962116Sjkh } 1972116Sjkh if(n==0) return(__ieee754_y0f(x)); 1982116Sjkh if(n==1) return(sign*__ieee754_y1f(x)); 1992116Sjkh if(ix==0x7f800000) return zero; 2002116Sjkh 2012116Sjkh a = __ieee754_y0f(x); 2022116Sjkh b = __ieee754_y1f(x); 2032116Sjkh /* quit if b is -inf */ 2042116Sjkh GET_FLOAT_WORD(ib,b); 2052116Sjkh for(i=1;i<n&&ib!=0xff800000;i++){ 2062116Sjkh temp = b; 2072116Sjkh b = ((float)(i+i)/x)*b - a; 2082116Sjkh GET_FLOAT_WORD(ib,b); 2092116Sjkh a = temp; 2102116Sjkh } 2112116Sjkh if(sign>0) return b; else return -b; 2122116Sjkh} 213