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