aboutsummaryrefslogtreecommitdiff
path: root/sysdeps
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps')
-rw-r--r--sysdeps/generic/math_private.h4
-rw-r--r--sysdeps/i386/fpu/libm-test-ulps80
-rw-r--r--sysdeps/ieee754/dbl-64/x2y2m1.c22
-rw-r--r--sysdeps/ieee754/dbl-64/x2y2m1f.c4
-rw-r--r--sysdeps/ieee754/ldbl-128/x2y2m1l.c22
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c24
-rw-r--r--sysdeps/ieee754/ldbl-96/x2y2m1.c4
-rw-r--r--sysdeps/ieee754/ldbl-96/x2y2m1l.c22
-rw-r--r--sysdeps/x86_64/fpu/libm-test-ulps80
9 files changed, 119 insertions, 143 deletions
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index a8f1a8e..cf1865d 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -365,8 +365,8 @@ extern double __slowpow (double __x, double __y, double __z);
extern void __docos (double __x, double __dx, double __v[]);
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
extern float __x2y2m1f (float x, float y);
extern double __x2y2m1 (double x, double y);
extern long double __x2y2m1l (long double x, long double y);
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 32f24d0..438a390 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -836,12 +836,12 @@ ildouble: 3
ldouble: 3
Function: Real part of "clog":
-double: 3
-float: 2
-idouble: 3
-ifloat: 2
-ildouble: 4
-ldouble: 4
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 3
+ldouble: 3
Function: Imaginary part of "clog":
double: 1
@@ -864,10 +864,10 @@ ildouble: 2
ldouble: 2
Function: Real part of "clog10_downward":
-double: 5
-float: 4
-idouble: 5
-ifloat: 4
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
ildouble: 8
ldouble: 8
@@ -876,14 +876,14 @@ double: 1
float: 1
idouble: 1
ifloat: 1
-ildouble: 2
-ldouble: 2
+ildouble: 3
+ldouble: 3
Function: Real part of "clog10_towardzero":
-double: 5
-float: 4
-idouble: 5
-ifloat: 4
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
ildouble: 8
ldouble: 8
@@ -896,12 +896,12 @@ ildouble: 3
ldouble: 3
Function: Real part of "clog10_upward":
-double: 5
-float: 5
-idouble: 5
-ifloat: 5
-ildouble: 6
-ldouble: 6
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
+ildouble: 7
+ldouble: 7
Function: Imaginary part of "clog10_upward":
double: 1
@@ -912,12 +912,12 @@ ildouble: 3
ldouble: 3
Function: Real part of "clog_downward":
-double: 5
-float: 5
-idouble: 5
-ifloat: 5
-ildouble: 7
-ldouble: 7
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
+ildouble: 5
+ldouble: 5
Function: Imaginary part of "clog_downward":
double: 1
@@ -928,12 +928,12 @@ ildouble: 1
ldouble: 1
Function: Real part of "clog_towardzero":
-double: 5
-float: 5
-idouble: 5
-ifloat: 5
-ildouble: 8
-ldouble: 8
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
+ildouble: 5
+ldouble: 5
Function: Imaginary part of "clog_towardzero":
double: 1
@@ -944,12 +944,12 @@ ildouble: 1
ldouble: 1
Function: Real part of "clog_upward":
-double: 5
-float: 5
-idouble: 5
-ifloat: 5
-ildouble: 6
-ldouble: 6
+double: 2
+float: 3
+idouble: 2
+ifloat: 3
+ildouble: 4
+ldouble: 4
Function: Imaginary part of "clog_upward":
double: 1
diff --git a/sysdeps/ieee754/dbl-64/x2y2m1.c b/sysdeps/ieee754/dbl-64/x2y2m1.c
index c96dae5..b040097 100644
--- a/sysdeps/ieee754/dbl-64/x2y2m1.c
+++ b/sysdeps/ieee754/dbl-64/x2y2m1.c
@@ -80,32 +80,26 @@ compare (const void *p, const void *q)
}
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
double
__x2y2m1 (double x, double y)
{
- double vals[4];
+ double vals[5];
SET_RESTORE_ROUND (FE_TONEAREST);
mul_split (&vals[1], &vals[0], x, x);
mul_split (&vals[3], &vals[2], y, y);
- if (x >= 0.75)
- vals[1] -= 1.0;
- else
- {
- vals[1] -= 0.5;
- vals[3] -= 0.5;
- }
- qsort (vals, 4, sizeof (double), compare);
+ vals[4] = -1.0;
+ qsort (vals, 5, sizeof (double), compare);
/* Add up the values so that each element of VALS has absolute value
at most equal to the last set bit of the next nonzero
element. */
- for (size_t i = 0; i <= 2; i++)
+ for (size_t i = 0; i <= 3; i++)
{
add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]);
- qsort (vals + i + 1, 3 - i, sizeof (double), compare);
+ qsort (vals + i + 1, 4 - i, sizeof (double), compare);
}
/* Now any error from this addition will be small. */
- return vals[3] + vals[2] + vals[1] + vals[0];
+ return vals[4] + vals[3] + vals[2] + vals[1] + vals[0];
}
diff --git a/sysdeps/ieee754/dbl-64/x2y2m1f.c b/sysdeps/ieee754/dbl-64/x2y2m1f.c
index 43a8acf..835f6a0 100644
--- a/sysdeps/ieee754/dbl-64/x2y2m1f.c
+++ b/sysdeps/ieee754/dbl-64/x2y2m1f.c
@@ -21,8 +21,8 @@
#include <float.h>
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
float
__x2y2m1f (float x, float y)
diff --git a/sysdeps/ieee754/ldbl-128/x2y2m1l.c b/sysdeps/ieee754/ldbl-128/x2y2m1l.c
index 11757c6..a0498c3 100644
--- a/sysdeps/ieee754/ldbl-128/x2y2m1l.c
+++ b/sysdeps/ieee754/ldbl-128/x2y2m1l.c
@@ -80,32 +80,26 @@ compare (const void *p, const void *q)
}
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
long double
__x2y2m1l (long double x, long double y)
{
- long double vals[4];
+ long double vals[5];
SET_RESTORE_ROUNDL (FE_TONEAREST);
mul_split (&vals[1], &vals[0], x, x);
mul_split (&vals[3], &vals[2], y, y);
- if (x >= 0.75L)
- vals[1] -= 1.0L;
- else
- {
- vals[1] -= 0.5L;
- vals[3] -= 0.5L;
- }
- qsort (vals, 4, sizeof (long double), compare);
+ vals[4] = -1.0L;
+ qsort (vals, 5, sizeof (long double), compare);
/* Add up the values so that each element of VALS has absolute value
at most equal to the last set bit of the next nonzero
element. */
- for (size_t i = 0; i <= 2; i++)
+ for (size_t i = 0; i <= 3; i++)
{
add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]);
- qsort (vals + i + 1, 3 - i, sizeof (long double), compare);
+ qsort (vals + i + 1, 4 - i, sizeof (long double), compare);
}
/* Now any error from this addition will be small. */
- return vals[3] + vals[2] + vals[1] + vals[0];
+ return vals[4] + vals[3] + vals[2] + vals[1] + vals[0];
}
diff --git a/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c b/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c
index a001b58..081fb98 100644
--- a/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c
+++ b/sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c
@@ -80,13 +80,13 @@ compare (const void *p, const void *q)
}
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
long double
__x2y2m1l (long double x, long double y)
{
- double vals[12];
+ double vals[13];
SET_RESTORE_ROUND (FE_TONEAREST);
union ibm_extended_long_double xu, yu;
xu.ld = x;
@@ -105,25 +105,19 @@ __x2y2m1l (long double x, long double y)
vals[8] *= 2.0;
vals[9] *= 2.0;
mul_split (&vals[11], &vals[10], yu.d[1].d, yu.d[1].d);
- if (xu.d[0].d >= 0.75)
- vals[1] -= 1.0;
- else
- {
- vals[1] -= 0.5;
- vals[7] -= 0.5;
- }
- qsort (vals, 12, sizeof (double), compare);
+ vals[12] = -1.0;
+ qsort (vals, 13, sizeof (double), compare);
/* Add up the values so that each element of VALS has absolute value
at most equal to the last set bit of the next nonzero
element. */
- for (size_t i = 0; i <= 10; i++)
+ for (size_t i = 0; i <= 11; i++)
{
add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]);
- qsort (vals + i + 1, 11 - i, sizeof (double), compare);
+ qsort (vals + i + 1, 12 - i, sizeof (double), compare);
}
/* Now any error from this addition will be small. */
- long double retval = (long double) vals[11];
- for (size_t i = 10; i != (size_t) -1; i--)
+ long double retval = (long double) vals[12];
+ for (size_t i = 11; i != (size_t) -1; i--)
retval += (long double) vals[i];
return retval;
}
diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1.c b/sysdeps/ieee754/ldbl-96/x2y2m1.c
index a6cc82c..2f6b0be 100644
--- a/sysdeps/ieee754/ldbl-96/x2y2m1.c
+++ b/sysdeps/ieee754/ldbl-96/x2y2m1.c
@@ -27,8 +27,8 @@
#else
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
double
__x2y2m1 (double x, double y)
diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1l.c b/sysdeps/ieee754/ldbl-96/x2y2m1l.c
index 11757c6..a0498c3 100644
--- a/sysdeps/ieee754/ldbl-96/x2y2m1l.c
+++ b/sysdeps/ieee754/ldbl-96/x2y2m1l.c
@@ -80,32 +80,26 @@ compare (const void *p, const void *q)
}
/* Return X^2 + Y^2 - 1, computed without large cancellation error.
- It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
- 0.75 or Y >= 0.5. */
+ It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >=
+ 0.5. */
long double
__x2y2m1l (long double x, long double y)
{
- long double vals[4];
+ long double vals[5];
SET_RESTORE_ROUNDL (FE_TONEAREST);
mul_split (&vals[1], &vals[0], x, x);
mul_split (&vals[3], &vals[2], y, y);
- if (x >= 0.75L)
- vals[1] -= 1.0L;
- else
- {
- vals[1] -= 0.5L;
- vals[3] -= 0.5L;
- }
- qsort (vals, 4, sizeof (long double), compare);
+ vals[4] = -1.0L;
+ qsort (vals, 5, sizeof (long double), compare);
/* Add up the values so that each element of VALS has absolute value
at most equal to the last set bit of the next nonzero
element. */
- for (size_t i = 0; i <= 2; i++)
+ for (size_t i = 0; i <= 3; i++)
{
add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]);
- qsort (vals + i + 1, 3 - i, sizeof (long double), compare);
+ qsort (vals + i + 1, 4 - i, sizeof (long double), compare);
}
/* Now any error from this addition will be small. */
- return vals[3] + vals[2] + vals[1] + vals[0];
+ return vals[4] + vals[3] + vals[2] + vals[1] + vals[0];
}
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 12c3dd1..ef3ab70 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -869,11 +869,11 @@ ldouble: 3
Function: Real part of "clog":
double: 3
-float: 2
+float: 3
idouble: 3
-ifloat: 2
-ildouble: 4
-ldouble: 4
+ifloat: 3
+ildouble: 3
+ldouble: 3
Function: Imaginary part of "clog":
float: 1
@@ -883,9 +883,9 @@ ldouble: 1
Function: Real part of "clog10":
double: 3
-float: 3
+float: 4
idouble: 3
-ifloat: 3
+ifloat: 4
ildouble: 4
ldouble: 4
@@ -898,10 +898,10 @@ ildouble: 2
ldouble: 2
Function: Real part of "clog10_downward":
-double: 6
-float: 6
-idouble: 6
-ifloat: 6
+double: 5
+float: 4
+idouble: 5
+ifloat: 4
ildouble: 8
ldouble: 8
@@ -910,14 +910,14 @@ double: 2
float: 4
idouble: 2
ifloat: 4
-ildouble: 2
-ldouble: 2
+ildouble: 3
+ldouble: 3
Function: Real part of "clog10_towardzero":
double: 5
-float: 4
+float: 5
idouble: 5
-ifloat: 4
+ifloat: 5
ildouble: 8
ldouble: 8
@@ -930,28 +930,28 @@ ildouble: 3
ldouble: 3
Function: Real part of "clog10_upward":
-double: 8
+double: 6
float: 5
-idouble: 8
+idouble: 6
ifloat: 5
-ildouble: 6
-ldouble: 6
+ildouble: 7
+ldouble: 7
Function: Imaginary part of "clog10_upward":
double: 2
-float: 3
+float: 4
idouble: 2
-ifloat: 3
+ifloat: 4
ildouble: 3
ldouble: 3
Function: Real part of "clog_downward":
-double: 7
-float: 5
-idouble: 7
-ifloat: 5
-ildouble: 7
-ldouble: 7
+double: 4
+float: 3
+idouble: 4
+ifloat: 3
+ildouble: 5
+ldouble: 5
Function: Imaginary part of "clog_downward":
double: 1
@@ -962,28 +962,28 @@ ildouble: 1
ldouble: 1
Function: Real part of "clog_towardzero":
-double: 7
-float: 5
-idouble: 7
-ifloat: 5
-ildouble: 8
-ldouble: 8
+double: 4
+float: 4
+idouble: 4
+ifloat: 4
+ildouble: 5
+ldouble: 5
Function: Imaginary part of "clog_towardzero":
double: 1
-float: 2
+float: 3
idouble: 1
-ifloat: 2
+ifloat: 3
ildouble: 1
ldouble: 1
Function: Real part of "clog_upward":
-double: 8
-float: 5
-idouble: 8
-ifloat: 5
-ildouble: 6
-ldouble: 6
+double: 4
+float: 3
+idouble: 4
+ifloat: 3
+ildouble: 4
+ldouble: 4
Function: Imaginary part of "clog_upward":
double: 1
generated by cgit v1.1 at 2024-10-31 21:06:52 +0000