LXR linux/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
  49int
  50dbl_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}
 401