]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/e_remainder.c
33cb48b8e3083eef0c06a82829807416c38bdcf9
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2019 Free Software Foundation, Inc.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /**************************************************************************/
20 /* MODULE_NAME urem.c */
22 /* FUNCTION: uremainder */
24 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
25 /* ,y it computes the correctly rounded (to nearest) value of remainder */
26 /* of dividing x by y. */
27 /* Assumption: Machine arithmetic operations are performed in */
28 /* round to nearest mode of IEEE 754 standard. */
30 /* ************************************************************************/
37 #include <math_private.h>
38 #include <fenv_private.h>
40 /**************************************************************************/
41 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
42 /* ,y it computes the correctly rounded (to nearest) value of remainder */
43 /**************************************************************************/
45 __ieee754_remainder (double x
, double y
)
48 int4 kx
, ky
, n
, nn
, n1
, m1
, l
;
49 mynumber u
, t
, w
= { { 0, 0 } }, v
= { { 0, 0 } }, ww
= { { 0, 0 } }, r
;
52 kx
= u
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign for x*/
53 t
.i
[HIGH_HALF
] &= 0x7fffffff; /*no sign for y */
55 /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/
56 if (kx
< 0x7fe00000 && ky
< 0x7ff00000 && ky
>= 0x03500000)
58 SET_RESTORE_ROUND_NOEX (FE_TONEAREST
);
59 if (kx
+ 0x00100000 < ky
)
61 if ((kx
- 0x01500000) < ky
)
64 v
.i
[HIGH_HALF
] = t
.i
[HIGH_HALF
];
65 d
= (z
+ big
.x
) - big
.x
;
66 xx
= (x
- d
* v
.x
) - d
* (t
.x
- v
.x
);
67 if (d
- z
!= 0.5 && d
- z
!= -0.5)
68 return (xx
!= 0) ? xx
: ((x
> 0) ? ZERO
.x
: nZERO
.x
);
71 if (fabs (xx
) > 0.5 * t
.x
)
72 return (z
> d
) ? xx
- t
.x
: xx
+ t
.x
;
76 } /* (kx<(ky+0x01500000)) */
81 nn
= (n
& 0x7ff00000) + 0x01400000;
84 l
= (kx
- nn
) & 0xfff00000;
89 r
.i
[HIGH_HALF
] = m1
- l
;
91 w
.i
[HIGH_HALF
] = n
+ l
;
92 ww
.i
[HIGH_HALF
] = (n1
) ? n1
+ l
: n1
;
93 d
= (z
+ big
.x
) - big
.x
;
94 u
.x
= (u
.x
- d
* w
.x
) - d
* ww
.x
;
95 l
= (u
.i
[HIGH_HALF
] & 0x7ff00000) - nn
;
101 d
= (z
+ big
.x
) - big
.x
;
102 u
.x
= (u
.x
- d
* w
.x
) - d
* ww
.x
;
103 if (fabs (u
.x
) < 0.5 * t
.x
)
104 return (u
.x
!= 0) ? u
.x
: ((x
> 0) ? ZERO
.x
: nZERO
.x
);
106 if (fabs (u
.x
) > 0.5 * t
.x
)
107 return (d
> z
) ? u
.x
+ t
.x
: u
.x
- t
.x
;
110 z
= u
.x
/ t
.x
; d
= (z
+ big
.x
) - big
.x
;
111 return ((u
.x
- d
* w
.x
) - d
* ww
.x
);
114 } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */
117 if (kx
< 0x7fe00000 && ky
< 0x7ff00000 && (ky
> 0 || t
.i
[LOW_HALF
] != 0))
119 y
= fabs (y
) * t128
.x
;
120 z
= __ieee754_remainder (x
, y
) * t128
.x
;
121 z
= __ieee754_remainder (z
, y
) * tm128
.x
;
126 if ((kx
& 0x7ff00000) == 0x7fe00000 && ky
< 0x7ff00000 &&
127 (ky
> 0 || t
.i
[LOW_HALF
] != 0))
130 z
= 2.0 * __ieee754_remainder (0.5 * x
, y
);
132 if (d
<= fabs (d
- y
))
137 return (z
> 0) ? z
- y
: z
+ y
;
139 else /* if x is too big */
141 if (ky
== 0 && t
.i
[LOW_HALF
] == 0) /* y = 0 */
142 return (x
* y
) / (x
* y
);
143 else if (kx
>= 0x7ff00000 /* x not finite */
144 || (ky
> 0x7ff00000 /* y is NaN */
145 || (ky
== 0x7ff00000 && t
.i
[LOW_HALF
] != 0)))
146 return (x
* y
) / (x
* y
);
153 strong_alias (__ieee754_remainder
, __remainder_finite
)