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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak __fmodf = fmodf
30
31 #include "libm.h"
32
33 /* INDENT OFF */
34 static const int
35 is = (int)0x80000000,
36 im = 0x007fffff,
37 ii = 0x7f800000,
38 iu = 0x00800000;
39 /* INDENT ON */
40
41 static const float zero = 0.0;
42
43 float
44 fmodf(float x, float y) {
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 /* INDENT OFF */
63 /*
64 * scale x,y to "normal" with
65 * ny = exponent of y
66 * nd = exponent of x minus exponent of y
67 */
68 /* INDENT ON */
69 ny = iy >> 23;
70 k = ix >> 23;
71
72 /* special case for subnormal y or x */
73 if (ny == 0) {
74 ny = 1;
75 while (iy < iu) {
76 ny -= 1;
77 iy += iy;
78 }
79 nd = k - ny;
80 if (k == 0) {
81 nd += 1;
82 while (ix < iu) {
83 nd -= 1;
84 ix += ix;
85 }
86 } else {
87 ix = iu | (ix & im);
88 }
89 } else {
90 nd = k - ny;
91 ix = iu | (ix & im);
92 iy = iu | (iy & im);
93 }
94
95 /* fix point fmod for normalized ix and iy */
96 /* INDENT OFF */
97 /*
98 * while (nd--) {
99 * iz = ix - iy;
100 * if (iz < 0)
101 * ix = ix + ix;
102 * else if (iz == 0) {
103 * *(int *) &w = is & hx;
104 * return w;
105 * }
106 * else
107 * ix = iz + iz;
108 * }
109 */
110 /* INDENT ON */
111 /* unroll the above loop 4 times to gain performance */
112 k = nd >> 2;
113 nd -= k << 2;
114 while (k--) {
115 iz = ix - iy;
116 if (iz >= 0)
117 ix = iz + iz;
118 else
119 ix += ix;
120 iz = ix - iy;
121 if (iz >= 0)
122 ix = iz + iz;
123 else
124 ix += ix;
125 iz = ix - iy;
126 if (iz >= 0)
127 ix = iz + iz;
128 else
129 ix += ix;
130 iz = ix - iy;
131 if (iz >= 0)
132 ix = iz + iz;
133 else
134 ix += ix;
135 if (iz == 0) {
136 *(int *)&w = is & hx;
137 return (w);
138 }
139 }
140 while (nd--) {
141 iz = ix - iy;
142 if (iz >= 0)
143 ix = iz + iz;
144 else
145 ix += ix;
146 }
147 /* end of unrolling */
148
149 iz = ix - iy;
150 if (iz >= 0)
151 ix = iz;
152
153 /* convert back to floating value and restore the sign */
154 if (ix == 0) {
155 *(int *)&w = is & hx;
156 return (w);
157 }
158 while (ix < iu) {
159 ix += ix;
160 ny -= 1;
161 }
162 while (ix > (iu + iu)) {
163 ny += 1;
164 ix >>= 1;
165 }
166 if (ny > 0) {
167 *(int *)&w = (is & hx) | (ix & im) | (ny << 23);
168 } else {
169 /* subnormal output */
170 k = -ny + 1;
171 ix >>= k;
172 *(int *)&w = (is & hx) | ix;
173 }
174 }
175 return (w);
176 }
|
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 }
|