]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/e_remainder.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 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 /* ************************************************************************/
36 #include <math_private.h>
38 /**************************************************************************/
39 /* An ultimate remainder routine. Given two IEEE double machine numbers x */
40 /* ,y it computes the correctly rounded (to nearest) value of remainder */
41 /**************************************************************************/
43 __ieee754_remainder (double x
, double y
)
46 int4 kx
, ky
, n
, nn
, n1
, m1
, l
;
47 mynumber u
, t
, w
= { { 0, 0 } }, v
= { { 0, 0 } }, ww
= { { 0, 0 } }, r
;
50 kx
= u
.i
[HIGH_HALF
] & 0x7fffffff; /* no sign for x*/
51 t
.i
[HIGH_HALF
] &= 0x7fffffff; /*no sign for y */
53 /*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/
54 if (kx
< 0x7fe00000 && ky
< 0x7ff00000 && ky
>= 0x03500000)
56 SET_RESTORE_ROUND_NOEX (FE_TONEAREST
);
57 if (kx
+ 0x00100000 < ky
)
59 if ((kx
- 0x01500000) < ky
)
62 v
.i
[HIGH_HALF
] = t
.i
[HIGH_HALF
];
63 d
= (z
+ big
.x
) - big
.x
;
64 xx
= (x
- d
* v
.x
) - d
* (t
.x
- v
.x
);
65 if (d
- z
!= 0.5 && d
- z
!= -0.5)
66 return (xx
!= 0) ? xx
: ((x
> 0) ? ZERO
.x
: nZERO
.x
);
69 if (ABS (xx
) > 0.5 * t
.x
)
70 return (z
> d
) ? xx
- t
.x
: xx
+ t
.x
;
74 } /* (kx<(ky+0x01500000)) */
79 nn
= (n
& 0x7ff00000) + 0x01400000;
82 l
= (kx
- nn
) & 0xfff00000;
87 r
.i
[HIGH_HALF
] = m1
- l
;
89 w
.i
[HIGH_HALF
] = n
+ l
;
90 ww
.i
[HIGH_HALF
] = (n1
) ? n1
+ l
: n1
;
91 d
= (z
+ big
.x
) - big
.x
;
92 u
.x
= (u
.x
- d
* w
.x
) - d
* ww
.x
;
93 l
= (u
.i
[HIGH_HALF
] & 0x7ff00000) - nn
;
99 d
= (z
+ big
.x
) - big
.x
;
100 u
.x
= (u
.x
- d
* w
.x
) - d
* ww
.x
;
101 if (ABS (u
.x
) < 0.5 * t
.x
)
102 return (u
.x
!= 0) ? u
.x
: ((x
> 0) ? ZERO
.x
: nZERO
.x
);
104 if (ABS (u
.x
) > 0.5 * t
.x
)
105 return (d
> z
) ? u
.x
+ t
.x
: u
.x
- t
.x
;
108 z
= u
.x
/ t
.x
; d
= (z
+ big
.x
) - big
.x
;
109 return ((u
.x
- d
* w
.x
) - d
* ww
.x
);
112 } /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */
115 if (kx
< 0x7fe00000 && ky
< 0x7ff00000 && (ky
> 0 || t
.i
[LOW_HALF
] != 0))
117 y
= ABS (y
) * t128
.x
;
118 z
= __ieee754_remainder (x
, y
) * t128
.x
;
119 z
= __ieee754_remainder (z
, y
) * tm128
.x
;
124 if ((kx
& 0x7ff00000) == 0x7fe00000 && ky
< 0x7ff00000 &&
125 (ky
> 0 || t
.i
[LOW_HALF
] != 0))
128 z
= 2.0 * __ieee754_remainder (0.5 * x
, y
);
130 if (d
<= ABS (d
- y
))
133 return (z
> 0) ? z
- y
: z
+ y
;
135 else /* if x is too big */
137 if (ky
== 0 && t
.i
[LOW_HALF
] == 0) /* y = 0 */
138 return (x
* y
) / (x
* y
);
139 else if (kx
>= 0x7ff00000 /* x not finite */
140 || (ky
> 0x7ff00000 /* y is NaN */
141 || (ky
== 0x7ff00000 && t
.i
[LOW_HALF
] != 0)))
142 return (x
* y
) / (x
* y
);
149 strong_alias (__ieee754_remainder
, __remainder_finite
)