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