sys/arch/hppa/spmath/dfsub.c

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