1
2 /* @(#)e_remainder.c 5.1 93/09/24 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14 /* remainder(x,p)
15 * Return :
16 * returns x REM p = x - [x/p]*p as if in infinite
17 * precise arithmetic, where [x/p] is the (infinite bit)
18 * integer nearest x/p (in half way case choose the even one).
19 * Method :
20 * Based on fmod() return x-[x/p]chopped*p exactlp.
21 */
22
23 #include "fdlibm.h"
24
25 #ifdef _NEED_FLOAT64
26
27 static const __float64 zero = _F_64(0.0);
28
29 __float64
remainder64(__float64 x,__float64 p)30 remainder64(__float64 x, __float64 p)
31 {
32 __int32_t hx, hp;
33 __uint32_t sx, lx, lp;
34 __float64 p_half;
35
36 EXTRACT_WORDS(hx, lx, x);
37 EXTRACT_WORDS(hp, lp, p);
38 sx = hx & 0x80000000;
39 hp &= 0x7fffffff;
40 hx &= 0x7fffffff;
41
42 /* purge off exception values */
43 if (isnan(x) || isnan(p))
44 return x + p;
45
46 if (isinf(x) || (hp | lp) == 0)
47 return __math_invalid(x);
48
49 if ((hp | lp) == 0)
50 return (x * p) / (x * p); /* p = 0 */
51
52 if (hp <= 0x7fdfffff)
53 x = fmod(x, p + p); /* now x < 2p */
54
55 if (((hx - hp) | (lx - lp)) == 0)
56 return zero * x;
57 x = fabs64(x);
58 p = fabs64(p);
59 if (hp < 0x00200000) {
60 if (x + x > p) {
61 x -= p;
62 if (x + x >= p)
63 x -= p;
64 }
65 } else {
66 p_half = _F_64(0.5) * p;
67 if (x > p_half) {
68 x -= p;
69 if (x >= p_half)
70 x -= p;
71 }
72 }
73 GET_HIGH_WORD(hx, x);
74 SET_HIGH_WORD(x, hx ^ sx);
75 return x;
76 }
77
78 _MATH_ALIAS_d_dd(remainder)
79
80 #endif /* _NEED_FLOAT64 */
81
