sys/arch/hppa/spmath/dfdiv.c

   1 /*      $NetBSD: dfdiv.c,v 1.3 2005/12/11 12:17:40 christos Exp $       */
   2
   3 /*      $OpenBSD: dfdiv.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: dfdiv.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 Precision Floating-point Divide
  52  */
  53
  54 int
  55 dbl_fdiv(srcptr1,srcptr2,dstptr,status)
  56
  57 dbl_floating_point *srcptr1, *srcptr2, *dstptr;
  58 unsigned int *status;
  59 {
  60         register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
  61         register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
  62         register int dest_exponent, count;
  63         register int inexact = false, guardbit = false, stickybit = false;
  64         int is_tiny;
  65
  66         Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
  67         Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
  68         /*
  69          * set sign bit of result
  70          */
  71         if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
  72                 Dbl_setnegativezerop1(resultp1);
  73         else Dbl_setzerop1(resultp1);
  74         /*
  75          * check first operand for NaN's or infinity
  76          */
  77         if (Dbl_isinfinity_exponent(opnd1p1)) {
  78                 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
  79                         if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
  80                                 if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
  81                                         /*
  82                                          * invalid since both operands
  83                                          * are infinity
  84                                          */
  85                                         if (Is_invalidtrap_enabled())
  86                                                 return(INVALIDEXCEPTION);
  87                                         Set_invalidflag();
  88                                         Dbl_makequietnan(resultp1,resultp2);
  89                                         Dbl_copytoptr(resultp1,resultp2,dstptr);
  90                                         return(NOEXCEPTION);
  91                                 }
  92                                 /*
  93                                  * return infinity
  94                                  */
  95                                 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
  96                                 Dbl_copytoptr(resultp1,resultp2,dstptr);
  97                                 return(NOEXCEPTION);
  98                         }
  99                 }
 100                 else {
 101                         /*
 102                          * is NaN; signaling or quiet?
 103                          */
 104                         if (Dbl_isone_signaling(opnd1p1)) {
 105                                 /* trap if INVALIDTRAP enabled */
 106                                 if (Is_invalidtrap_enabled())
 107                                         return(INVALIDEXCEPTION);
 108                                 /* make NaN quiet */
 109                                 Set_invalidflag();
 110                                 Dbl_set_quiet(opnd1p1);
 111                         }
 112                         /*
 113                          * is second operand a signaling NaN?
 114                          */
 115                         else if (Dbl_is_signalingnan(opnd2p1)) {
 116                                 /* trap if INVALIDTRAP enabled */
 117                                 if (Is_invalidtrap_enabled())
 118                                         return(INVALIDEXCEPTION);
 119                                 /* make NaN quiet */
 120                                 Set_invalidflag();
 121                                 Dbl_set_quiet(opnd2p1);
 122                                 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
 123                                 return(NOEXCEPTION);
 124                         }
 125                         /*
 126                          * return quiet NaN
 127                          */
 128                         Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
 129                         return(NOEXCEPTION);
 130                 }
 131         }
 132         /*
 133          * check second operand for NaN's or infinity
 134          */
 135         if (Dbl_isinfinity_exponent(opnd2p1)) {
 136                 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
 137                         /*
 138                          * return zero
 139                          */
 140                         Dbl_setzero_exponentmantissa(resultp1,resultp2);
 141                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 142                         return(NOEXCEPTION);
 143                 }
 144                 /*
 145                  * is NaN; signaling or quiet?
 146                  */
 147                 if (Dbl_isone_signaling(opnd2p1)) {
 148                         /* trap if INVALIDTRAP enabled */
 149                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 150                         /* make NaN quiet */
 151                         Set_invalidflag();
 152                         Dbl_set_quiet(opnd2p1);
 153                 }
 154                 /*
 155                  * return quiet NaN
 156                  */
 157                 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
 158                 return(NOEXCEPTION);
 159         }
 160         /*
 161          * check for division by zero
 162          */
 163         if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
 164                 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
 165                         /* invalid since both operands are zero */
 166                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 167                         Set_invalidflag();
 168                         Dbl_makequietnan(resultp1,resultp2);
 169                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 170                         return(NOEXCEPTION);
 171                 }
 172                 if (Is_divisionbyzerotrap_enabled())
 173                         return(DIVISIONBYZEROEXCEPTION);
 174                 Set_divisionbyzeroflag();
 175                 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
 176                 Dbl_copytoptr(resultp1,resultp2,dstptr);
 177                 return(NOEXCEPTION);
 178         }
 179         /*
 180          * Generate exponent
 181          */
 182         dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
 183
 184         /*
 185          * Generate mantissa
 186          */
 187         if (Dbl_isnotzero_exponent(opnd1p1)) {
 188                 /* set hidden bit */
 189                 Dbl_clear_signexponent_set_hidden(opnd1p1);
 190         }
 191         else {
 192                 /* check for zero */
 193                 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
 194                         Dbl_setzero_exponentmantissa(resultp1,resultp2);
 195                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 196                         return(NOEXCEPTION);
 197                 }
 198                 /* is denormalized, want to normalize */
 199                 Dbl_clear_signexponent(opnd1p1);
 200                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
 201                 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
 202         }
 203         /* opnd2 needs to have hidden bit set with msb in hidden bit */
 204         if (Dbl_isnotzero_exponent(opnd2p1)) {
 205                 Dbl_clear_signexponent_set_hidden(opnd2p1);
 206         }
 207         else {
 208                 /* is denormalized; want to normalize */
 209                 Dbl_clear_signexponent(opnd2p1);
 210                 Dbl_leftshiftby1(opnd2p1,opnd2p2);
 211                 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
 212                         dest_exponent+=8;
 213                         Dbl_leftshiftby8(opnd2p1,opnd2p2);
 214                 }
 215                 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
 216                         dest_exponent+=4;
 217                         Dbl_leftshiftby4(opnd2p1,opnd2p2);
 218                 }
 219                 while (Dbl_iszero_hidden(opnd2p1)) {
 220                         dest_exponent++;
 221                         Dbl_leftshiftby1(opnd2p1,opnd2p2);
 222                 }
 223         }
 224
 225         /* Divide the source mantissas */
 226
 227         /*
 228          * A non-restoring divide algorithm is used.
 229          */
 230         Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
 231         Dbl_setzero(opnd3p1,opnd3p2);
 232         for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
 233                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
 234                 Dbl_leftshiftby1(opnd3p1,opnd3p2);
 235                 if (Dbl_iszero_sign(opnd1p1)) {
 236                         Dbl_setone_lowmantissap2(opnd3p2);
 237                         Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
 238                 }
 239                 else {
 240                         Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
 241                 }
 242         }
 243         if (count <= DBL_P) {
 244                 Dbl_leftshiftby1(opnd3p1,opnd3p2);
 245                 Dbl_setone_lowmantissap2(opnd3p2);
 246                 Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
 247                 if (Dbl_iszero_hidden(opnd3p1)) {
 248                         Dbl_leftshiftby1(opnd3p1,opnd3p2);
 249                         dest_exponent--;
 250                 }
 251         }
 252         else {
 253                 if (Dbl_iszero_hidden(opnd3p1)) {
 254                         /* need to get one more bit of result */
 255                         Dbl_leftshiftby1(opnd1p1,opnd1p2);
 256                         Dbl_leftshiftby1(opnd3p1,opnd3p2);
 257                         if (Dbl_iszero_sign(opnd1p1)) {
 258                                 Dbl_setone_lowmantissap2(opnd3p2);
 259                                 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
 260                         }
 261                         else {
 262                                 Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
 263                         }
 264                         dest_exponent--;
 265                 }
 266                 if (Dbl_iszero_sign(opnd1p1)) guardbit = true;
 267                 stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
 268         }
 269         inexact = guardbit | stickybit;
 270
 271         /*
 272          * round result
 273          */
 274         if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
 275                 Dbl_clear_signexponent(opnd3p1);
 276                 switch (Rounding_mode()) {
 277                         case ROUNDPLUS:
 278                                 if (Dbl_iszero_sign(resultp1))
 279                                         Dbl_increment(opnd3p1,opnd3p2);
 280                                 break;
 281                         case ROUNDMINUS:
 282                                 if (Dbl_isone_sign(resultp1))
 283                                         Dbl_increment(opnd3p1,opnd3p2);
 284                                 break;
 285                         case ROUNDNEAREST:
 286                                 if (guardbit && (stickybit ||
 287                                     Dbl_isone_lowmantissap2(opnd3p2))) {
 288                                         Dbl_increment(opnd3p1,opnd3p2);
 289                                 }
 290                 }
 291                 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
 292         }
 293         Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
 294
 295         /*
 296          * Test for overflow
 297          */
 298         if (dest_exponent >= DBL_INFINITY_EXPONENT) {
 299                 /* trap if OVERFLOWTRAP enabled */
 300                 if (Is_overflowtrap_enabled()) {
 301                         /*
 302                          * Adjust bias of result
 303                          */
 304                         Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
 305                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 306                         if (inexact) {
 307                                 if (Is_inexacttrap_enabled())
 308                                         return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
 309                                 else
 310                                         Set_inexactflag();
 311                         }
 312                         return(OVERFLOWEXCEPTION);
 313                 }
 314                 Set_overflowflag();
 315                 /* set result to infinity or largest number */
 316                 Dbl_setoverflow(resultp1,resultp2);
 317                 inexact = true;
 318         }
 319         /*
 320          * Test for underflow
 321          */
 322         else if (dest_exponent <= 0) {
 323                 /* trap if UNDERFLOWTRAP enabled */
 324                 if (Is_underflowtrap_enabled()) {
 325                         /*
 326                          * Adjust bias of result
 327                          */
 328                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
 329                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 330                         if (inexact) {
 331                                 if (Is_inexacttrap_enabled())
 332                                         return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
 333                                 else
 334                                         Set_inexactflag();
 335                         }
 336                         return(UNDERFLOWEXCEPTION);
 337                 }
 338
 339                 /* Determine if should set underflow flag */
 340                 is_tiny = true;
 341                 if (dest_exponent == 0 && inexact) {
 342                         switch (Rounding_mode()) {
 343                         case ROUNDPLUS:
 344                                 if (Dbl_iszero_sign(resultp1)) {
 345                                         Dbl_increment(opnd3p1,opnd3p2);
 346                                         if (Dbl_isone_hiddenoverflow(opnd3p1))
 347                                                 is_tiny = false;
 348                                         Dbl_decrement(opnd3p1,opnd3p2);
 349                                 }
 350                                 break;
 351                         case ROUNDMINUS:
 352                                 if (Dbl_isone_sign(resultp1)) {
 353                                         Dbl_increment(opnd3p1,opnd3p2);
 354                                         if (Dbl_isone_hiddenoverflow(opnd3p1))
 355                                                 is_tiny = false;
 356                                         Dbl_decrement(opnd3p1,opnd3p2);
 357                                 }
 358                                 break;
 359                         case ROUNDNEAREST:
 360                                 if (guardbit && (stickybit ||
 361                                     Dbl_isone_lowmantissap2(opnd3p2))) {
 362                                         Dbl_increment(opnd3p1,opnd3p2);
 363                                         if (Dbl_isone_hiddenoverflow(opnd3p1))
 364                                                 is_tiny = false;
 365                                         Dbl_decrement(opnd3p1,opnd3p2);
 366                                 }
 367                                 break;
 368                         }
 369                 }
 370
 371                 /*
 372                  * denormalize result or set to signed zero
 373                  */
 374                 stickybit = inexact;
 375                 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
 376                  stickybit,inexact);
 377
 378                 /* return rounded number */
 379                 if (inexact) {
 380                         switch (Rounding_mode()) {
 381                         case ROUNDPLUS:
 382                                 if (Dbl_iszero_sign(resultp1)) {
 383                                         Dbl_increment(opnd3p1,opnd3p2);
 384                                 }
 385                                 break;
 386                         case ROUNDMINUS:
 387                                 if (Dbl_isone_sign(resultp1)) {
 388                                         Dbl_increment(opnd3p1,opnd3p2);
 389                                 }
 390                                 break;
 391                         case ROUNDNEAREST:
 392                                 if (guardbit && (stickybit ||
 393                                     Dbl_isone_lowmantissap2(opnd3p2))) {
 394                                         Dbl_increment(opnd3p1,opnd3p2);
 395                                 }
 396                                 break;
 397                         }
 398                         if (is_tiny)
 399                                 Set_underflowflag();
 400                 }
 401                 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
 402         }
 403         else Dbl_set_exponent(resultp1,dest_exponent);
 404         Dbl_copytoptr(resultp1,resultp2,dstptr);
 405
 406         /* check for inexact */
 407         if (inexact) {
 408                 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
 409                 else Set_inexactflag();
 410         }
 411         return(NOEXCEPTION);
 412 }