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