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
  49 int
  50 dbl_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 }