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