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 }