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) <[email protected]>
   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
  36 int
  37 dbl_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 }