diff options
Diffstat (limited to 'sysdeps/ieee754')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/x2y2m1.c | 22 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/x2y2m1f.c | 4 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128/x2y2m1l.c | 22 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c | 24 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/x2y2m1.c | 4 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/x2y2m1l.c | 22 |
6 files changed, 37 insertions, 61 deletions
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]; } |