1
2 /* @(#)e_fmod.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 /*
15 * fmod(x,y)
16 * Return x mod y in exact arithmetic
17 * Method: shift and subtract
18 */
19
20 #include "fdlibm.h"
21
22 #ifdef _NEED_FLOAT64
23
24 static const __float64
25 one = _F_64(1.0),
26 Zero[] = { _F_64(0.0), _F_64(-0.0) };
27
28 __float64
fmod64(__float64 x,__float64 y)29 fmod64(__float64 x, __float64 y)
30 {
31 __int32_t n, hx, hy, hz, ix, iy, sx, i;
32 __uint32_t lx, ly, lz;
33
34 EXTRACT_WORDS(hx, lx, x);
35 EXTRACT_WORDS(hy, ly, y);
36 sx = hx & 0x80000000; /* sign of x */
37 hx ^= sx; /* |x| */
38 hy &= 0x7fffffff; /* |y| */
39
40 /* purge off exception values */
41 if (isnan(x) || isnan(y)) /* x or y nan, return nan */
42 return x + y;
43
44 if (isinf(x)) /* x == inf, domain error */
45 return __math_invalid(x);
46
47 if ((hy | ly) == 0) /* y=0, domain error */
48 return __math_invalid(y);
49
50 if (hx <= hy) {
51 if ((hx < hy) || (lx < ly))
52 return x; /* |x|<|y| return x */
53 if (lx == ly)
54 return Zero[(__uint32_t)sx >> 31]; /* |x|=|y| return x*0*/
55 }
56
57 /* determine ix = ilogb(x) */
58 if (hx < 0x00100000) { /* subnormal x */
59 if (hx == 0) {
60 for (ix = -1043, i = lx; i > 0; i <<= 1)
61 ix -= 1;
62 } else {
63 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
64 ix -= 1;
65 }
66 } else
67 ix = (hx >> 20) - 1023;
68
69 /* determine iy = ilogb(y) */
70 if (hy < 0x00100000) { /* subnormal y */
71 if (hy == 0) {
72 for (iy = -1043, i = ly; i > 0; i <<= 1)
73 iy -= 1;
74 } else {
75 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
76 iy -= 1;
77 }
78 } else
79 iy = (hy >> 20) - 1023;
80
81 /* set up {hx,lx}, {hy,ly} and align y to x */
82 if (ix >= -1022)
83 hx = 0x00100000 | (0x000fffff & hx);
84 else { /* subnormal x, shift x to normal */
85 n = -1022 - ix;
86 if (n <= 31) {
87 hx = (hx << n) | (lx >> (32 - n));
88 lx <<= n;
89 } else {
90 hx = lx << (n - 32);
91 lx = 0;
92 }
93 }
94 if (iy >= -1022)
95 hy = 0x00100000 | (0x000fffff & hy);
96 else { /* subnormal y, shift y to normal */
97 n = -1022 - iy;
98 if (n <= 31) {
99 hy = (hy << n) | (ly >> (32 - n));
100 ly <<= n;
101 } else {
102 hy = ly << (n - 32);
103 ly = 0;
104 }
105 }
106
107 /* fix point fmod */
108 n = ix - iy;
109 while (n--) {
110 hz = hx - hy;
111 lz = lx - ly;
112 if (lx < ly)
113 hz -= 1;
114 if (hz < 0) {
115 hx = hx + hx + (lx >> 31);
116 lx = lx + lx;
117 } else {
118 if ((hz | lz) == 0) /* return sign(x)*0 */
119 return Zero[(__uint32_t)sx >> 31];
120 hx = hz + hz + (lz >> 31);
121 lx = lz + lz;
122 }
123 }
124 hz = hx - hy;
125 lz = lx - ly;
126 if (lx < ly)
127 hz -= 1;
128 if (hz >= 0) {
129 hx = hz;
130 lx = lz;
131 }
132
133 /* convert back to floating value and restore the sign */
134 if ((hx | lx) == 0) /* return sign(x)*0 */
135 return Zero[(__uint32_t)sx >> 31];
136 while (hx < 0x00100000) { /* normalize x */
137 hx = hx + hx + (lx >> 31);
138 lx = lx + lx;
139 iy -= 1;
140 }
141 if (iy >= -1022) { /* normalize output */
142 hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
143 INSERT_WORDS(x, hx | sx, lx);
144 } else { /* subnormal output */
145 n = -1022 - iy;
146 if (n <= 20) {
147 lx = (lx >> n) | ((__uint32_t)hx << (32 - n));
148 hx >>= n;
149 } else if (n <= 31) {
150 lx = (hx << (32 - n)) | (lx >> n);
151 hx = sx;
152 } else {
153 lx = hx >> (n - 32);
154 hx = sx;
155 }
156 INSERT_WORDS(x, hx | sx, lx);
157 x *= one; /* create necessary signal */
158 }
159 return x; /* exact output */
160 }
161
162 _MATH_ALIAS_d_dd(fmod)
163
164 #endif /* _NEED_FLOAT64 */
165
