1 /* $OpenBSD: dfsub.c,v 1.7 2002/11/29 09:27:34 deraadt 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 /* @(#)dfsub.c: Revision: 2.8.88.1 Date: 93/12/07 15:05:48 */
16
17 #include "float.h"
18 #include "dbl_float.h"
19
20 /*
21 * Double_subtract: subtract two double precision values.
22 */
23 int
dbl_fsub(leftptr,rightptr,dstptr,status)24 dbl_fsub(leftptr, rightptr, dstptr, status)
25 dbl_floating_point *leftptr, *rightptr, *dstptr;
26 unsigned int *status;
27 {
28 register unsigned int signless_upper_left, signless_upper_right, save;
29 register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
30 register unsigned int resultp1 = 0, resultp2 = 0;
31
32 register int result_exponent, right_exponent, diff_exponent;
33 register int sign_save, jumpsize;
34 register int inexact = FALSE, underflowtrap;
35
36 /* Create local copies of the numbers */
37 Dbl_copyfromptr(leftptr,leftp1,leftp2);
38 Dbl_copyfromptr(rightptr,rightp1,rightp2);
39
40 /* A zero "save" helps discover equal operands (for later), *
41 * and is used in swapping operands (if needed). */
42 Dbl_xortointp1(leftp1,rightp1,/*to*/save);
43
44 /*
45 * check first operand for NaN's or infinity
46 */
47 if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
48 {
49 if (Dbl_iszero_mantissa(leftp1,leftp2))
50 {
51 if (Dbl_isnotnan(rightp1,rightp2))
52 {
53 if (Dbl_isinfinity(rightp1,rightp2) && save==0)
54 {
55 /*
56 * invalid since operands are same signed infinity's
57 */
58 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
59 Set_invalidflag();
60 Dbl_makequietnan(resultp1,resultp2);
61 Dbl_copytoptr(resultp1,resultp2,dstptr);
62 return(NOEXCEPTION);
63 }
64 /*
65 * return infinity
66 */
67 Dbl_copytoptr(leftp1,leftp2,dstptr);
68 return(NOEXCEPTION);
69 }
70 }
71 else
72 {
73 /*
74 * is NaN; signaling or quiet?
75 */
76 if (Dbl_isone_signaling(leftp1))
77 {
78 /* trap if INVALIDTRAP enabled */
79 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
80 /* make NaN quiet */
81 Set_invalidflag();
82 Dbl_set_quiet(leftp1);
83 }
84 /*
85 * is second operand a signaling NaN?
86 */
87 else if (Dbl_is_signalingnan(rightp1))
88 {
89 /* trap if INVALIDTRAP enabled */
90 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
91 /* make NaN quiet */
92 Set_invalidflag();
93 Dbl_set_quiet(rightp1);
94 Dbl_copytoptr(rightp1,rightp2,dstptr);
95 return(NOEXCEPTION);
96 }
97 /*
98 * return quiet NaN
99 */
100 Dbl_copytoptr(leftp1,leftp2,dstptr);
101 return(NOEXCEPTION);
102 }
103 } /* End left NaN or Infinity processing */
104 /*
105 * check second operand for NaN's or infinity
106 */
107 if (Dbl_isinfinity_exponent(rightp1))
108 {
109 if (Dbl_iszero_mantissa(rightp1,rightp2))
110 {
111 /* return infinity */
112 Dbl_invert_sign(rightp1);
113 Dbl_copytoptr(rightp1,rightp2,dstptr);
114 return(NOEXCEPTION);
115 }
116 /*
117 * is NaN; signaling or quiet?
118 */
119 if (Dbl_isone_signaling(rightp1))
120 {
121 /* trap if INVALIDTRAP enabled */
122 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
123 /* make NaN quiet */
124 Set_invalidflag();
125 Dbl_set_quiet(rightp1);
126 }
127 /*
128 * return quiet NaN
129 */
130 Dbl_copytoptr(rightp1,rightp2,dstptr);
131 return(NOEXCEPTION);
132 } /* End right NaN or Infinity processing */
133
134 /* Invariant: Must be dealing with finite numbers */
135
136 /* Compare operands by removing the sign */
137 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
138 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);
139
140 /* sign difference selects add or sub operation. */
141 if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
142 {
143 /* Set the left operand to the larger one by XOR swap *
144 * First finish the first word using "save" */
145 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
146 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
147 Dbl_swap_lower(leftp2,rightp2);
148 result_exponent = Dbl_exponent(leftp1);
149 Dbl_invert_sign(leftp1);
150 }
151 /* Invariant: left is not smaller than right. */
152
153 if((right_exponent = Dbl_exponent(rightp1)) == 0)
154 {
155 /* Denormalized operands. First look for zeroes */
156 if(Dbl_iszero_mantissa(rightp1,rightp2))
157 {
158 /* right is zero */
159 if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
160 {
161 /* Both operands are zeros */
162 Dbl_invert_sign(rightp1);
163 if(Is_rounding_mode(ROUNDMINUS))
164 {
165 Dbl_or_signs(leftp1,/*with*/rightp1);
166 }
167 else
168 {
169 Dbl_and_signs(leftp1,/*with*/rightp1);
170 }
171 }
172 else
173 {
174 /* Left is not a zero and must be the result. Trapped
175 * underflows are signaled if left is denormalized. Result
176 * is always exact. */
177 if( (result_exponent == 0) && Is_underflowtrap_enabled() )
178 {
179 /* need to normalize results mantissa */
180 sign_save = Dbl_signextendedsign(leftp1);
181 Dbl_leftshiftby1(leftp1,leftp2);
182 Dbl_normalize(leftp1,leftp2,result_exponent);
183 Dbl_set_sign(leftp1,/*using*/sign_save);
184 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
185 Dbl_copytoptr(leftp1,leftp2,dstptr);
186 /* inexact = FALSE */
187 return(UNDERFLOWEXCEPTION);
188 }
189 }
190 Dbl_copytoptr(leftp1,leftp2,dstptr);
191 return(NOEXCEPTION);
192 }
193
194 /* Neither are zeroes */
195 Dbl_clear_sign(rightp1); /* Exponent is already cleared */
196 if(result_exponent == 0 )
197 {
198 /* Both operands are denormalized. The result must be exact
199 * and is simply calculated. A sum could become normalized and a
200 * difference could cancel to a true zero. */
201 if( (/*signed*/int) save >= 0 )
202 {
203 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
204 /*into*/resultp1,resultp2);
205 if(Dbl_iszero_mantissa(resultp1,resultp2))
206 {
207 if(Is_rounding_mode(ROUNDMINUS))
208 {
209 Dbl_setone_sign(resultp1);
210 }
211 else
212 {
213 Dbl_setzero_sign(resultp1);
214 }
215 Dbl_copytoptr(resultp1,resultp2,dstptr);
216 return(NOEXCEPTION);
217 }
218 }
219 else
220 {
221 Dbl_addition(leftp1,leftp2,rightp1,rightp2,
222 /*into*/resultp1,resultp2);
223 if(Dbl_isone_hidden(resultp1))
224 {
225 Dbl_copytoptr(resultp1,resultp2,dstptr);
226 return(NOEXCEPTION);
227 }
228 }
229 if(Is_underflowtrap_enabled())
230 {
231 /* need to normalize result */
232 sign_save = Dbl_signextendedsign(resultp1);
233 Dbl_leftshiftby1(resultp1,resultp2);
234 Dbl_normalize(resultp1,resultp2,result_exponent);
235 Dbl_set_sign(resultp1,/*using*/sign_save);
236 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
237 Dbl_copytoptr(resultp1,resultp2,dstptr);
238 /* inexact = FALSE */
239 return(UNDERFLOWEXCEPTION);
240 }
241 Dbl_copytoptr(resultp1,resultp2,dstptr);
242 return(NOEXCEPTION);
243 }
244 right_exponent = 1; /* Set exponent to reflect different bias
245 * with denomalized numbers. */
246 }
247 else
248 {
249 Dbl_clear_signexponent_set_hidden(rightp1);
250 }
251 Dbl_clear_exponent_set_hidden(leftp1);
252 diff_exponent = result_exponent - right_exponent;
253
254 /*
255 * Special case alignment of operands that would force alignment
256 * beyond the extent of the extension. A further optimization
257 * could special case this but only reduces the path length for this
258 * infrequent case.
259 */
260 if(diff_exponent > DBL_THRESHOLD)
261 {
262 diff_exponent = DBL_THRESHOLD;
263 }
264
265 /* Align right operand by shifting to right */
266 Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
267 /*and lower to*/extent);
268
269 /* Treat sum and difference of the operands separately. */
270 if( (/*signed*/int) save >= 0 )
271 {
272 /*
273 * Difference of the two operands. Their can be no overflow. A
274 * borrow can occur out of the hidden bit and force a post
275 * normalization phase.
276 */
277 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
278 /*with*/extent,/*into*/resultp1,resultp2);
279 if(Dbl_iszero_hidden(resultp1))
280 {
281 /* Handle normalization */
282 /* A straight forward algorithm would now shift the result
283 * and extension left until the hidden bit becomes one. Not
284 * all of the extension bits need participate in the shift.
285 * Only the two most significant bits (round and guard) are
286 * needed. If only a single shift is needed then the guard
287 * bit becomes a significant low order bit and the extension
288 * must participate in the rounding. If more than a single
289 * shift is needed, then all bits to the right of the guard
290 * bit are zeros, and the guard bit may or may not be zero. */
291 sign_save = Dbl_signextendedsign(resultp1);
292 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);
293
294 /* Need to check for a zero result. The sign and exponent
295 * fields have already been zeroed. The more efficient test
296 * of the full object can be used.
297 */
298 if(Dbl_iszero(resultp1,resultp2))
299 /* Must have been "x-x" or "x+(-x)". */
300 {
301 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
302 Dbl_copytoptr(resultp1,resultp2,dstptr);
303 return(NOEXCEPTION);
304 }
305 result_exponent--;
306 /* Look to see if normalization is finished. */
307 if(Dbl_isone_hidden(resultp1)) {
308 if(result_exponent==0) {
309 /* Denormalized, exponent should be zero. Left operand *
310 * was normalized, so extent (guard, round) was zero */
311 goto underflow;
312 } else {
313 /* No further normalization is needed. */
314 Dbl_set_sign(resultp1,/*using*/sign_save);
315 Ext_leftshiftby1(extent);
316 goto round;
317 }
318 }
319
320 /* Check for denormalized, exponent should be zero. Left *
321 * operand was normalized, so extent (guard, round) was zero */
322 if(!(underflowtrap = Is_underflowtrap_enabled()) &&
323 result_exponent==0) goto underflow;
324
325 /* Shift extension to complete one bit of normalization and
326 * update exponent. */
327 Ext_leftshiftby1(extent);
328
329 /* Discover first one bit to determine shift amount. Use a
330 * modified binary search. We have already shifted the result
331 * one position right and still not found a one so the remainder
332 * of the extension must be zero and simplifies rounding. */
333 /* Scan bytes */
334 while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
335 {
336 Dbl_leftshiftby8(resultp1,resultp2);
337 if((result_exponent -= 8) <= 0 && !underflowtrap)
338 goto underflow;
339 }
340 /* Now narrow it down to the nibble */
341 if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
342 {
343 /* The lower nibble contains the normalizing one */
344 Dbl_leftshiftby4(resultp1,resultp2);
345 if((result_exponent -= 4) <= 0 && !underflowtrap)
346 goto underflow;
347 }
348 /* Select case were first bit is set (already normalized)
349 * otherwise select the proper shift. */
350 if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
351 {
352 /* Already normalized */
353 if(result_exponent <= 0) goto underflow;
354 Dbl_set_sign(resultp1,/*using*/sign_save);
355 Dbl_set_exponent(resultp1,/*using*/result_exponent);
356 Dbl_copytoptr(resultp1,resultp2,dstptr);
357 return(NOEXCEPTION);
358 }
359 Dbl_sethigh4bits(resultp1,/*using*/sign_save);
360 switch(jumpsize)
361 {
362 case 1:
363 {
364 Dbl_leftshiftby3(resultp1,resultp2);
365 result_exponent -= 3;
366 break;
367 }
368 case 2:
369 case 3:
370 {
371 Dbl_leftshiftby2(resultp1,resultp2);
372 result_exponent -= 2;
373 break;
374 }
375 case 4:
376 case 5:
377 case 6:
378 case 7:
379 {
380 Dbl_leftshiftby1(resultp1,resultp2);
381 result_exponent -= 1;
382 break;
383 }
384 }
385 if(result_exponent > 0)
386 {
387 Dbl_set_exponent(resultp1,/*using*/result_exponent);
388 Dbl_copytoptr(resultp1,resultp2,dstptr);
389 return(NOEXCEPTION); /* Sign bit is already set */
390 }
391 /* Fixup potential underflows */
392 underflow:
393 if(Is_underflowtrap_enabled())
394 {
395 Dbl_set_sign(resultp1,sign_save);
396 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
397 Dbl_copytoptr(resultp1,resultp2,dstptr);
398 /* inexact = FALSE */
399 return(UNDERFLOWEXCEPTION);
400 }
401 /*
402 * Since we cannot get an inexact denormalized result,
403 * we can now return.
404 */
405 Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
406 Dbl_clear_signexponent(resultp1);
407 Dbl_set_sign(resultp1,sign_save);
408 Dbl_copytoptr(resultp1,resultp2,dstptr);
409 return(NOEXCEPTION);
410 } /* end if(hidden...)... */
411 /* Fall through and round */
412 } /* end if(save >= 0)... */
413 else
414 {
415 /* Subtract magnitudes */
416 Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
417 if(Dbl_isone_hiddenoverflow(resultp1))
418 {
419 /* Prenormalization required. */
420 Dbl_rightshiftby1_withextent(resultp2,extent,extent);
421 Dbl_arithrightshiftby1(resultp1,resultp2);
422 result_exponent++;
423 } /* end if hiddenoverflow... */
424 } /* end else ...subtract magnitudes... */
425
426 /* Round the result. If the extension is all zeros,then the result is
427 * exact. Otherwise round in the correct direction. No underflow is
428 * possible. If a postnormalization is necessary, then the mantissa is
429 * all zeros so no shift is needed. */
430 round:
431 if(Ext_isnotzero(extent))
432 {
433 inexact = TRUE;
434 switch(Rounding_mode())
435 {
436 case ROUNDNEAREST: /* The default. */
437 if(Ext_isone_sign(extent))
438 {
439 /* at least 1/2 ulp */
440 if(Ext_isnotzero_lower(extent) ||
441 Dbl_isone_lowmantissap2(resultp2))
442 {
443 /* either exactly half way and odd or more than 1/2ulp */
444 Dbl_increment(resultp1,resultp2);
445 }
446 }
447 break;
448
449 case ROUNDPLUS:
450 if(Dbl_iszero_sign(resultp1))
451 {
452 /* Round up positive results */
453 Dbl_increment(resultp1,resultp2);
454 }
455 break;
456
457 case ROUNDMINUS:
458 if(Dbl_isone_sign(resultp1))
459 {
460 /* Round down negative results */
461 Dbl_increment(resultp1,resultp2);
462 }
463
464 case ROUNDZERO:;
465 /* truncate is simple */
466 } /* end switch... */
467 if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
468 }
469 if(result_exponent == DBL_INFINITY_EXPONENT)
470 {
471 /* Overflow */
472 if(Is_overflowtrap_enabled())
473 {
474 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
475 Dbl_copytoptr(resultp1,resultp2,dstptr);
476 if (inexact) {
477 if (Is_inexacttrap_enabled())
478 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
479 else
480 Set_inexactflag();
481 }
482 return(OVERFLOWEXCEPTION);
483 }
484 else
485 {
486 inexact = TRUE;
487 Set_overflowflag();
488 Dbl_setoverflow(resultp1,resultp2);
489 }
490 }
491 else Dbl_set_exponent(resultp1,result_exponent);
492 Dbl_copytoptr(resultp1,resultp2,dstptr);
493 if(inexact) {
494 if(Is_inexacttrap_enabled())
495 return(INEXACTEXCEPTION);
496 else
497 Set_inexactflag();
498 }
499 return(NOEXCEPTION);
500 }
501