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