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 118870Srgrimes * software is freely granted, provided that this notice 122116Sjkh * is preserved. 132116Sjkh * ==================================================== 142116Sjkh */ 152116Sjkh 16176451Sdas#include <sys/cdefs.h> 17176451Sdas__FBSDID("$FreeBSD: releng/11.0/lib/msun/src/e_jnf.c 279856 2015-03-10 17:10:54Z kargl $"); 182116Sjkh 19279856Skargl/* 20279856Skargl * See e_jn.c for complete comments. 21279856Skargl */ 22279856Skargl 232116Sjkh#include "math.h" 242116Sjkh#include "math_private.h" 252116Sjkh 26279856Skarglstatic const volatile float vone = 1, vzero = 0; 27279856Skargl 282116Sjkhstatic const float 292116Sjkhtwo = 2.0000000000e+00, /* 0x40000000 */ 302116Sjkhone = 1.0000000000e+00; /* 0x3F800000 */ 312116Sjkh 322116Sjkhstatic const float zero = 0.0000000000e+00; 332116Sjkh 3497413Salfredfloat 3597413Salfred__ieee754_jnf(int n, float x) 362116Sjkh{ 372116Sjkh int32_t i,hx,ix, sgn; 382116Sjkh float a, b, temp, di; 392116Sjkh float z, w; 402116Sjkh 412116Sjkh /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) 422116Sjkh * Thus, J(-n,x) = J(n,-x) 432116Sjkh */ 442116Sjkh GET_FLOAT_WORD(hx,x); 452116Sjkh ix = 0x7fffffff&hx; 462116Sjkh /* if J(n,NaN) is NaN */ 472116Sjkh if(ix>0x7f800000) return x+x; 488870Srgrimes if(n<0){ 492116Sjkh n = -n; 502116Sjkh x = -x; 512116Sjkh hx ^= 0x80000000; 522116Sjkh } 532116Sjkh if(n==0) return(__ieee754_j0f(x)); 542116Sjkh if(n==1) return(__ieee754_j1f(x)); 552116Sjkh sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ 562116Sjkh x = fabsf(x); 572116Sjkh if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */ 582116Sjkh b = zero; 598870Srgrimes else if((float)n<=x) { 602116Sjkh /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ 612116Sjkh a = __ieee754_j0f(x); 622116Sjkh b = __ieee754_j1f(x); 632116Sjkh for(i=1;i<n;i++){ 642116Sjkh temp = b; 652116Sjkh b = b*((float)(i+i)/x) - a; /* avoid underflow */ 662116Sjkh a = temp; 672116Sjkh } 682116Sjkh } else { 692116Sjkh if(ix<0x30800000) { /* x < 2**-29 */ 708870Srgrimes /* x is tiny, return the first Taylor expansion of J(n,x) 712116Sjkh * J(n,x) = 1/n!*(x/2)^n - ... 722116Sjkh */ 732116Sjkh if(n>33) /* underflow */ 742116Sjkh b = zero; 752116Sjkh else { 762116Sjkh temp = x*(float)0.5; b = temp; 772116Sjkh for (a=one,i=2;i<=n;i++) { 782116Sjkh a *= (float)i; /* a = n! */ 792116Sjkh b *= temp; /* b = (x/2)^n */ 802116Sjkh } 812116Sjkh b = b/a; 822116Sjkh } 832116Sjkh } else { 842116Sjkh /* use backward recurrence */ 858870Srgrimes /* x x^2 x^2 862116Sjkh * J(n,x)/J(n-1,x) = ---- ------ ------ ..... 872116Sjkh * 2n - 2(n+1) - 2(n+2) 882116Sjkh * 898870Srgrimes * 1 1 1 902116Sjkh * (for large x) = ---- ------ ------ ..... 912116Sjkh * 2n 2(n+1) 2(n+2) 928870Srgrimes * -- - ------ - ------ - 932116Sjkh * x x x 942116Sjkh * 952116Sjkh * Let w = 2n/x and h=2/x, then the above quotient 962116Sjkh * is equal to the continued fraction: 972116Sjkh * 1 982116Sjkh * = ----------------------- 992116Sjkh * 1 1002116Sjkh * w - ----------------- 1012116Sjkh * 1 1022116Sjkh * w+h - --------- 1032116Sjkh * w+2h - ... 1042116Sjkh * 1052116Sjkh * To determine how many terms needed, let 1062116Sjkh * Q(0) = w, Q(1) = w(w+h) - 1, 1072116Sjkh * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), 1088870Srgrimes * When Q(k) > 1e4 good for single 1098870Srgrimes * When Q(k) > 1e9 good for double 1108870Srgrimes * When Q(k) > 1e17 good for quadruple 1112116Sjkh */ 1122116Sjkh /* determine k */ 1132116Sjkh float t,v; 1142116Sjkh float q0,q1,h,tmp; int32_t k,m; 1152116Sjkh w = (n+n)/(float)x; h = (float)2.0/(float)x; 1162116Sjkh q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1; 1172116Sjkh while(q1<(float)1.0e9) { 1182116Sjkh k += 1; z += h; 1192116Sjkh tmp = z*q1 - q0; 1202116Sjkh q0 = q1; 1212116Sjkh q1 = tmp; 1222116Sjkh } 1232116Sjkh m = n+n; 1242116Sjkh for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); 1252116Sjkh a = t; 1262116Sjkh b = one; 1272116Sjkh /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) 1282116Sjkh * Hence, if n*(log(2n/x)) > ... 1292116Sjkh * single 8.8722839355e+01 1302116Sjkh * double 7.09782712893383973096e+02 1312116Sjkh * long double 1.1356523406294143949491931077970765006170e+04 1328870Srgrimes * then recurrent value may overflow and the result is 1332116Sjkh * likely underflow to zero 1342116Sjkh */ 1352116Sjkh tmp = n; 1362116Sjkh v = two/x; 1372116Sjkh tmp = tmp*__ieee754_logf(fabsf(v*tmp)); 1382116Sjkh if(tmp<(float)8.8721679688e+01) { 1392116Sjkh for(i=n-1,di=(float)(i+i);i>0;i--){ 1402116Sjkh temp = b; 1412116Sjkh b *= di; 1422116Sjkh b = b/x - a; 1432116Sjkh a = temp; 1442116Sjkh di -= two; 1452116Sjkh } 1462116Sjkh } else { 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 /* scale b to avoid spurious overflow */ 1542116Sjkh if(b>(float)1e10) { 1552116Sjkh a /= b; 1562116Sjkh t /= b; 1572116Sjkh b = one; 1582116Sjkh } 1592116Sjkh } 1602116Sjkh } 161215237Suqs z = __ieee754_j0f(x); 162215237Suqs w = __ieee754_j1f(x); 163215237Suqs if (fabsf(z) >= fabsf(w)) 164215237Suqs b = (t*z/b); 165215237Suqs else 166215237Suqs b = (t*w/a); 1672116Sjkh } 1682116Sjkh } 1692116Sjkh if(sgn==1) return -b; else return b; 1702116Sjkh} 1712116Sjkh 17297413Salfredfloat 17397413Salfred__ieee754_ynf(int n, float x) 1742116Sjkh{ 1752116Sjkh int32_t i,hx,ix,ib; 1762116Sjkh int32_t sign; 1772116Sjkh float a, b, temp; 1782116Sjkh 1792116Sjkh GET_FLOAT_WORD(hx,x); 1802116Sjkh ix = 0x7fffffff&hx; 1812116Sjkh if(ix>0x7f800000) return x+x; 182279856Skargl if(ix==0) return -one/vzero; 183279856Skargl if(hx<0) return vzero/vzero; 1842116Sjkh sign = 1; 1852116Sjkh if(n<0){ 1862116Sjkh n = -n; 1877658Sbde sign = 1 - ((n&1)<<1); 1882116Sjkh } 1892116Sjkh if(n==0) return(__ieee754_y0f(x)); 1902116Sjkh if(n==1) return(sign*__ieee754_y1f(x)); 1912116Sjkh if(ix==0x7f800000) return zero; 1922116Sjkh 1932116Sjkh a = __ieee754_y0f(x); 1942116Sjkh b = __ieee754_y1f(x); 1952116Sjkh /* quit if b is -inf */ 1962116Sjkh GET_FLOAT_WORD(ib,b); 1978870Srgrimes for(i=1;i<n&&ib!=0xff800000;i++){ 1982116Sjkh temp = b; 1992116Sjkh b = ((float)(i+i)/x)*b - a; 2002116Sjkh GET_FLOAT_WORD(ib,b); 2012116Sjkh a = temp; 2022116Sjkh } 2032116Sjkh if(sign>0) return b; else return -b; 2042116Sjkh} 205