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 denormalized 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 | } |
509 | |