diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/x2y2m1.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/x2y2m1.c | 24 |
1 files changed, 9 insertions, 15 deletions
diff --git a/sysdeps/ieee754/dbl-64/x2y2m1.c b/sysdeps/ieee754/dbl-64/x2y2m1.c index c96dae59e1..96078888a7 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2015 Free Software Foundation, Inc. + Copyright (C) 2012-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -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]; } |