arch/parisc/math-emu/sfadd.c

   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 denormalized 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     }