Skip to content

Commit 4c8740a

Browse files
committed
1 parent c29e8b1 commit 4c8740a

3 files changed

Lines changed: 27 additions & 39 deletions

File tree

src/e_j0f.c

Lines changed: 12 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -14,12 +14,9 @@
1414
*/
1515

1616
#include <assert.h>
17-
1817
#include "cdefs-compat.h"
19-
//__FBSDID("$FreeBSD: src/lib/msun/src/e_j0f.c,v 1.8 2008/02/22 02:30:35 das Exp $");
2018

2119
#include <openlibm_math.h>
22-
2320
#include "math_private.h"
2421

2522
static float pzerof(float), qzerof(float);
@@ -65,17 +62,17 @@ __ieee754_j0f(float x)
6562
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
6663
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
6764
*/
68-
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
65+
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */
6966
else {
7067
u = pzerof(x); v = qzerof(x);
7168
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
7269
}
7370
return z;
7471
}
75-
if(ix<0x39000000) { /* |x| < 2**-13 */
72+
if(ix<0x3b000000) { /* |x| < 2**-9 */
7673
if(huge+x>one) { /* raise inexact if x != 0 */
77-
if(ix<0x32000000) return one; /* |x|<2**-27 */
78-
else return one - (float)0.25*x*x;
74+
if(ix<0x39800000) return one; /* |x|<2**-12 */
75+
else return one - x*x/4;
7976
}
8077
}
8178
z = x*x;
@@ -139,14 +136,14 @@ __ieee754_y0f(float x)
139136
if ((s*c)<zero) cc = z/ss;
140137
else ss = z/cc;
141138
}
142-
if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
139+
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
143140
else {
144141
u = pzerof(x); v = qzerof(x);
145142
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
146143
}
147144
return z;
148145
}
149-
if(ix<=0x32000000) { /* x < 2**-27 */
146+
if(ix<=0x39000000) { /* x < 2**-13 */
150147
return(u00 + tpi*__ieee754_logf(x));
151148
}
152149
z = x*x;
@@ -227,19 +224,17 @@ static const float pS2[5] = {
227224
1.4657617569e+01, /* 0x416a859a */
228225
};
229226

230-
/* Note: This function is only called for ix>=0x40000000 (see above) */
231227
static float pzerof(float x)
232228
{
233229
const float *p,*q;
234230
float z,r,s;
235231
int32_t ix;
236232
GET_FLOAT_WORD(ix,x);
237233
ix &= 0x7fffffff;
238-
assert(ix>=0x40000000 && ix<=0x48000000);
239234
if(ix>=0x41000000) {p = pR8; q= pS8;}
240-
else if(ix>=0x40f71c58){p = pR5; q= pS5;}
241-
else if(ix>=0x4036db68){p = pR3; q= pS3;}
242-
else {p = pR2; q= pS2;}
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 */
243238
z = one/(x*x);
244239
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
245240
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@@ -324,19 +319,17 @@ static const float qS2[6] = {
324319
-5.3109550476e+00, /* 0xc0a9f358 */
325320
};
326321

327-
/* Note: This function is only called for ix>=0x40000000 (see above) */
328322
static float qzerof(float x)
329323
{
330324
const float *p,*q;
331325
float s,r,z;
332326
int32_t ix;
333327
GET_FLOAT_WORD(ix,x);
334328
ix &= 0x7fffffff;
335-
assert(ix>=0x40000000 && ix<=0x48000000);
336329
if(ix>=0x41000000) {p = qR8; q= qS8;}
337-
else if(ix>=0x40f71c58){p = qR5; q= qS5;}
338-
else if(ix>=0x4036db68){p = qR3; q= qS3;}
339-
else {p = qR2; q= qS2;}
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 */
340333
z = one/(x*x);
341334
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
342335
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

src/e_j1f.c

Lines changed: 11 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -14,12 +14,8 @@
1414
*/
1515

1616
#include <assert.h>
17-
1817
#include "cdefs-compat.h"
19-
//__FBSDID("$FreeBSD: src/lib/msun/src/e_j1f.c,v 1.8 2008/02/22 02:30:35 das Exp $");
20-
2118
#include <openlibm_math.h>
22-
2319
#include "math_private.h"
2420

2521
static float ponef(float), qonef(float);
@@ -66,15 +62,15 @@ __ieee754_j1f(float x)
6662
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
6763
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
6864
*/
69-
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
65+
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */
7066
else {
7167
u = ponef(y); v = qonef(y);
7268
z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
7369
}
7470
if(hx<0) return -z;
7571
else return z;
7672
}
77-
if(ix<0x32000000) { /* |x|<2**-27 */
73+
if(ix<0x39000000) { /* |x|<2**-13 */
7874
if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
7975
}
8076
z = x*x;
@@ -132,14 +128,14 @@ __ieee754_y1f(float x)
132128
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
133129
* to compute the worse one.
134130
*/
135-
if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
131+
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
136132
else {
137133
u = ponef(x); v = qonef(x);
138134
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
139135
}
140136
return z;
141137
}
142-
if(ix<=0x24800000) { /* x < 2**-54 */
138+
if(ix<=0x33000000) { /* x < 2**-25 */
143139
return(-tpi/x);
144140
}
145141
z = x*x;
@@ -222,19 +218,17 @@ static const float ps2[5] = {
222218
8.3646392822e+00, /* 0x4105d590 */
223219
};
224220

225-
/* Note: This function is only called for ix>=0x40000000 (see above) */
226221
static float ponef(float x)
227222
{
228223
const float *p,*q;
229224
float z,r,s;
230225
int32_t ix;
231226
GET_FLOAT_WORD(ix,x);
232227
ix &= 0x7fffffff;
233-
assert(ix>=0x40000000 && ix<=0x48000000);
234228
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;}
237-
else {p = pr2; q= ps2;}
229+
else if(ix>=0x409173eb){p = pr5; q= ps5;}
230+
else if(ix>=0x4036d917){p = pr3; q= ps3;}
231+
else {p = pr2; q= ps2;} /* ix>=0x40000000 */
238232
z = one/(x*x);
239233
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
240234
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@@ -320,19 +314,17 @@ static const float qs2[6] = {
320314
-4.9594988823e+00, /* 0xc09eb437 */
321315
};
322316

323-
/* Note: This function is only called for ix>=0x40000000 (see above) */
324317
static float qonef(float x)
325318
{
326319
const float *p,*q;
327320
float s,r,z;
328321
int32_t ix;
329322
GET_FLOAT_WORD(ix,x);
330323
ix &= 0x7fffffff;
331-
assert(ix>=0x40000000 && ix<=0x48000000);
332-
if(ix>=0x40200000) {p = qr8; q= qs8;}
333-
else if(ix>=0x40f71c58){p = qr5; q= qs5;}
334-
else if(ix>=0x4036db68){p = qr3; q= qs3;}
335-
else {p = qr2; q= qs2;}
324+
if(ix>=0x41000000) {p = qr8; q= qs8;}
325+
else if(ix>=0x409173eb){p = qr5; q= qs5;}
326+
else if(ix>=0x4036d917){p = qr3; q= qs3;}
327+
else {p = qr2; q= qs2;} /* ix>=0x40000000 */
336328
z = one/(x*x);
337329
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
338330
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

test/Makefile

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,10 @@ OPENLIBM_HOME=$(abspath ..)
22
include ../Make.inc
33

44
# Set rpath of tests to builddir for loading shared library
5-
OPENLIBM_LIB = -L.. -lopenlibm -Wl,-rpath=$(OPENLIBM_HOME)
5+
OPENLIBM_LIB = -L.. -lopenlibm
6+
ifeq ($(OS),Linux)
7+
OPENLIBM_LIB += -Wl,-rpath=$(OPENLIBM_HOME)
8+
endif
69

710
all: test-double test-float # test-double-system test-float-system
811

0 commit comments

Comments
 (0)