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