Deleted Added
sdiff udiff text old ( 279240 ) new ( 279491 )
full compact
e_j1f.c (279240) e_j1f.c (279491)
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>
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 279240 2015-02-24 16:45:16Z pfg $");
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 */
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>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
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 }
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<0x32000000) { /* |x|<2**-27 */
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 */
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>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
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 }
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<=0x24800000) { /* x < 2**-54 */
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;}
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>=0x40f71c58){p = pr5; q= ps5;}
231 else if(ix>=0x4036db68){p = pr3; q= ps3;}
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;
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>=0x40200000) {p = qr8; q= qs8;}
326 else if(ix>=0x40f71c58){p = qr5; q= qs5;}
327 else if(ix>=0x4036db68){p = qr3; q= qs3;}
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}
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}