1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __fmodf = fmodf
32
33 #include "libm.h"
34
35 static const int is = (int)0x80000000,
36 im = 0x007fffff,
37 ii = 0x7f800000,
38 iu = 0x00800000;
39
40 static const float zero = 0.0;
41
42 float
43 fmodf(float x, float y)
44 {
45 float w;
46 int hx, ix, iy, iz, k, ny, nd;
47
48 hx = *(int *)&x;
49 ix = hx & 0x7fffffff;
50 iy = *(int *)&y & 0x7fffffff;
51
52 /* purge off exception values */
53 if (ix >= ii || iy > ii || iy == 0) {
54 w = x * y;
55 w = w / w;
56 } else if (ix <= iy) {
57 if (ix < iy)
58 w = x; /* return x if |x|<|y| */
59 else
60 w = zero * x; /* return sign(x)*0.0 */
61 } else {
62
63 /*
64 * scale x,y to "normal" with
65 * ny = exponent of y
66 * nd = exponent of x minus exponent of y
67 */
68 ny = iy >> 23;
69 k = ix >> 23;
70
71 /* special case for subnormal y or x */
72 if (ny == 0) {
73 ny = 1;
74
75 while (iy < iu) {
76 ny -= 1;
77 iy += iy;
78 }
79
80 nd = k - ny;
81
82 if (k == 0) {
83 nd += 1;
84
85 while (ix < iu) {
86 nd -= 1;
87 ix += ix;
88 }
89 } else {
90 ix = iu | (ix & im);
91 }
92 } else {
93 nd = k - ny;
94 ix = iu | (ix & im);
95 iy = iu | (iy & im);
96 }
97
98 /*
99 * fix point fmod for normalized ix and iy
100 */
101
102 /* BEGIN CSTYLED */
103 /*
104 * while (nd--) {
105 * iz = ix - iy;
106 * if (iz < 0)
107 * ix = ix + ix;
108 * else if (iz == 0) {
109 * *(int *) &w = is & hx;
110 * return w;
111 * }
112 * else
113 * ix = iz + iz;
114 * }
115 */
116 /* END CSTYLED */
117
118 /*
119 * unroll the above loop 4 times to gain performance
120 */
121 k = nd >> 2;
122 nd -= k << 2;
123
124 while (k--) {
125 iz = ix - iy;
126
127 if (iz >= 0)
128 ix = iz + iz;
129 else
130 ix += ix;
131
132 iz = ix - iy;
133
134 if (iz >= 0)
135 ix = iz + iz;
136 else
137 ix += ix;
138
139 iz = ix - iy;
140
141 if (iz >= 0)
142 ix = iz + iz;
143 else
144 ix += ix;
145
146 iz = ix - iy;
147
148 if (iz >= 0)
149 ix = iz + iz;
150 else
151 ix += ix;
152
153 if (iz == 0) {
154 *(int *)&w = is & hx;
155 return (w);
156 }
157 }
158
159 while (nd--) {
160 iz = ix - iy;
161
162 if (iz >= 0)
163 ix = iz + iz;
164 else
165 ix += ix;
166 }
167
168 /* end of unrolling */
169
170 iz = ix - iy;
171
172 if (iz >= 0)
173 ix = iz;
174
175 /* convert back to floating value and restore the sign */
176 if (ix == 0) {
177 *(int *)&w = is & hx;
178 return (w);
179 }
180
181 while (ix < iu) {
182 ix += ix;
183 ny -= 1;
184 }
185
186 while (ix > (iu + iu)) {
187 ny += 1;
188 ix >>= 1;
189 }
190
191 if (ny > 0) {
192 *(int *)&w = (is & hx) | (ix & im) | (ny << 23);
193 } else {
194 /* subnormal output */
195 k = -ny + 1;
196 ix >>= k;
197 *(int *)&w = (is & hx) | ix;
198 }
199 }
200
201 return (w);
202 }