]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ieee754/dbl-64/e_fmod.c
Use glibc_likely instead __builtin_expect.
[thirdparty/glibc.git] / sysdeps / ieee754 / dbl-64 / e_fmod.c
CommitLineData
f7eac6eb
RM
1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
0ac5ae23 7 * software is freely granted, provided that this notice
f7eac6eb
RM
8 * is preserved.
9 * ====================================================
10 */
11
0ac5ae23 12/*
f7eac6eb
RM
13 * __ieee754_fmod(x,y)
14 * Return x mod y in exact arithmetic
15 * Method: shift and subtract
16 */
17
1ed0291c
RH
18#include <math.h>
19#include <math_private.h>
f7eac6eb 20
c5d5d574 21static const double one = 1.0, Zero[] = { 0.0, -0.0, };
f7eac6eb 22
0ac5ae23
UD
23double
24__ieee754_fmod (double x, double y)
f7eac6eb 25{
c5d5d574
OB
26 int32_t n, hx, hy, hz, ix, iy, sx, i;
27 u_int32_t lx, ly, lz;
f7eac6eb 28
c5d5d574
OB
29 EXTRACT_WORDS (hx, lx, x);
30 EXTRACT_WORDS (hy, ly, y);
31 sx = hx & 0x80000000; /* sign of x */
32 hx ^= sx; /* |x| */
33 hy &= 0x7fffffff; /* |y| */
f7eac6eb 34
c5d5d574
OB
35 /* purge off exception values */
36 if ((hy | ly) == 0 || (hx >= 0x7ff00000) || /* y=0,or x not finite */
37 ((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */
38 return (x * y) / (x * y);
39 if (hx <= hy)
40 {
41 if ((hx < hy) || (lx < ly))
42 return x; /* |x|<|y| return x */
43 if (lx == ly)
44 return Zero[(u_int32_t) sx >> 31]; /* |x|=|y| return x*0*/
45 }
f7eac6eb 46
c5d5d574 47 /* determine ix = ilogb(x) */
a1ffb40e 48 if (__glibc_unlikely (hx < 0x00100000)) /* subnormal x */
c5d5d574
OB
49 {
50 if (hx == 0)
51 {
52 for (ix = -1043, i = lx; i > 0; i <<= 1)
53 ix -= 1;
54 }
55 else
56 {
57 for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
58 ix -= 1;
59 }
60 }
61 else
62 ix = (hx >> 20) - 1023;
f7eac6eb 63
c5d5d574 64 /* determine iy = ilogb(y) */
a1ffb40e 65 if (__glibc_unlikely (hy < 0x00100000)) /* subnormal y */
c5d5d574
OB
66 {
67 if (hy == 0)
68 {
69 for (iy = -1043, i = ly; i > 0; i <<= 1)
70 iy -= 1;
71 }
72 else
73 {
74 for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
75 iy -= 1;
76 }
77 }
78 else
79 iy = (hy >> 20) - 1023;
f7eac6eb 80
c5d5d574 81 /* set up {hx,lx}, {hy,ly} and align y to x */
a1ffb40e 82 if (__glibc_likely (ix >= -1022))
c5d5d574
OB
83 hx = 0x00100000 | (0x000fffff & hx);
84 else /* subnormal x, shift x to normal */
85 {
86 n = -1022 - ix;
87 if (n <= 31)
88 {
89 hx = (hx << n) | (lx >> (32 - n));
90 lx <<= n;
91 }
92 else
93 {
94 hx = lx << (n - 32);
95 lx = 0;
96 }
97 }
a1ffb40e 98 if (__glibc_likely (iy >= -1022))
c5d5d574
OB
99 hy = 0x00100000 | (0x000fffff & hy);
100 else /* subnormal y, shift y to normal */
101 {
102 n = -1022 - iy;
103 if (n <= 31)
104 {
105 hy = (hy << n) | (ly >> (32 - n));
106 ly <<= n;
f7eac6eb 107 }
c5d5d574
OB
108 else
109 {
110 hy = ly << (n - 32);
111 ly = 0;
f7eac6eb 112 }
c5d5d574 113 }
f7eac6eb 114
c5d5d574
OB
115 /* fix point fmod */
116 n = ix - iy;
117 while (n--)
118 {
119 hz = hx - hy; lz = lx - ly; if (lx < ly)
120 hz -= 1;
121 if (hz < 0)
122 {
123 hx = hx + hx + (lx >> 31); lx = lx + lx;
f7eac6eb 124 }
c5d5d574
OB
125 else
126 {
127 if ((hz | lz) == 0) /* return sign(x)*0 */
128 return Zero[(u_int32_t) sx >> 31];
129 hx = hz + hz + (lz >> 31); lx = lz + lz;
130 }
131 }
132 hz = hx - hy; lz = lx - ly; if (lx < ly)
133 hz -= 1;
134 if (hz >= 0)
135 {
136 hx = hz; lx = lz;
137 }
f7eac6eb 138
c5d5d574
OB
139 /* convert back to floating value and restore the sign */
140 if ((hx | lx) == 0) /* return sign(x)*0 */
141 return Zero[(u_int32_t) sx >> 31];
142 while (hx < 0x00100000) /* normalize x */
143 {
144 hx = hx + hx + (lx >> 31); lx = lx + lx;
145 iy -= 1;
146 }
a1ffb40e 147 if (__glibc_likely (iy >= -1022)) /* normalize output */
c5d5d574
OB
148 {
149 hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
150 INSERT_WORDS (x, hx | sx, lx);
151 }
152 else /* subnormal output */
153 {
154 n = -1022 - iy;
155 if (n <= 20)
156 {
157 lx = (lx >> n) | ((u_int32_t) hx << (32 - n));
158 hx >>= n;
159 }
160 else if (n <= 31)
161 {
162 lx = (hx << (32 - n)) | (lx >> n); hx = sx;
f7eac6eb 163 }
c5d5d574
OB
164 else
165 {
166 lx = hx >> (n - 32); hx = sx;
f7eac6eb 167 }
c5d5d574
OB
168 INSERT_WORDS (x, hx | sx, lx);
169 x *= one; /* create necessary signal */
170 }
171 return x; /* exact output */
f7eac6eb 172}
0ac5ae23 173strong_alias (__ieee754_fmod, __fmod_finite)