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