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 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 */
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 }
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 }
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;}
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;}
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}
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 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 */
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 }
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 }
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;}
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;}
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}