1/* e_j1f.c -- float version of e_j1.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_j1f.c 279491 2015-03-01 20:26:03Z kargl $"); |
18 19#include "math.h" 20#include "math_private.h" 21 22static float ponef(float), qonef(float); 23 24static const float 25huge = 1e30, --- 32 unchanged lines hidden (view full) --- 58 z = cosf(y+y); 59 if ((s*c)>zero) cc = z/ss; 60 else ss = z/cc; 61 } 62 /* 63 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) 64 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) 65 */ |
66 if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */ |
67 else { 68 u = ponef(y); v = qonef(y); 69 z = invsqrtpi*(u*cc-v*ss)/sqrtf(y); 70 } 71 if(hx<0) return -z; 72 else return z; 73 } |
74 if(ix<0x39000000) { /* |x|<2**-13 */ |
75 if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */ 76 } 77 z = x*x; 78 r = z*(r00+z*(r01+z*(r02+z*r03))); 79 s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); 80 r *= x; 81 return(x*(float)0.5+r/s); 82} --- 41 unchanged lines hidden (view full) --- 124 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) 125 * = 1/sqrt(2) * (sin(x) - cos(x)) 126 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) 127 * = -1/sqrt(2) * (cos(x) + sin(x)) 128 * To avoid cancellation, use 129 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) 130 * to compute the worse one. 131 */ |
132 if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ |
133 else { 134 u = ponef(x); v = qonef(x); 135 z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); 136 } 137 return z; 138 } |
139 if(ix<=0x33000000) { /* x < 2**-25 */ |
140 return(-tpi/x); 141 } 142 z = x*x; 143 u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); 144 v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); 145 return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x)); 146} 147 --- 74 unchanged lines hidden (view full) --- 222 static float ponef(float x) 223{ 224 const float *p,*q; 225 float z,r,s; 226 int32_t ix; 227 GET_FLOAT_WORD(ix,x); 228 ix &= 0x7fffffff; 229 if(ix>=0x41000000) {p = pr8; q= ps8;} |
230 else if(ix>=0x409173eb){p = pr5; q= ps5;} 231 else if(ix>=0x4036d917){p = pr3; q= ps3;} |
232 else {p = pr2; q= ps2;} /* ix>=0x40000000 */ 233 z = one/(x*x); 234 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 235 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); 236 return one+ r/s; 237} 238 239 --- 77 unchanged lines hidden (view full) --- 317 318 static float qonef(float x) 319{ 320 const float *p,*q; 321 float s,r,z; 322 int32_t ix; 323 GET_FLOAT_WORD(ix,x); 324 ix &= 0x7fffffff; |
325 if(ix>=0x41000000) {p = qr8; q= qs8;} 326 else if(ix>=0x409173eb){p = qr5; q= qs5;} 327 else if(ix>=0x4036d917){p = qr3; q= qs3;} |
328 else {p = qr2; q= qs2;} /* ix>=0x40000000 */ 329 z = one/(x*x); 330 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); 331 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); 332 return ((float).375 + r/s)/x; 333} |