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}
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}