1 /* ef_fmod.c -- float version of e_fmod.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16 /*
17 * fmodf(x,y)
18 * Return x mod y in exact arithmetic
19 * Method: shift and subtract
20 */
21
22 #include "fdlibm.h"
23
24 static const float one = 1.0, Zero[] = {
25 0.0,
26 -0.0,
27 };
28
29 float
fmodf(float x,float y)30 fmodf(float x, float y)
31 {
32 __int32_t n, hx, hy, hz, ix, iy, sx, i;
33
34 GET_FLOAT_WORD(hx, x);
35 GET_FLOAT_WORD(hy, 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_invalidf(x);
46
47 if (hy == 0) /* y=0, domain error */
48 return __math_invalidf(y);
49
50 if (hx < hy)
51 return x; /* |x|<|y| return x */
52 if (hx == hy)
53 return Zero[(__uint32_t)sx >> 31]; /* |x|=|y| return x*0*/
54
55 /* Note: y cannot be zero if we reach here. */
56
57 /* determine ix = ilogb(x) */
58 if (FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */
59 for (ix = -126, i = (hx << 8); i > 0; i <<= 1)
60 ix -= 1;
61 } else
62 ix = (hx >> 23) - 127;
63
64 /* determine iy = ilogb(y) */
65 if (FLT_UWORD_IS_SUBNORMAL(hy)) { /* subnormal y */
66 for (iy = -126, i = (hy << 8); i >= 0; i <<= 1)
67 iy -= 1;
68 } else
69 iy = (hy >> 23) - 127;
70
71 /* set up {hx,lx}, {hy,ly} and align y to x */
72 if (ix >= -126)
73 hx = 0x00800000 | (0x007fffff & hx);
74 else { /* subnormal x, shift x to normal */
75 n = -126 - ix;
76 hx = hx << n;
77 }
78 if (iy >= -126)
79 hy = 0x00800000 | (0x007fffff & hy);
80 else { /* subnormal y, shift y to normal */
81 n = -126 - iy;
82 hy = hy << n;
83 }
84
85 /* fix point fmod */
86 n = ix - iy;
87 while (n--) {
88 hz = hx - hy;
89 if (hz < 0) {
90 hx = hx + hx;
91 } else {
92 if (hz == 0) /* return sign(x)*0 */
93 return Zero[(__uint32_t)sx >> 31];
94 hx = hz + hz;
95 }
96 }
97 hz = hx - hy;
98 if (hz >= 0) {
99 hx = hz;
100 }
101
102 /* convert back to floating value and restore the sign */
103 if (hx == 0) /* return sign(x)*0 */
104 return Zero[(__uint32_t)sx >> 31];
105 while (hx < 0x00800000) { /* normalize x */
106 hx = hx + hx;
107 iy -= 1;
108 }
109 if (iy >= -126) { /* normalize output */
110 hx = ((hx - 0x00800000) | ((iy + 127) << 23));
111 SET_FLOAT_WORD(x, hx | sx);
112 } else { /* subnormal output */
113 /* If denormals are not supported, this code will generate a
114 zero representation. */
115 n = -126 - iy;
116 hx >>= n;
117 SET_FLOAT_WORD(x, hx | sx);
118 x *= one; /* create necessary signal */
119 }
120 return x; /* exact output */
121 }
122
123 _MATH_ALIAS_f_ff(fmod)
124
