1/* $NetBSD: dfadd.c,v 1.7 2023/07/31 20:48:04 andvar Exp $ */ 2 3/* $OpenBSD: dfadd.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 * pmk1.1 27 */ 28/* 29 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY 30 * 31 * To anyone who acknowledges that this file is provided "AS IS" 32 * without any express or implied warranty: 33 * permission to use, copy, modify, and distribute this file 34 * for any purpose is hereby granted without fee, provided that 35 * the above copyright notice and this notice appears in all 36 * copies, and that the name of Hewlett-Packard Company not be 37 * used in advertising or publicity pertaining to distribution 38 * of the software without specific, written prior permission. 39 * Hewlett-Packard Company makes no representations about the 40 * suitability of this software for any purpose. 41 */ 42 43#include <sys/cdefs.h> 44__KERNEL_RCSID(0, "$NetBSD: dfadd.c,v 1.7 2023/07/31 20:48:04 andvar Exp $"); 45 46#include "../spmath/float.h" 47#include "../spmath/dbl_float.h" 48 49/* 50 * Double_add: add two double precision values. 51 */ 52int 53dbl_fadd(dbl_floating_point *leftptr, dbl_floating_point *rightptr, 54 dbl_floating_point *dstptr, unsigned int *status) 55{ 56 register unsigned int signless_upper_left, signless_upper_right, save; 57 register unsigned int leftp1, leftp2, rightp1, rightp2, extent; 58 register unsigned int resultp1 = 0, resultp2 = 0; 59 60 register int result_exponent, right_exponent, diff_exponent; 61 register int sign_save, jumpsize; 62 register int inexact = false; 63 register int underflowtrap; 64 65 /* Create local copies of the numbers */ 66 Dbl_copyfromptr(leftptr,leftp1,leftp2); 67 Dbl_copyfromptr(rightptr,rightp1,rightp2); 68 69 /* A zero "save" helps discover equal operands (for later), * 70 * and is used in swapping operands (if needed). */ 71 Dbl_xortointp1(leftp1,rightp1,/*to*/save); 72 73 /* 74 * check first operand for NaN's or infinity 75 */ 76 if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT) 77 { 78 if (Dbl_iszero_mantissa(leftp1,leftp2)) 79 { 80 if (Dbl_isnotnan(rightp1,rightp2)) 81 { 82 if (Dbl_isinfinity(rightp1,rightp2) && save!=0) 83 { 84 /* 85 * invalid since operands are opposite signed infinity's 86 */ 87 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 88 Set_invalidflag(); 89 Dbl_makequietnan(resultp1,resultp2); 90 Dbl_copytoptr(resultp1,resultp2,dstptr); 91 return(NOEXCEPTION); 92 } 93 /* 94 * return infinity 95 */ 96 Dbl_copytoptr(leftp1,leftp2,dstptr); 97 return(NOEXCEPTION); 98 } 99 } 100 else 101 { 102 /* 103 * is NaN; signaling or quiet? 104 */ 105 if (Dbl_isone_signaling(leftp1)) 106 { 107 /* trap if INVALIDTRAP enabled */ 108 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 109 /* make NaN quiet */ 110 Set_invalidflag(); 111 Dbl_set_quiet(leftp1); 112 } 113 /* 114 * is second operand a signaling NaN? 115 */ 116 else if (Dbl_is_signalingnan(rightp1)) 117 { 118 /* trap if INVALIDTRAP enabled */ 119 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 120 /* make NaN quiet */ 121 Set_invalidflag(); 122 Dbl_set_quiet(rightp1); 123 Dbl_copytoptr(rightp1,rightp2,dstptr); 124 return(NOEXCEPTION); 125 } 126 /* 127 * return quiet NaN 128 */ 129 Dbl_copytoptr(leftp1,leftp2,dstptr); 130 return(NOEXCEPTION); 131 } 132 } /* End left NaN or Infinity processing */ 133 /* 134 * check second operand for NaN's or infinity 135 */ 136 if (Dbl_isinfinity_exponent(rightp1)) 137 { 138 if (Dbl_iszero_mantissa(rightp1,rightp2)) 139 { 140 /* return infinity */ 141 Dbl_copytoptr(rightp1,rightp2,dstptr); 142 return(NOEXCEPTION); 143 } 144 /* 145 * is NaN; signaling or quiet? 146 */ 147 if (Dbl_isone_signaling(rightp1)) 148 { 149 /* trap if INVALIDTRAP enabled */ 150 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 151 /* make NaN quiet */ 152 Set_invalidflag(); 153 Dbl_set_quiet(rightp1); 154 } 155 /* 156 * return quiet NaN 157 */ 158 Dbl_copytoptr(rightp1,rightp2,dstptr); 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 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left); 166 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right); 167 168 /* sign difference selects add or sub operation. */ 169 if(Dbl_ismagnitudeless(leftp2,rightp2,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 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1); 174 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1); 175 Dbl_swap_lower(leftp2,rightp2); 176 result_exponent = Dbl_exponent(leftp1); 177 } 178 /* Invariant: left is not smaller than right. */ 179 180 if((right_exponent = Dbl_exponent(rightp1)) == 0) 181 { 182 /* Denormalized operands. First look for zeroes */ 183 if(Dbl_iszero_mantissa(rightp1,rightp2)) 184 { 185 /* right is zero */ 186 if(Dbl_iszero_exponentmantissa(leftp1,leftp2)) 187 { 188 /* Both operands are zeros */ 189 if(Is_rounding_mode(ROUNDMINUS)) 190 { 191 Dbl_or_signs(leftp1,/*with*/rightp1); 192 } 193 else 194 { 195 Dbl_and_signs(leftp1,/*with*/rightp1); 196 } 197 } 198 else 199 { 200 /* Left is not a zero and must be the result. Trapped 201 * underflows are signaled if left is denormalized. Result 202 * is always exact. */ 203 if( (result_exponent == 0) && Is_underflowtrap_enabled() ) 204 { 205 /* need to normalize results mantissa */ 206 sign_save = Dbl_signextendedsign(leftp1); 207 Dbl_leftshiftby1(leftp1,leftp2); 208 Dbl_normalize(leftp1,leftp2,result_exponent); 209 Dbl_set_sign(leftp1,/*using*/sign_save); 210 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl); 211 Dbl_copytoptr(leftp1,leftp2,dstptr); 212 /* inexact = false */ 213 return(UNDERFLOWEXCEPTION); 214 } 215 } 216 Dbl_copytoptr(leftp1,leftp2,dstptr); 217 return(NOEXCEPTION); 218 } 219 220 /* Neither are zeroes */ 221 Dbl_clear_sign(rightp1); /* Exponent is already cleared */ 222 if(result_exponent == 0 ) 223 { 224 /* Both operands are denormalized. The result must be exact 225 * and is simply calculated. A sum could become normalized and a 226 * difference could cancel to a true zero. */ 227 if( (/*signed*/int) save < 0 ) 228 { 229 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2, 230 /*into*/resultp1,resultp2); 231 if(Dbl_iszero_mantissa(resultp1,resultp2)) 232 { 233 if(Is_rounding_mode(ROUNDMINUS)) 234 { 235 Dbl_setone_sign(resultp1); 236 } 237 else 238 { 239 Dbl_setzero_sign(resultp1); 240 } 241 Dbl_copytoptr(resultp1,resultp2,dstptr); 242 return(NOEXCEPTION); 243 } 244 } 245 else 246 { 247 Dbl_addition(leftp1,leftp2,rightp1,rightp2, 248 /*into*/resultp1,resultp2); 249 if(Dbl_isone_hidden(resultp1)) 250 { 251 Dbl_copytoptr(resultp1,resultp2,dstptr); 252 return(NOEXCEPTION); 253 } 254 } 255 if(Is_underflowtrap_enabled()) 256 { 257 /* need to normalize result */ 258 sign_save = Dbl_signextendedsign(resultp1); 259 Dbl_leftshiftby1(resultp1,resultp2); 260 Dbl_normalize(resultp1,resultp2,result_exponent); 261 Dbl_set_sign(resultp1,/*using*/sign_save); 262 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl); 263 Dbl_copytoptr(resultp1,resultp2,dstptr); 264 /* inexact = false */ 265 return(UNDERFLOWEXCEPTION); 266 } 267 Dbl_copytoptr(resultp1,resultp2,dstptr); 268 return(NOEXCEPTION); 269 } 270 right_exponent = 1; /* Set exponent to reflect different bias 271 * with denormalized numbers. */ 272 } 273 else 274 { 275 Dbl_clear_signexponent_set_hidden(rightp1); 276 } 277 Dbl_clear_exponent_set_hidden(leftp1); 278 diff_exponent = result_exponent - right_exponent; 279 280 /* 281 * Special case alignment of operands that would force alignment 282 * beyond the extent of the extension. A further optimization 283 * could special case this but only reduces the path length for this 284 * infrequent case. 285 */ 286 if(diff_exponent > DBL_THRESHOLD) 287 { 288 diff_exponent = DBL_THRESHOLD; 289 } 290 291 /* Align right operand by shifting to right */ 292 Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent, 293 /*and lower to*/extent); 294 295 /* Treat sum and difference of the operands separately. */ 296 if( (/*signed*/int) save < 0 ) 297 { 298 /* 299 * Difference of the two operands. Their can be no overflow. A 300 * borrow can occur out of the hidden bit and force a post 301 * normalization phase. 302 */ 303 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2, 304 /*with*/extent,/*into*/resultp1,resultp2); 305 if(Dbl_iszero_hidden(resultp1)) 306 { 307 /* Handle normalization */ 308 /* A straight forward algorithm would now shift the result 309 * and extension left until the hidden bit becomes one. Not 310 * all of the extension bits need participate in the shift. 311 * Only the two most significant bits (round and guard) are 312 * needed. If only a single shift is needed then the guard 313 * bit becomes a significant low order bit and the extension 314 * must participate in the rounding. If more than a single 315 * shift is needed, then all bits to the right of the guard 316 * bit are zeros, and the guard bit may or may not be zero. */ 317 sign_save = Dbl_signextendedsign(resultp1); 318 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2); 319 320 /* Need to check for a zero result. The sign and exponent 321 * fields have already been zeroed. The more efficient test 322 * of the full object can be used. 323 */ 324 if(Dbl_iszero(resultp1,resultp2)) 325 /* Must have been "x-x" or "x+(-x)". */ 326 { 327 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1); 328 Dbl_copytoptr(resultp1,resultp2,dstptr); 329 return(NOEXCEPTION); 330 } 331 result_exponent--; 332 /* Look to see if normalization is finished. */ 333 if(Dbl_isone_hidden(resultp1)) 334 { 335 if(result_exponent==0) 336 { 337 /* Denormalized, exponent should be zero. Left operand * 338 * was normalized, so extent (guard, round) was zero */ 339 goto underflow; 340 } 341 else 342 { 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 /* Add 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 ...add 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