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