1 /* $OpenBSD: dfdiv.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 /* @(#)dfdiv.c: Revision: 2.11.88.1 Date: 93/12/07 15:05:39 */
16
17 #include "float.h"
18 #include "dbl_float.h"
19
20 /*
21 * Double Precision Floating-point Divide
22 */
23
24 int
dbl_fdiv(srcptr1,srcptr2,dstptr,status)25 dbl_fdiv(srcptr1,srcptr2,dstptr,status)
26 dbl_floating_point *srcptr1, *srcptr2, *dstptr;
27 unsigned int *status;
28 {
29 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
30 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
31 register int dest_exponent, count;
32 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
33 int is_tiny;
34
35 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
36 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
37 /*
38 * set sign bit of result
39 */
40 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
41 Dbl_setnegativezerop1(resultp1);
42 else Dbl_setzerop1(resultp1);
43 /*
44 * check first operand for NaN's or infinity
45 */
46 if (Dbl_isinfinity_exponent(opnd1p1)) {
47 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
48 if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
49 if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
50 /*
51 * invalid since both operands
52 * are infinity
53 */
54 if (Is_invalidtrap_enabled())
55 return(INVALIDEXCEPTION);
56 Set_invalidflag();
57 Dbl_makequietnan(resultp1,resultp2);
58 Dbl_copytoptr(resultp1,resultp2,dstptr);
59 return(NOEXCEPTION);
60 }
61 /*
62 * return infinity
63 */
64 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
65 Dbl_copytoptr(resultp1,resultp2,dstptr);
66 return(NOEXCEPTION);
67 }
68 }
69 else {
70 /*
71 * is NaN; signaling or quiet?
72 */
73 if (Dbl_isone_signaling(opnd1p1)) {
74 /* trap if INVALIDTRAP enabled */
75 if (Is_invalidtrap_enabled())
76 return(INVALIDEXCEPTION);
77 /* make NaN quiet */
78 Set_invalidflag();
79 Dbl_set_quiet(opnd1p1);
80 }
81 /*
82 * is second operand a signaling NaN?
83 */
84 else if (Dbl_is_signalingnan(opnd2p1)) {
85 /* trap if INVALIDTRAP enabled */
86 if (Is_invalidtrap_enabled())
87 return(INVALIDEXCEPTION);
88 /* make NaN quiet */
89 Set_invalidflag();
90 Dbl_set_quiet(opnd2p1);
91 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
92 return(NOEXCEPTION);
93 }
94 /*
95 * return quiet NaN
96 */
97 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
98 return(NOEXCEPTION);
99 }
100 }
101 /*
102 * check second operand for NaN's or infinity
103 */
104 if (Dbl_isinfinity_exponent(opnd2p1)) {
105 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
106 /*
107 * return zero
108 */
109 Dbl_setzero_exponentmantissa(resultp1,resultp2);
110 Dbl_copytoptr(resultp1,resultp2,dstptr);
111 return(NOEXCEPTION);
112 }
113 /*
114 * is NaN; signaling or quiet?
115 */
116 if (Dbl_isone_signaling(opnd2p1)) {
117 /* trap if INVALIDTRAP enabled */
118 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
119 /* make NaN quiet */
120 Set_invalidflag();
121 Dbl_set_quiet(opnd2p1);
122 }
123 /*
124 * return quiet NaN
125 */
126 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
127 return(NOEXCEPTION);
128 }
129 /*
130 * check for division by zero
131 */
132 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
133 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
134 /* invalid since both operands are zero */
135 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
136 Set_invalidflag();
137 Dbl_makequietnan(resultp1,resultp2);
138 Dbl_copytoptr(resultp1,resultp2,dstptr);
139 return(NOEXCEPTION);
140 }
141 if (Is_divisionbyzerotrap_enabled())
142 return(DIVISIONBYZEROEXCEPTION);
143 Set_divisionbyzeroflag();
144 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
145 Dbl_copytoptr(resultp1,resultp2,dstptr);
146 return(NOEXCEPTION);
147 }
148 /*
149 * Generate exponent
150 */
151 dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
152
153 /*
154 * Generate mantissa
155 */
156 if (Dbl_isnotzero_exponent(opnd1p1)) {
157 /* set hidden bit */
158 Dbl_clear_signexponent_set_hidden(opnd1p1);
159 }
160 else {
161 /* check for zero */
162 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
163 Dbl_setzero_exponentmantissa(resultp1,resultp2);
164 Dbl_copytoptr(resultp1,resultp2,dstptr);
165 return(NOEXCEPTION);
166 }
167 /* is denormalized, want to normalize */
168 Dbl_clear_signexponent(opnd1p1);
169 Dbl_leftshiftby1(opnd1p1,opnd1p2);
170 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
171 }
172 /* opnd2 needs to have hidden bit set with msb in hidden bit */
173 if (Dbl_isnotzero_exponent(opnd2p1)) {
174 Dbl_clear_signexponent_set_hidden(opnd2p1);
175 }
176 else {
177 /* is denormalized; want to normalize */
178 Dbl_clear_signexponent(opnd2p1);
179 Dbl_leftshiftby1(opnd2p1,opnd2p2);
180 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
181 dest_exponent+=8;
182 Dbl_leftshiftby8(opnd2p1,opnd2p2);
183 }
184 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
185 dest_exponent+=4;
186 Dbl_leftshiftby4(opnd2p1,opnd2p2);
187 }
188 while (Dbl_iszero_hidden(opnd2p1)) {
189 dest_exponent++;
190 Dbl_leftshiftby1(opnd2p1,opnd2p2);
191 }
192 }
193
194 /* Divide the source mantissas */
195
196 /*
197 * A non-restoring divide algorithm is used.
198 */
199 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
200 Dbl_setzero(opnd3p1,opnd3p2);
201 for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
202 Dbl_leftshiftby1(opnd1p1,opnd1p2);
203 Dbl_leftshiftby1(opnd3p1,opnd3p2);
204 if (Dbl_iszero_sign(opnd1p1)) {
205 Dbl_setone_lowmantissap2(opnd3p2);
206 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
207 }
208 else {
209 Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
210 }
211 }
212 if (count <= DBL_P) {
213 Dbl_leftshiftby1(opnd3p1,opnd3p2);
214 Dbl_setone_lowmantissap2(opnd3p2);
215 Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
216 if (Dbl_iszero_hidden(opnd3p1)) {
217 Dbl_leftshiftby1(opnd3p1,opnd3p2);
218 dest_exponent--;
219 }
220 }
221 else {
222 if (Dbl_iszero_hidden(opnd3p1)) {
223 /* need to get one more bit of result */
224 Dbl_leftshiftby1(opnd1p1,opnd1p2);
225 Dbl_leftshiftby1(opnd3p1,opnd3p2);
226 if (Dbl_iszero_sign(opnd1p1)) {
227 Dbl_setone_lowmantissap2(opnd3p2);
228 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
229 }
230 else {
231 Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
232 }
233 dest_exponent--;
234 }
235 if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
236 stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
237 }
238 inexact = guardbit | stickybit;
239
240 /*
241 * round result
242 */
243 if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
244 Dbl_clear_signexponent(opnd3p1);
245 switch (Rounding_mode()) {
246 case ROUNDPLUS:
247 if (Dbl_iszero_sign(resultp1))
248 Dbl_increment(opnd3p1,opnd3p2);
249 break;
250 case ROUNDMINUS:
251 if (Dbl_isone_sign(resultp1))
252 Dbl_increment(opnd3p1,opnd3p2);
253 break;
254 case ROUNDNEAREST:
255 if (guardbit && (stickybit ||
256 Dbl_isone_lowmantissap2(opnd3p2))) {
257 Dbl_increment(opnd3p1,opnd3p2);
258 }
259 }
260 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
261 }
262 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
263
264 /*
265 * Test for overflow
266 */
267 if (dest_exponent >= DBL_INFINITY_EXPONENT) {
268 /* trap if OVERFLOWTRAP enabled */
269 if (Is_overflowtrap_enabled()) {
270 /*
271 * Adjust bias of result
272 */
273 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
274 Dbl_copytoptr(resultp1,resultp2,dstptr);
275 if (inexact) {
276 if (Is_inexacttrap_enabled())
277 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
278 else
279 Set_inexactflag();
280 }
281 return(OVERFLOWEXCEPTION);
282 }
283 Set_overflowflag();
284 /* set result to infinity or largest number */
285 Dbl_setoverflow(resultp1,resultp2);
286 inexact = TRUE;
287 }
288 /*
289 * Test for underflow
290 */
291 else if (dest_exponent <= 0) {
292 /* trap if UNDERFLOWTRAP enabled */
293 if (Is_underflowtrap_enabled()) {
294 /*
295 * Adjust bias of result
296 */
297 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
298 Dbl_copytoptr(resultp1,resultp2,dstptr);
299 if (inexact) {
300 if (Is_inexacttrap_enabled())
301 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
302 else
303 Set_inexactflag();
304 }
305 return(UNDERFLOWEXCEPTION);
306 }
307
308 /* Determine if should set underflow flag */
309 is_tiny = TRUE;
310 if (dest_exponent == 0 && inexact) {
311 switch (Rounding_mode()) {
312 case ROUNDPLUS:
313 if (Dbl_iszero_sign(resultp1)) {
314 Dbl_increment(opnd3p1,opnd3p2);
315 if (Dbl_isone_hiddenoverflow(opnd3p1))
316 is_tiny = FALSE;
317 Dbl_decrement(opnd3p1,opnd3p2);
318 }
319 break;
320 case ROUNDMINUS:
321 if (Dbl_isone_sign(resultp1)) {
322 Dbl_increment(opnd3p1,opnd3p2);
323 if (Dbl_isone_hiddenoverflow(opnd3p1))
324 is_tiny = FALSE;
325 Dbl_decrement(opnd3p1,opnd3p2);
326 }
327 break;
328 case ROUNDNEAREST:
329 if (guardbit && (stickybit ||
330 Dbl_isone_lowmantissap2(opnd3p2))) {
331 Dbl_increment(opnd3p1,opnd3p2);
332 if (Dbl_isone_hiddenoverflow(opnd3p1))
333 is_tiny = FALSE;
334 Dbl_decrement(opnd3p1,opnd3p2);
335 }
336 break;
337 }
338 }
339
340 /*
341 * denormalize result or set to signed zero
342 */
343 stickybit = inexact;
344 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
345 stickybit,inexact);
346
347 /* return rounded number */
348 if (inexact) {
349 switch (Rounding_mode()) {
350 case ROUNDPLUS:
351 if (Dbl_iszero_sign(resultp1)) {
352 Dbl_increment(opnd3p1,opnd3p2);
353 }
354 break;
355 case ROUNDMINUS:
356 if (Dbl_isone_sign(resultp1)) {
357 Dbl_increment(opnd3p1,opnd3p2);
358 }
359 break;
360 case ROUNDNEAREST:
361 if (guardbit && (stickybit ||
362 Dbl_isone_lowmantissap2(opnd3p2))) {
363 Dbl_increment(opnd3p1,opnd3p2);
364 }
365 break;
366 }
367 if (is_tiny)
368 Set_underflowflag();
369 }
370 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
371 }
372 else Dbl_set_exponent(resultp1,dest_exponent);
373 Dbl_copytoptr(resultp1,resultp2,dstptr);
374
375 /* check for inexact */
376 if (inexact) {
377 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
378 else Set_inexactflag();
379 }
380 return(NOEXCEPTION);
381 }
382