LXR linux/arch/parisc/math-emu/dfmpy.c

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