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