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