diff options
Diffstat (limited to 'REORG.TODO/sysdeps/ieee754/dbl-64/e_remainder.c')
| -rw-r--r-- | REORG.TODO/sysdeps/ieee754/dbl-64/e_remainder.c | 152 |
1 files changed, 152 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/ieee754/dbl-64/e_remainder.c b/REORG.TODO/sysdeps/ieee754/dbl-64/e_remainder.c new file mode 100644 index 0000000000..1a2eeed2e1 --- /dev/null +++ b/REORG.TODO/sysdeps/ieee754/dbl-64/e_remainder.c @@ -0,0 +1,152 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001-2017 Free Software Foundation, Inc. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, see <http://www.gnu.org/licenses/>. + */ +/**************************************************************************/ +/* MODULE_NAME urem.c */ +/* */ +/* FUNCTION: uremainder */ +/* */ +/* An ultimate remainder routine. Given two IEEE double machine numbers x */ +/* ,y it computes the correctly rounded (to nearest) value of remainder */ +/* of dividing x by y. */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/* ************************************************************************/ + +#include "endian.h" +#include "mydefs.h" +#include "urem.h" +#include "MathLib.h" +#include <math.h> +#include <math_private.h> + +/**************************************************************************/ +/* An ultimate remainder routine. Given two IEEE double machine numbers x */ +/* ,y it computes the correctly rounded (to nearest) value of remainder */ +/**************************************************************************/ +double +__ieee754_remainder (double x, double y) +{ + double z, d, xx; + int4 kx, ky, n, nn, n1, m1, l; + mynumber u, t, w = { { 0, 0 } }, v = { { 0, 0 } }, ww = { { 0, 0 } }, r; + u.x = x; + t.x = y; + kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign for x*/ + t.i[HIGH_HALF] &= 0x7fffffff; /*no sign for y */ + ky = t.i[HIGH_HALF]; + /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/ + if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000) + { + SET_RESTORE_ROUND_NOEX (FE_TONEAREST); + if (kx + 0x00100000 < ky) + return x; + if ((kx - 0x01500000) < ky) + { + z = x / t.x; + v.i[HIGH_HALF] = t.i[HIGH_HALF]; + d = (z + big.x) - big.x; + xx = (x - d * v.x) - d * (t.x - v.x); + if (d - z != 0.5 && d - z != -0.5) + return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x); + else + { + if (fabs (xx) > 0.5 * t.x) + return (z > d) ? xx - t.x : xx + t.x; + else + return xx; + } + } /* (kx<(ky+0x01500000)) */ + else + { + r.x = 1.0 / t.x; + n = t.i[HIGH_HALF]; + nn = (n & 0x7ff00000) + 0x01400000; + w.i[HIGH_HALF] = n; + ww.x = t.x - w.x; + l = (kx - nn) & 0xfff00000; + n1 = ww.i[HIGH_HALF]; + m1 = r.i[HIGH_HALF]; + while (l > 0) + { + r.i[HIGH_HALF] = m1 - l; + z = u.x * r.x; + w.i[HIGH_HALF] = n + l; + ww.i[HIGH_HALF] = (n1) ? n1 + l : n1; + d = (z + big.x) - big.x; + u.x = (u.x - d * w.x) - d * ww.x; + l = (u.i[HIGH_HALF] & 0x7ff00000) - nn; + } + r.i[HIGH_HALF] = m1; + w.i[HIGH_HALF] = n; + ww.i[HIGH_HALF] = n1; + z = u.x * r.x; + d = (z + big.x) - big.x; + u.x = (u.x - d * w.x) - d * ww.x; + if (fabs (u.x) < 0.5 * t.x) + return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x); + else + if (fabs (u.x) > 0.5 * t.x) + return (d > z) ? u.x + t.x : u.x - t.x; + else + { + z = u.x / t.x; d = (z + big.x) - big.x; + return ((u.x - d * w.x) - d * ww.x); + } + } + } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */ + else + { + if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) + { + y = fabs (y) * t128.x; + z = __ieee754_remainder (x, y) * t128.x; + z = __ieee754_remainder (z, y) * tm128.x; + return z; + } + else + { + if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 && + (ky > 0 || t.i[LOW_HALF] != 0)) + { + y = fabs (y); + z = 2.0 * __ieee754_remainder (0.5 * x, y); + d = fabs (z); + if (d <= fabs (d - y)) + return z; + else if (d == y) + return 0.0 * x; + else + return (z > 0) ? z - y : z + y; + } + else /* if x is too big */ + { + if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */ + return (x * y) / (x * y); + else if (kx >= 0x7ff00000 /* x not finite */ + || (ky > 0x7ff00000 /* y is NaN */ + || (ky == 0x7ff00000 && t.i[LOW_HALF] != 0))) + return (x * y) / (x * y); + else + return x; + } + } + } +} +strong_alias (__ieee754_remainder, __remainder_finite) |