1 /* $NetBSD: sfrem.c,v 1.5 2012/02/04 17:03:10 skrll Exp $ */
2
3 /* $OpenBSD: sfrem.c,v 1.4 2001/03/29 03:58:19 mickey Exp $ */
4
5 /*
6 * Copyright 1996 1995 by Open Software Foundation, Inc.
7 * All Rights Reserved
8 *
9 * Permission to use, copy, modify, and distribute this software and
10 * its documentation for any purpose and without fee is hereby granted,
11 * provided that the above copyright notice appears in all copies and
12 * that both the copyright notice and this permission notice appear in
13 * supporting documentation.
14 *
15 * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE
16 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
17 * FOR A PARTICULAR PURPOSE.
18 *
19 * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR
20 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
21 * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT,
22 * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
23 * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
24 *
25 */
26 /*
27 * pmk1.1
28 */
29 /*
30 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
31 *
32 * To anyone who acknowledges that this file is provided "AS IS"
33 * without any express or implied warranty:
34 * permission to use, copy, modify, and distribute this file
35 * for any purpose is hereby granted without fee, provided that
36 * the above copyright notice and this notice appears in all
37 * copies, and that the name of Hewlett-Packard Company not be
38 * used in advertising or publicity pertaining to distribution
39 * of the software without specific, written prior permission.
40 * Hewlett-Packard Company makes no representations about the
41 * suitability of this software for any purpose.
42 */
43
44 #include <sys/cdefs.h>
45 __KERNEL_RCSID(0, "$NetBSD: sfrem.c,v 1.5 2012/02/04 17:03:10 skrll Exp $");
46
47 #include "../spmath/float.h"
48 #include "../spmath/sgl_float.h"
49
50 /*
51 * Single Precision Floating-point Remainder
52 */
53 int
sgl_frem(sgl_floating_point * srcptr1,sgl_floating_point * srcptr2,sgl_floating_point * dstptr,unsigned int * status)54 sgl_frem(sgl_floating_point *srcptr1, sgl_floating_point *srcptr2,
55 sgl_floating_point *dstptr, unsigned int *status)
56 {
57 register unsigned int opnd1, opnd2, result;
58 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
59 register int roundup = false;
60
61 opnd1 = *srcptr1;
62 opnd2 = *srcptr2;
63 /*
64 * check first operand for NaN's or infinity
65 */
66 if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
67 if (Sgl_iszero_mantissa(opnd1)) {
68 if (Sgl_isnotnan(opnd2)) {
69 /* invalid since first operand is infinity */
70 if (Is_invalidtrap_enabled())
71 return(INVALIDEXCEPTION);
72 Set_invalidflag();
73 Sgl_makequietnan(result);
74 *dstptr = result;
75 return(NOEXCEPTION);
76 }
77 }
78 else {
79 /*
80 * is NaN; signaling or quiet?
81 */
82 if (Sgl_isone_signaling(opnd1)) {
83 /* trap if INVALIDTRAP enabled */
84 if (Is_invalidtrap_enabled())
85 return(INVALIDEXCEPTION);
86 /* make NaN quiet */
87 Set_invalidflag();
88 Sgl_set_quiet(opnd1);
89 }
90 /*
91 * is second operand a signaling NaN?
92 */
93 else if (Sgl_is_signalingnan(opnd2)) {
94 /* trap if INVALIDTRAP enabled */
95 if (Is_invalidtrap_enabled())
96 return(INVALIDEXCEPTION);
97 /* make NaN quiet */
98 Set_invalidflag();
99 Sgl_set_quiet(opnd2);
100 *dstptr = opnd2;
101 return(NOEXCEPTION);
102 }
103 /*
104 * return quiet NaN
105 */
106 *dstptr = opnd1;
107 return(NOEXCEPTION);
108 }
109 }
110 /*
111 * check second operand for NaN's or infinity
112 */
113 if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
114 if (Sgl_iszero_mantissa(opnd2)) {
115 /*
116 * return first operand
117 */
118 *dstptr = opnd1;
119 return(NOEXCEPTION);
120 }
121 /*
122 * is NaN; signaling or quiet?
123 */
124 if (Sgl_isone_signaling(opnd2)) {
125 /* trap if INVALIDTRAP enabled */
126 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
127 /* make NaN quiet */
128 Set_invalidflag();
129 Sgl_set_quiet(opnd2);
130 }
131 /*
132 * return quiet NaN
133 */
134 *dstptr = opnd2;
135 return(NOEXCEPTION);
136 }
137 /*
138 * check second operand for zero
139 */
140 if (Sgl_iszero_exponentmantissa(opnd2)) {
141 /* invalid since second operand is zero */
142 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
143 Set_invalidflag();
144 Sgl_makequietnan(result);
145 *dstptr = result;
146 return(NOEXCEPTION);
147 }
148
149 /*
150 * get sign of result
151 */
152 result = opnd1;
153
154 /*
155 * check for denormalized operands
156 */
157 if (opnd1_exponent == 0) {
158 /* check for zero */
159 if (Sgl_iszero_mantissa(opnd1)) {
160 *dstptr = opnd1;
161 return(NOEXCEPTION);
162 }
163 /* normalize, then continue */
164 opnd1_exponent = 1;
165 Sgl_normalize(opnd1,opnd1_exponent);
166 }
167 else {
168 Sgl_clear_signexponent_set_hidden(opnd1);
169 }
170 if (opnd2_exponent == 0) {
171 /* normalize, then continue */
172 opnd2_exponent = 1;
173 Sgl_normalize(opnd2,opnd2_exponent);
174 }
175 else {
176 Sgl_clear_signexponent_set_hidden(opnd2);
177 }
178
179 /* find result exponent and divide step loop count */
180 dest_exponent = opnd2_exponent - 1;
181 stepcount = opnd1_exponent - opnd2_exponent;
182
183 /*
184 * check for opnd1/opnd2 < 1
185 */
186 if (stepcount < 0) {
187 /*
188 * check for opnd1/opnd2 > 1/2
189 *
190 * In this case n will round to 1, so
191 * r = opnd1 - opnd2
192 */
193 if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
194 Sgl_all(result) = ~Sgl_all(result); /* set sign */
195 /* align opnd2 with opnd1 */
196 Sgl_leftshiftby1(opnd2);
197 Sgl_subtract(opnd2,opnd1,opnd2);
198 /* now normalize */
199 while (Sgl_iszero_hidden(opnd2)) {
200 Sgl_leftshiftby1(opnd2);
201 dest_exponent--;
202 }
203 Sgl_set_exponentmantissa(result,opnd2);
204 goto testforunderflow;
205 }
206 /*
207 * opnd1/opnd2 <= 1/2
208 *
209 * In this case n will round to zero, so
210 * r = opnd1
211 */
212 Sgl_set_exponentmantissa(result,opnd1);
213 dest_exponent = opnd1_exponent;
214 goto testforunderflow;
215 }
216
217 /*
218 * Generate result
219 *
220 * Do iterative subtract until remainder is less than operand 2.
221 */
222 while (stepcount-- > 0 && Sgl_all(opnd1)) {
223 if (Sgl_isnotlessthan(opnd1,opnd2))
224 Sgl_subtract(opnd1,opnd2,opnd1);
225 Sgl_leftshiftby1(opnd1);
226 }
227 /*
228 * Do last subtract, then determine which way to round if remainder
229 * is exactly 1/2 of opnd2
230 */
231 if (Sgl_isnotlessthan(opnd1,opnd2)) {
232 Sgl_subtract(opnd1,opnd2,opnd1);
233 roundup = true;
234 }
235 if (stepcount > 0 || Sgl_iszero(opnd1)) {
236 /* division is exact, remainder is zero */
237 Sgl_setzero_exponentmantissa(result);
238 *dstptr = result;
239 return(NOEXCEPTION);
240 }
241
242 /*
243 * Check for cases where opnd1/opnd2 < n
244 *
245 * In this case the result's sign will be opposite that of
246 * opnd1. The mantissa also needs some correction.
247 */
248 Sgl_leftshiftby1(opnd1);
249 if (Sgl_isgreaterthan(opnd1,opnd2)) {
250 Sgl_invert_sign(result);
251 Sgl_subtract((opnd2<<1),opnd1,opnd1);
252 }
253 /* check for remainder being exactly 1/2 of opnd2 */
254 else if (Sgl_isequal(opnd1,opnd2) && roundup) {
255 Sgl_invert_sign(result);
256 }
257
258 /* normalize result's mantissa */
259 while (Sgl_iszero_hidden(opnd1)) {
260 dest_exponent--;
261 Sgl_leftshiftby1(opnd1);
262 }
263 Sgl_set_exponentmantissa(result,opnd1);
264
265 /*
266 * Test for underflow
267 */
268 testforunderflow:
269 if (dest_exponent <= 0) {
270 /* trap if UNDERFLOWTRAP enabled */
271 if (Is_underflowtrap_enabled()) {
272 /*
273 * Adjust bias of result
274 */
275 Sgl_setwrapped_exponent(result,dest_exponent,unfl);
276 *dstptr = result;
277 /* frem is always exact */
278 return(UNDERFLOWEXCEPTION);
279 }
280 /*
281 * denormalize result or set to signed zero
282 */
283 if (dest_exponent >= (1 - SGL_P)) {
284 Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
285 } else {
286 Sgl_setzero_exponentmantissa(result);
287 }
288 }
289 else Sgl_set_exponent(result,dest_exponent);
290 *dstptr = result;
291 return(NOEXCEPTION);
292 }
293