arch/parisc/math-emu/sfmpy.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/sfmpy.c               $Revision: 1.1 $
  26  *
  27  *  Purpose:
  28  *      Single Precision Floating-point Multiply
  29  *
  30  *  External Interfaces:
  31  *      sgl_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 "sgl_float.h"
  44
  45 /*
  46  *  Single Precision Floating-point Multiply
  47  */
  48
  49 int
  50 sgl_fmpy(
  51     sgl_floating_point *srcptr1,
  52     sgl_floating_point *srcptr2,
  53     sgl_floating_point *dstptr,
  54     unsigned int *status)
  55 {
  56         register unsigned int opnd1, opnd2, opnd3, result;
  57         register int dest_exponent, count;
  58         register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
  59         boolean is_tiny;
  60
  61         opnd1 = *srcptr1;
  62         opnd2 = *srcptr2;
  63         /*
  64          * set sign bit of result
  65          */
  66         if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
  67         else Sgl_setzero(result);
  68         /*
  69          * check first operand for NaN's or infinity
  70          */
  71         if (Sgl_isinfinity_exponent(opnd1)) {
  72                 if (Sgl_iszero_mantissa(opnd1)) {
  73                         if (Sgl_isnotnan(opnd2)) {
  74                                 if (Sgl_iszero_exponentmantissa(opnd2)) {
  75                                         /*
  76                                          * invalid since operands are infinity
  77                                          * and zero
  78                                          */
  79                                         if (Is_invalidtrap_enabled())
  80                                                 return(INVALIDEXCEPTION);
  81                                         Set_invalidflag();
  82                                         Sgl_makequietnan(result);
  83                                         *dstptr = result;
  84                                         return(NOEXCEPTION);
  85                                 }
  86                                 /*
  87                                  * return infinity
  88                                  */
  89                                 Sgl_setinfinity_exponentmantissa(result);
  90                                 *dstptr = result;
  91                                 return(NOEXCEPTION);
  92                         }
  93                 }
  94                 else {
  95                         /*
  96                          * is NaN; signaling or quiet?
  97                          */
  98                         if (Sgl_isone_signaling(opnd1)) {
  99                                 /* trap if INVALIDTRAP enabled */
 100                                 if (Is_invalidtrap_enabled())
 101                                         return(INVALIDEXCEPTION);
 102                                 /* make NaN quiet */
 103                                 Set_invalidflag();
 104                                 Sgl_set_quiet(opnd1);
 105                         }
 106                         /*
 107                          * is second operand a signaling NaN?
 108                          */
 109                         else if (Sgl_is_signalingnan(opnd2)) {
 110                                 /* trap if INVALIDTRAP enabled */
 111                                 if (Is_invalidtrap_enabled())
 112                                         return(INVALIDEXCEPTION);
 113                                 /* make NaN quiet */
 114                                 Set_invalidflag();
 115                                 Sgl_set_quiet(opnd2);
 116                                 *dstptr = opnd2;
 117                                 return(NOEXCEPTION);
 118                         }
 119                         /*
 120                          * return quiet NaN
 121                          */
 122                         *dstptr = opnd1;
 123                         return(NOEXCEPTION);
 124                 }
 125         }
 126         /*
 127          * check second operand for NaN's or infinity
 128          */
 129         if (Sgl_isinfinity_exponent(opnd2)) {
 130                 if (Sgl_iszero_mantissa(opnd2)) {
 131                         if (Sgl_iszero_exponentmantissa(opnd1)) {
 132                                 /* invalid since operands are zero & infinity */
 133                                 if (Is_invalidtrap_enabled())
 134                                         return(INVALIDEXCEPTION);
 135                                 Set_invalidflag();
 136                                 Sgl_makequietnan(opnd2);
 137                                 *dstptr = opnd2;
 138                                 return(NOEXCEPTION);
 139                         }
 140                         /*
 141                          * return infinity
 142                          */
 143                         Sgl_setinfinity_exponentmantissa(result);
 144                         *dstptr = result;
 145                         return(NOEXCEPTION);
 146                 }
 147                 /*
 148                  * is NaN; signaling or quiet?
 149                  */
 150                 if (Sgl_isone_signaling(opnd2)) {
 151                         /* trap if INVALIDTRAP enabled */
 152                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 153
 154                         /* make NaN quiet */
 155                         Set_invalidflag();
 156                         Sgl_set_quiet(opnd2);
 157                 }
 158                 /*
 159                  * return quiet NaN
 160                  */
 161                 *dstptr = opnd2;
 162                 return(NOEXCEPTION);
 163         }
 164         /*
 165          * Generate exponent
 166          */
 167         dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
 168
 169         /*
 170          * Generate mantissa
 171          */
 172         if (Sgl_isnotzero_exponent(opnd1)) {
 173                 /* set hidden bit */
 174                 Sgl_clear_signexponent_set_hidden(opnd1);
 175         }
 176         else {
 177                 /* check for zero */
 178                 if (Sgl_iszero_mantissa(opnd1)) {
 179                         Sgl_setzero_exponentmantissa(result);
 180                         *dstptr = result;
 181                         return(NOEXCEPTION);
 182                 }
 183                 /* is denormalized, adjust exponent */
 184                 Sgl_clear_signexponent(opnd1);
 185                 Sgl_leftshiftby1(opnd1);
 186                 Sgl_normalize(opnd1,dest_exponent);
 187         }
 188         /* opnd2 needs to have hidden bit set with msb in hidden bit */
 189         if (Sgl_isnotzero_exponent(opnd2)) {
 190                 Sgl_clear_signexponent_set_hidden(opnd2);
 191         }
 192         else {
 193                 /* check for zero */
 194                 if (Sgl_iszero_mantissa(opnd2)) {
 195                         Sgl_setzero_exponentmantissa(result);
 196                         *dstptr = result;
 197                         return(NOEXCEPTION);
 198                 }
 199                 /* is denormalized; want to normalize */
 200                 Sgl_clear_signexponent(opnd2);
 201                 Sgl_leftshiftby1(opnd2);
 202                 Sgl_normalize(opnd2,dest_exponent);
 203         }
 204
 205         /* Multiply two source mantissas together */
 206
 207         Sgl_leftshiftby4(opnd2);     /* make room for guard bits */
 208         Sgl_setzero(opnd3);
 209         /*
 210          * Four bits at a time are inspected in each loop, and a
 211          * simple shift and add multiply algorithm is used.
 212          */
 213         for (count=1;count<SGL_P;count+=4) {
 214                 stickybit |= Slow4(opnd3);
 215                 Sgl_rightshiftby4(opnd3);
 216                 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3);
 217                 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2);
 218                 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1);
 219                 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2);
 220                 Sgl_rightshiftby4(opnd1);
 221         }
 222         /* make sure result is left-justified */
 223         if (Sgl_iszero_sign(opnd3)) {
 224                 Sgl_leftshiftby1(opnd3);
 225         }
 226         else {
 227                 /* result mantissa >= 2. */
 228                 dest_exponent++;
 229         }
 230         /* check for denormalized result */
 231         while (Sgl_iszero_sign(opnd3)) {
 232                 Sgl_leftshiftby1(opnd3);
 233                 dest_exponent--;
 234         }
 235         /*
 236          * check for guard, sticky and inexact bits
 237          */
 238         stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1);
 239         guardbit = Sbit24(opnd3);
 240         inexact = guardbit | stickybit;
 241
 242         /* re-align mantissa */
 243         Sgl_rightshiftby8(opnd3);
 244
 245         /*
 246          * round result
 247          */
 248         if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
 249                 Sgl_clear_signexponent(opnd3);
 250                 switch (Rounding_mode()) {
 251                         case ROUNDPLUS:
 252                                 if (Sgl_iszero_sign(result))
 253                                         Sgl_increment(opnd3);
 254                                 break;
 255                         case ROUNDMINUS:
 256                                 if (Sgl_isone_sign(result))
 257                                         Sgl_increment(opnd3);
 258                                 break;
 259                         case ROUNDNEAREST:
 260                                 if (guardbit) {
 261                                 if (stickybit || Sgl_isone_lowmantissa(opnd3))
 262                                 Sgl_increment(opnd3);
 263                                 }
 264                 }
 265                 if (Sgl_isone_hidden(opnd3)) dest_exponent++;
 266         }
 267         Sgl_set_mantissa(result,opnd3);
 268
 269         /*
 270          * Test for overflow
 271          */
 272         if (dest_exponent >= SGL_INFINITY_EXPONENT) {
 273                 /* trap if OVERFLOWTRAP enabled */
 274                 if (Is_overflowtrap_enabled()) {
 275                         /*
 276                          * Adjust bias of result
 277                          */
 278                         Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
 279                         *dstptr = result;
 280                         if (inexact)
 281                             if (Is_inexacttrap_enabled())
 282                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
 283                             else Set_inexactflag();
 284                         return(OVERFLOWEXCEPTION);
 285                 }
 286                 inexact = TRUE;
 287                 Set_overflowflag();
 288                 /* set result to infinity or largest number */
 289                 Sgl_setoverflow(result);
 290         }
 291         /*
 292          * Test for underflow
 293          */
 294         else if (dest_exponent <= 0) {
 295                 /* trap if UNDERFLOWTRAP enabled */
 296                 if (Is_underflowtrap_enabled()) {
 297                         /*
 298                          * Adjust bias of result
 299                          */
 300                         Sgl_setwrapped_exponent(result,dest_exponent,unfl);
 301                         *dstptr = result;
 302                         if (inexact)
 303                             if (Is_inexacttrap_enabled())
 304                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
 305                             else Set_inexactflag();
 306                         return(UNDERFLOWEXCEPTION);
 307                 }
 308
 309                 /* Determine if should set underflow flag */
 310                 is_tiny = TRUE;
 311                 if (dest_exponent == 0 && inexact) {
 312                         switch (Rounding_mode()) {
 313                         case ROUNDPLUS:
 314                                 if (Sgl_iszero_sign(result)) {
 315                                         Sgl_increment(opnd3);
 316                                         if (Sgl_isone_hiddenoverflow(opnd3))
 317                                             is_tiny = FALSE;
 318                                         Sgl_decrement(opnd3);
 319                                 }
 320                                 break;
 321                         case ROUNDMINUS:
 322                                 if (Sgl_isone_sign(result)) {
 323                                         Sgl_increment(opnd3);
 324                                         if (Sgl_isone_hiddenoverflow(opnd3))
 325                                             is_tiny = FALSE;
 326                                         Sgl_decrement(opnd3);
 327                                 }
 328                                 break;
 329                         case ROUNDNEAREST:
 330                                 if (guardbit && (stickybit ||
 331                                     Sgl_isone_lowmantissa(opnd3))) {
 332                                         Sgl_increment(opnd3);
 333                                         if (Sgl_isone_hiddenoverflow(opnd3))
 334                                             is_tiny = FALSE;
 335                                         Sgl_decrement(opnd3);
 336                                 }
 337                                 break;
 338                         }
 339                 }
 340
 341                 /*
 342                  * denormalize result or set to signed zero
 343                  */
 344                 stickybit = inexact;
 345                 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
 346
 347                 /* return zero or smallest number */
 348                 if (inexact) {
 349                         switch (Rounding_mode()) {
 350                         case ROUNDPLUS:
 351                                 if (Sgl_iszero_sign(result)) {
 352                                         Sgl_increment(opnd3);
 353                                 }
 354                                 break;
 355                         case ROUNDMINUS:
 356                                 if (Sgl_isone_sign(result)) {
 357                                         Sgl_increment(opnd3);
 358                                 }
 359                                 break;
 360                         case ROUNDNEAREST:
 361                                 if (guardbit && (stickybit ||
 362                                     Sgl_isone_lowmantissa(opnd3))) {
 363                                         Sgl_increment(opnd3);
 364                                 }
 365                                 break;
 366                         }
 367                 if (is_tiny) Set_underflowflag();
 368                 }
 369                 Sgl_set_exponentmantissa(result,opnd3);
 370         }
 371         else Sgl_set_exponent(result,dest_exponent);
 372         *dstptr = result;
 373
 374         /* check for inexact */
 375         if (inexact) {
 376                 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
 377                 else Set_inexactflag();
 378         }
 379         return(NOEXCEPTION);
 380 }