1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
From 8cba1dc46c8f29261aa441e70bac798f2c2c0f58 Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Wed, 15 Mar 2017 02:55:49 +0100
Subject: [PATCH] fix threshold constants in j0f, y0f, j1f, y1f
partly following freebsd rev 279491
https://svnweb.freebsd.org/base?view=revision&revision=279491
(musl had some of the fixes before freebsd).
the change should not matter much for j0f, y0f, but it improves
j1f and y1f in [2.5,~3.75] (that is [0x40200000,~0x40700000]).
near roots (e.g. around 3.8317 for j1f) there are still large
ulp errors.
dropped code that tried to raise inexact.
---
src/math/j0f.c | 8 ++++----
src/math/j1f.c | 17 ++++++++---------
2 files changed, 12 insertions(+), 13 deletions(-)
diff --git a/src/math/j0f.c b/src/math/j0f.c
index 45883dc4..fab554a3 100644
--- a/src/math/j0f.c
+++ b/src/math/j0f.c
@@ -208,8 +208,8 @@ static float pzerof(float x)
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = pR8; q = pS8;}
- else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
- else if (ix >= 0x4036db68){p = pR3; q = pS3;}
+ else if (ix >= 0x409173eb){p = pR5; q = pS5;}
+ else if (ix >= 0x4036d917){p = pR3; q = pS3;}
else /*ix >= 0x40000000*/ {p = pR2; q = pS2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
@@ -304,8 +304,8 @@ static float qzerof(float x)
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = qR8; q = qS8;}
- else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
- else if (ix >= 0x4036db68){p = qR3; q = qS3;}
+ else if (ix >= 0x409173eb){p = qR5; q = qS5;}
+ else if (ix >= 0x4036d917){p = qR3; q = qS3;}
else /*ix >= 0x40000000*/ {p = qR2; q = qS2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
diff --git a/src/math/j1f.c b/src/math/j1f.c
index 58875af9..3434c53d 100644
--- a/src/math/j1f.c
+++ b/src/math/j1f.c
@@ -74,14 +74,13 @@ float j1f(float x)
return 1/(x*x);
if (ix >= 0x40000000) /* |x| >= 2 */
return common(ix, fabsf(x), 0, sign);
- if (ix >= 0x32000000) { /* |x| >= 2**-27 */
+ if (ix >= 0x39000000) { /* |x| >= 2**-13 */
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
z = 0.5f + r/s;
} else
- /* raise inexact if x!=0 */
- z = 0.5f + x;
+ z = 0.5f;
return z*x;
}
@@ -114,7 +113,7 @@ float y1f(float x)
return 1/x;
if (ix >= 0x40000000) /* |x| >= 2.0 */
return common(ix,x,1,0);
- if (ix < 0x32000000) /* x < 2**-27 */
+ if (ix < 0x33000000) /* x < 2**-25 */
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
@@ -205,8 +204,8 @@ static float ponef(float x)
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = pr8; q = ps8;}
- else if (ix >= 0x40f71c58){p = pr5; q = ps5;}
- else if (ix >= 0x4036db68){p = pr3; q = ps3;}
+ else if (ix >= 0x409173eb){p = pr5; q = ps5;}
+ else if (ix >= 0x4036d917){p = pr3; q = ps3;}
else /*ix >= 0x40000000*/ {p = pr2; q = ps2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
@@ -300,9 +299,9 @@ static float qonef(float x)
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
- if (ix >= 0x40200000){p = qr8; q = qs8;}
- else if (ix >= 0x40f71c58){p = qr5; q = qs5;}
- else if (ix >= 0x4036db68){p = qr3; q = qs3;}
+ if (ix >= 0x41000000){p = qr8; q = qs8;}
+ else if (ix >= 0x409173eb){p = qr5; q = qs5;}
+ else if (ix >= 0x4036d917){p = qr3; q = qs3;}
else /*ix >= 0x40000000*/ {p = qr2; q = qs2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
--
2.11.1
|
