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