arch/parisc/math-emu/sfadd.c

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