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