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