| e_j0f.c (279240) | e_j0f.c (279491) |
|---|---|
| 1/* e_j0f.c -- float version of e_j0.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 */ 4 5/* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16#include <sys/cdefs.h> | 1/* e_j0f.c -- float version of e_j0.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 */ 4 5/* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16#include <sys/cdefs.h> |
| 17__FBSDID("$FreeBSD: head/lib/msun/src/e_j0f.c 279240 2015-02-24 16:45:16Z pfg $"); | 17__FBSDID("$FreeBSD: head/lib/msun/src/e_j0f.c 279491 2015-03-01 20:26:03Z kargl $"); |
| 18 19#include "math.h" 20#include "math_private.h" 21 22static float pzerof(float), qzerof(float); 23 24static const float 25huge = 1e30, --- 31 unchanged lines hidden (view full) --- 57 z = -cosf(x+x); 58 if ((s*c)<zero) cc = z/ss; 59 else ss = z/cc; 60 } 61 /* 62 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 63 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 64 */ | 18 19#include "math.h" 20#include "math_private.h" 21 22static float pzerof(float), qzerof(float); 23 24static const float 25huge = 1e30, --- 31 unchanged lines hidden (view full) --- 57 z = -cosf(x+x); 58 if ((s*c)<zero) cc = z/ss; 59 else ss = z/cc; 60 } 61 /* 62 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 63 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 64 */ |
| 65 if(ix>0x54000000) z = (invsqrtpi*cc)/sqrtf(x); | 65 if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */ |
| 66 else { 67 u = pzerof(x); v = qzerof(x); 68 z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); 69 } 70 return z; 71 } 72 if(ix<0x3b000000) { /* |x| < 2**-9 */ 73 if(huge+x>one) { /* raise inexact if x != 0 */ --- 57 unchanged lines hidden (view full) --- 131 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 132 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 133 */ 134 if(ix<0x7f000000) { /* make sure x+x not overflow */ 135 z = -cosf(x+x); 136 if ((s*c)<zero) cc = z/ss; 137 else ss = z/cc; 138 } | 66 else { 67 u = pzerof(x); v = qzerof(x); 68 z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); 69 } 70 return z; 71 } 72 if(ix<0x3b000000) { /* |x| < 2**-9 */ 73 if(huge+x>one) { /* raise inexact if x != 0 */ --- 57 unchanged lines hidden (view full) --- 131 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) 132 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) 133 */ 134 if(ix<0x7f000000) { /* make sure x+x not overflow */ 135 z = -cosf(x+x); 136 if ((s*c)<zero) cc = z/ss; 137 else ss = z/cc; 138 } |
| 139 if(ix>0x54800000) z = (invsqrtpi*ss)/sqrtf(x); | 139 if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ |
| 140 else { 141 u = pzerof(x); v = qzerof(x); 142 z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); 143 } 144 return z; 145 } | 140 else { 141 u = pzerof(x); v = qzerof(x); 142 z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); 143 } 144 return z; 145 } |
| 146 if(ix<=0x32000000) { /* x < 2**-27 */ | 146 if(ix<=0x39000000) { /* x < 2**-13 */ |
| 147 return(u00 + tpi*__ieee754_logf(x)); 148 } 149 z = x*x; 150 u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); 151 v = one+z*(v01+z*(v02+z*(v03+z*v04))); 152 return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); 153} 154 --- 72 unchanged lines hidden (view full) --- 227 static float pzerof(float x) 228{ 229 const float *p,*q; 230 float z,r,s; 231 int32_t ix; 232 GET_FLOAT_WORD(ix,x); 233 ix &= 0x7fffffff; 234 if(ix>=0x41000000) {p = pR8; q= pS8;} | 147 return(u00 + tpi*__ieee754_logf(x)); 148 } 149 z = x*x; 150 u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); 151 v = one+z*(v01+z*(v02+z*(v03+z*v04))); 152 return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); 153} 154 --- 72 unchanged lines hidden (view full) --- 227 static float pzerof(float x) 228{ 229 const float *p,*q; 230 float z,r,s; 231 int32_t ix; 232 GET_FLOAT_WORD(ix,x); 233 ix &= 0x7fffffff; 234 if(ix>=0x41000000) {p = pR8; q= pS8;} |
| 235 else if(ix>=0x40f71c58){p = pR5; q= pS5;} 236 else if(ix>=0x4036db68){p = pR3; q= pS3;} | 235 else if(ix>=0x409173eb){p = pR5; q= pS5;} 236 else if(ix>=0x4036d917){p = pR3; q= pS3;} |
| 237 else {p = pR2; q= pS2;} /* ix>=0x40000000 */ 238 z = one/(x*x); 239 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 240 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 241 return one+ r/s; 242} 243 244 --- 77 unchanged lines hidden (view full) --- 322 static float qzerof(float x) 323{ 324 const float *p,*q; 325 float s,r,z; 326 int32_t ix; 327 GET_FLOAT_WORD(ix,x); 328 ix &= 0x7fffffff; 329 if(ix>=0x41000000) {p = qR8; q= qS8;} | 237 else {p = pR2; q= pS2;} /* ix>=0x40000000 */ 238 z = one/(x*x); 239 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 240 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 241 return one+ r/s; 242} 243 244 --- 77 unchanged lines hidden (view full) --- 322 static float qzerof(float x) 323{ 324 const float *p,*q; 325 float s,r,z; 326 int32_t ix; 327 GET_FLOAT_WORD(ix,x); 328 ix &= 0x7fffffff; 329 if(ix>=0x41000000) {p = qR8; q= qS8;} |
| 330 else if(ix>=0x40f71c58){p = qR5; q= qS5;} 331 else if(ix>=0x4036db68){p = qR3; q= qS3;} | 330 else if(ix>=0x409173eb){p = qR5; q= qS5;} 331 else if(ix>=0x4036d917){p = qR3; q= qS3;} |
| 332 else {p = qR2; q= qS2;} /* ix>=0x40000000 */ 333 z = one/(x*x); 334 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 335 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 336 return (-(float).125 + r/s)/x; 337} | 332 else {p = qR2; q= qS2;} /* ix>=0x40000000 */ 333 z = one/(x*x); 334 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 335 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 336 return (-(float).125 + r/s)/x; 337} |