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