1 /* $OpenBSD: sfdiv.c,v 1.5 2002/05/07 22:19:30 mickey Exp $ */
2 /*
3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4 To anyone who acknowledges that this file is provided "AS IS"
5 without any express or implied warranty:
6 permission to use, copy, modify, and distribute this file
7 for any purpose is hereby granted without fee, provided that
8 the above copyright notice and this notice appears in all
9 copies, and that the name of Hewlett-Packard Company not be
10 used in advertising or publicity pertaining to distribution
11 of the software without specific, written prior permission.
12 Hewlett-Packard Company makes no representations about the
13 suitability of this software for any purpose.
14 */
15 /* @(#)sfdiv.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:05 */
16
17 #include "float.h"
18 #include "sgl_float.h"
19
20 /*
21 * Single Precision Floating-point Divide
22 */
23 int
sgl_fdiv(srcptr1,srcptr2,dstptr,status)24 sgl_fdiv(srcptr1,srcptr2,dstptr,status)
25 sgl_floating_point *srcptr1, *srcptr2, *dstptr;
26 unsigned int *status;
27 {
28 register unsigned int opnd1, opnd2, opnd3, result;
29 register int dest_exponent, count;
30 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
31 int is_tiny;
32
33 opnd1 = *srcptr1;
34 opnd2 = *srcptr2;
35 /*
36 * set sign bit of result
37 */
38 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
39 else Sgl_setzero(result);
40 /*
41 * check first operand for NaN's or infinity
42 */
43 if (Sgl_isinfinity_exponent(opnd1)) {
44 if (Sgl_iszero_mantissa(opnd1)) {
45 if (Sgl_isnotnan(opnd2)) {
46 if (Sgl_isinfinity(opnd2)) {
47 /*
48 * invalid since both operands
49 * are infinity
50 */
51 if (Is_invalidtrap_enabled())
52 return(INVALIDEXCEPTION);
53 Set_invalidflag();
54 Sgl_makequietnan(result);
55 *dstptr = result;
56 return(NOEXCEPTION);
57 }
58 /*
59 * return infinity
60 */
61 Sgl_setinfinity_exponentmantissa(result);
62 *dstptr = result;
63 return(NOEXCEPTION);
64 }
65 }
66 else {
67 /*
68 * is NaN; signaling or quiet?
69 */
70 if (Sgl_isone_signaling(opnd1)) {
71 /* trap if INVALIDTRAP enabled */
72 if (Is_invalidtrap_enabled())
73 return(INVALIDEXCEPTION);
74 /* make NaN quiet */
75 Set_invalidflag();
76 Sgl_set_quiet(opnd1);
77 }
78 /*
79 * is second operand a signaling NaN?
80 */
81 else if (Sgl_is_signalingnan(opnd2)) {
82 /* trap if INVALIDTRAP enabled */
83 if (Is_invalidtrap_enabled())
84 return(INVALIDEXCEPTION);
85 /* make NaN quiet */
86 Set_invalidflag();
87 Sgl_set_quiet(opnd2);
88 *dstptr = opnd2;
89 return(NOEXCEPTION);
90 }
91 /*
92 * return quiet NaN
93 */
94 *dstptr = opnd1;
95 return(NOEXCEPTION);
96 }
97 }
98 /*
99 * check second operand for NaN's or infinity
100 */
101 if (Sgl_isinfinity_exponent(opnd2)) {
102 if (Sgl_iszero_mantissa(opnd2)) {
103 /*
104 * return zero
105 */
106 Sgl_setzero_exponentmantissa(result);
107 *dstptr = result;
108 return(NOEXCEPTION);
109 }
110 /*
111 * is NaN; signaling or quiet?
112 */
113 if (Sgl_isone_signaling(opnd2)) {
114 /* trap if INVALIDTRAP enabled */
115 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
116 /* make NaN quiet */
117 Set_invalidflag();
118 Sgl_set_quiet(opnd2);
119 }
120 /*
121 * return quiet NaN
122 */
123 *dstptr = opnd2;
124 return(NOEXCEPTION);
125 }
126 /*
127 * check for division by zero
128 */
129 if (Sgl_iszero_exponentmantissa(opnd2)) {
130 if (Sgl_iszero_exponentmantissa(opnd1)) {
131 /* invalid since both operands are zero */
132 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
133 Set_invalidflag();
134 Sgl_makequietnan(result);
135 *dstptr = result;
136 return(NOEXCEPTION);
137 }
138 if (Is_divisionbyzerotrap_enabled())
139 return(DIVISIONBYZEROEXCEPTION);
140 Set_divisionbyzeroflag();
141 Sgl_setinfinity_exponentmantissa(result);
142 *dstptr = result;
143 return(NOEXCEPTION);
144 }
145 /*
146 * Generate exponent
147 */
148 dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
149
150 /*
151 * Generate mantissa
152 */
153 if (Sgl_isnotzero_exponent(opnd1)) {
154 /* set hidden bit */
155 Sgl_clear_signexponent_set_hidden(opnd1);
156 }
157 else {
158 /* check for zero */
159 if (Sgl_iszero_mantissa(opnd1)) {
160 Sgl_setzero_exponentmantissa(result);
161 *dstptr = result;
162 return(NOEXCEPTION);
163 }
164 /* is denormalized; want to normalize */
165 Sgl_clear_signexponent(opnd1);
166 Sgl_leftshiftby1(opnd1);
167 Sgl_normalize(opnd1,dest_exponent);
168 }
169 /* opnd2 needs to have hidden bit set with msb in hidden bit */
170 if (Sgl_isnotzero_exponent(opnd2)) {
171 Sgl_clear_signexponent_set_hidden(opnd2);
172 }
173 else {
174 /* is denormalized; want to normalize */
175 Sgl_clear_signexponent(opnd2);
176 Sgl_leftshiftby1(opnd2);
177 while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
178 Sgl_leftshiftby8(opnd2);
179 dest_exponent += 8;
180 }
181 if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
182 Sgl_leftshiftby4(opnd2);
183 dest_exponent += 4;
184 }
185 while(Sgl_iszero_hidden(opnd2)) {
186 Sgl_leftshiftby1(opnd2);
187 dest_exponent += 1;
188 }
189 }
190
191 /* Divide the source mantissas */
192
193 /*
194 * A non_restoring divide algorithm is used.
195 */
196 Sgl_subtract(opnd1,opnd2,opnd1);
197 Sgl_setzero(opnd3);
198 for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
199 Sgl_leftshiftby1(opnd1);
200 Sgl_leftshiftby1(opnd3);
201 if (Sgl_iszero_sign(opnd1)) {
202 Sgl_setone_lowmantissa(opnd3);
203 Sgl_subtract(opnd1,opnd2,opnd1);
204 }
205 else Sgl_addition(opnd1,opnd2,opnd1);
206 }
207 if (count <= SGL_P) {
208 Sgl_leftshiftby1(opnd3);
209 Sgl_setone_lowmantissa(opnd3);
210 Sgl_leftshift(opnd3,SGL_P-count);
211 if (Sgl_iszero_hidden(opnd3)) {
212 Sgl_leftshiftby1(opnd3);
213 dest_exponent--;
214 }
215 }
216 else {
217 if (Sgl_iszero_hidden(opnd3)) {
218 /* need to get one more bit of result */
219 Sgl_leftshiftby1(opnd1);
220 Sgl_leftshiftby1(opnd3);
221 if (Sgl_iszero_sign(opnd1)) {
222 Sgl_setone_lowmantissa(opnd3);
223 Sgl_subtract(opnd1,opnd2,opnd1);
224 }
225 else Sgl_addition(opnd1,opnd2,opnd1);
226 dest_exponent--;
227 }
228 if (Sgl_iszero_sign(opnd1)) guardbit = TRUE;
229 stickybit = Sgl_all(opnd1);
230 }
231 inexact = guardbit | stickybit;
232
233 /*
234 * round result
235 */
236 if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
237 Sgl_clear_signexponent(opnd3);
238 switch (Rounding_mode()) {
239 case ROUNDPLUS:
240 if (Sgl_iszero_sign(result))
241 Sgl_increment_mantissa(opnd3);
242 break;
243 case ROUNDMINUS:
244 if (Sgl_isone_sign(result))
245 Sgl_increment_mantissa(opnd3);
246 break;
247 case ROUNDNEAREST:
248 if (guardbit &&
249 (stickybit || Sgl_isone_lowmantissa(opnd3)))
250 Sgl_increment_mantissa(opnd3);
251 }
252 if (Sgl_isone_hidden(opnd3)) dest_exponent++;
253 }
254 Sgl_set_mantissa(result,opnd3);
255
256 /*
257 * Test for overflow
258 */
259 if (dest_exponent >= SGL_INFINITY_EXPONENT) {
260 /* trap if OVERFLOWTRAP enabled */
261 if (Is_overflowtrap_enabled()) {
262 /*
263 * Adjust bias of result
264 */
265 Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
266 *dstptr = result;
267 if (inexact) {
268 if (Is_inexacttrap_enabled())
269 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
270 else Set_inexactflag();
271 }
272 return(OVERFLOWEXCEPTION);
273 }
274 Set_overflowflag();
275 /* set result to infinity or largest number */
276 Sgl_setoverflow(result);
277 inexact = TRUE;
278 }
279 /*
280 * Test for underflow
281 */
282 else if (dest_exponent <= 0) {
283 /* trap if UNDERFLOWTRAP enabled */
284 if (Is_underflowtrap_enabled()) {
285 /*
286 * Adjust bias of result
287 */
288 Sgl_setwrapped_exponent(result,dest_exponent,unfl);
289 *dstptr = result;
290 if (inexact) {
291 if (Is_inexacttrap_enabled())
292 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
293 else Set_inexactflag();
294 }
295 return(UNDERFLOWEXCEPTION);
296 }
297
298 /* Determine if should set underflow flag */
299 is_tiny = TRUE;
300 if (dest_exponent == 0 && inexact) {
301 switch (Rounding_mode()) {
302 case ROUNDPLUS:
303 if (Sgl_iszero_sign(result)) {
304 Sgl_increment(opnd3);
305 if (Sgl_isone_hiddenoverflow(opnd3))
306 is_tiny = FALSE;
307 Sgl_decrement(opnd3);
308 }
309 break;
310 case ROUNDMINUS:
311 if (Sgl_isone_sign(result)) {
312 Sgl_increment(opnd3);
313 if (Sgl_isone_hiddenoverflow(opnd3))
314 is_tiny = FALSE;
315 Sgl_decrement(opnd3);
316 }
317 break;
318 case ROUNDNEAREST:
319 if (guardbit && (stickybit ||
320 Sgl_isone_lowmantissa(opnd3))) {
321 Sgl_increment(opnd3);
322 if (Sgl_isone_hiddenoverflow(opnd3))
323 is_tiny = FALSE;
324 Sgl_decrement(opnd3);
325 }
326 break;
327 }
328 }
329
330 /*
331 * denormalize result or set to signed zero
332 */
333 stickybit = inexact;
334 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
335
336 /* return rounded number */
337 if (inexact) {
338 switch (Rounding_mode()) {
339 case ROUNDPLUS:
340 if (Sgl_iszero_sign(result)) {
341 Sgl_increment(opnd3);
342 }
343 break;
344 case ROUNDMINUS:
345 if (Sgl_isone_sign(result)) {
346 Sgl_increment(opnd3);
347 }
348 break;
349 case ROUNDNEAREST:
350 if (guardbit && (stickybit ||
351 Sgl_isone_lowmantissa(opnd3))) {
352 Sgl_increment(opnd3);
353 }
354 break;
355 }
356 if (is_tiny) Set_underflowflag();
357 }
358 Sgl_set_exponentmantissa(result,opnd3);
359 }
360 else Sgl_set_exponent(result,dest_exponent);
361 *dstptr = result;
362 /* check for inexact */
363 if (inexact) {
364 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
365 else Set_inexactflag();
366 }
367 return(NOEXCEPTION);
368 }
369