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
sgl_fsub(leftptr,rightptr,dstptr,status)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 }
497