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