![[BACK]](/icons/back.gif){kind=icon}
Return to dfadd.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath|
Return to dfadd.c CVS log | Up to [local] / sys / arch / hppa / spmath |
1.1 nbrk 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
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: }