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