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