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