sys/arch/hppa/spmath/sfmpy.c

   1 /*      $NetBSD: sfmpy.c,v 1.3 2005/12/11 12:17:40 christos Exp $       */
   2
   3 /*      $OpenBSD: sfmpy.c,v 1.4 2001/03/29 03:58:19 mickey Exp $        */
   4
   5 /*
   6  * Copyright 1996 1995 by Open Software Foundation, Inc.
   7  *              All Rights Reserved
   8  *
   9  * Permission to use, copy, modify, and distribute this software and
  10  * its documentation for any purpose and without fee is hereby granted,
  11  * provided that the above copyright notice appears in all copies and
  12  * that both the copyright notice and this permission notice appear in
  13  * supporting documentation.
  14  *
  15  * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE
  16  * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  17  * FOR A PARTICULAR PURPOSE.
  18  *
  19  * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR
  20  * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
  21  * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT,
  22  * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
  23  * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  24  *
  25  */
  26 /*
  27  * pmk1.1
  28  */
  29 /*
  30  * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  31  *
  32  * To anyone who acknowledges that this file is provided "AS IS"
  33  * without any express or implied warranty:
  34  *     permission to use, copy, modify, and distribute this file
  35  * for any purpose is hereby granted without fee, provided that
  36  * the above copyright notice and this notice appears in all
  37  * copies, and that the name of Hewlett-Packard Company not be
  38  * used in advertising or publicity pertaining to distribution
  39  * of the software without specific, written prior permission.
  40  * Hewlett-Packard Company makes no representations about the
  41  * suitability of this software for any purpose.
  42  */
  43
  44 #include <sys/cdefs.h>
  45 __KERNEL_RCSID(0, "$NetBSD: sfmpy.c,v 1.3 2005/12/11 12:17:40 christos Exp $");
  46
  47 #include "../spmath/float.h"
  48 #include "../spmath/sgl_float.h"
  49
  50 /*
  51  *  Single Precision Floating-point Multiply
  52  */
  53 int
  54 sgl_fmpy(srcptr1,srcptr2,dstptr,status)
  55
  56 sgl_floating_point *srcptr1, *srcptr2, *dstptr;
  57 unsigned int *status;
  58 {
  59         register unsigned int opnd1, opnd2, opnd3, result;
  60         register int dest_exponent, count;
  61         register int inexact = false, guardbit = false, stickybit = false;
  62         int is_tiny;
  63
  64         opnd1 = *srcptr1;
  65         opnd2 = *srcptr2;
  66         /*
  67          * set sign bit of result
  68          */
  69         if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);
  70         else Sgl_setzero(result);
  71         /*
  72          * check first operand for NaN's or infinity
  73          */
  74         if (Sgl_isinfinity_exponent(opnd1)) {
  75                 if (Sgl_iszero_mantissa(opnd1)) {
  76                         if (Sgl_isnotnan(opnd2)) {
  77                                 if (Sgl_iszero_exponentmantissa(opnd2)) {
  78                                         /*
  79                                          * invalid since operands are infinity
  80                                          * and zero
  81                                          */
  82                                         if (Is_invalidtrap_enabled())
  83                                                 return(INVALIDEXCEPTION);
  84                                         Set_invalidflag();
  85                                         Sgl_makequietnan(result);
  86                                         *dstptr = result;
  87                                         return(NOEXCEPTION);
  88                                 }
  89                                 /*
  90                                  * return infinity
  91                                  */
  92                                 Sgl_setinfinity_exponentmantissa(result);
  93                                 *dstptr = result;
  94                                 return(NOEXCEPTION);
  95                         }
  96                 }
  97                 else {
  98                         /*
  99                          * is NaN; signaling or quiet?
 100                          */
 101                         if (Sgl_isone_signaling(opnd1)) {
 102                                 /* trap if INVALIDTRAP enabled */
 103                                 if (Is_invalidtrap_enabled())
 104                                         return(INVALIDEXCEPTION);
 105                                 /* make NaN quiet */
 106                                 Set_invalidflag();
 107                                 Sgl_set_quiet(opnd1);
 108                         }
 109                         /*
 110                          * is second operand a signaling NaN?
 111                          */
 112                         else if (Sgl_is_signalingnan(opnd2)) {
 113                                 /* trap if INVALIDTRAP enabled */
 114                                 if (Is_invalidtrap_enabled())
 115                                         return(INVALIDEXCEPTION);
 116                                 /* make NaN quiet */
 117                                 Set_invalidflag();
 118                                 Sgl_set_quiet(opnd2);
 119                                 *dstptr = opnd2;
 120                                 return(NOEXCEPTION);
 121                         }
 122                         /*
 123                          * return quiet NaN
 124                          */
 125                         *dstptr = opnd1;
 126                         return(NOEXCEPTION);
 127                 }
 128         }
 129         /*
 130          * check second operand for NaN's or infinity
 131          */
 132         if (Sgl_isinfinity_exponent(opnd2)) {
 133                 if (Sgl_iszero_mantissa(opnd2)) {
 134                         if (Sgl_iszero_exponentmantissa(opnd1)) {
 135                                 /* invalid since operands are zero & infinity */
 136                                 if (Is_invalidtrap_enabled())
 137                                         return(INVALIDEXCEPTION);
 138                                 Set_invalidflag();
 139                                 Sgl_makequietnan(opnd2);
 140                                 *dstptr = opnd2;
 141                                 return(NOEXCEPTION);
 142                         }
 143                         /*
 144                          * return infinity
 145                          */
 146                         Sgl_setinfinity_exponentmantissa(result);
 147                         *dstptr = result;
 148                         return(NOEXCEPTION);
 149                 }
 150                 /*
 151                  * is NaN; signaling or quiet?
 152                  */
 153                 if (Sgl_isone_signaling(opnd2)) {
 154                         /* trap if INVALIDTRAP enabled */
 155                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 156
 157                         /* make NaN quiet */
 158                         Set_invalidflag();
 159                         Sgl_set_quiet(opnd2);
 160                 }
 161                 /*
 162                  * return quiet NaN
 163                  */
 164                 *dstptr = opnd2;
 165                 return(NOEXCEPTION);
 166         }
 167         /*
 168          * Generate exponent
 169          */
 170         dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
 171
 172         /*
 173          * Generate mantissa
 174          */
 175         if (Sgl_isnotzero_exponent(opnd1)) {
 176                 /* set hidden bit */
 177                 Sgl_clear_signexponent_set_hidden(opnd1);
 178         }
 179         else {
 180                 /* check for zero */
 181                 if (Sgl_iszero_mantissa(opnd1)) {
 182                         Sgl_setzero_exponentmantissa(result);
 183                         *dstptr = result;
 184                         return(NOEXCEPTION);
 185                 }
 186                 /* is denormalized, adjust exponent */
 187                 Sgl_clear_signexponent(opnd1);
 188                 Sgl_leftshiftby1(opnd1);
 189                 Sgl_normalize(opnd1,dest_exponent);
 190         }
 191         /* opnd2 needs to have hidden bit set with msb in hidden bit */
 192         if (Sgl_isnotzero_exponent(opnd2)) {
 193                 Sgl_clear_signexponent_set_hidden(opnd2);
 194         }
 195         else {
 196                 /* check for zero */
 197                 if (Sgl_iszero_mantissa(opnd2)) {
 198                         Sgl_setzero_exponentmantissa(result);
 199                         *dstptr = result;
 200                         return(NOEXCEPTION);
 201                 }
 202                 /* is denormalized; want to normalize */
 203                 Sgl_clear_signexponent(opnd2);
 204                 Sgl_leftshiftby1(opnd2);
 205                 Sgl_normalize(opnd2,dest_exponent);
 206         }
 207
 208         /* Multiply two source mantissas together */
 209
 210         Sgl_leftshiftby4(opnd2);     /* make room for guard bits */
 211         Sgl_setzero(opnd3);
 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<SGL_P;count+=4) {
 217                 stickybit |= Slow4(opnd3);
 218                 Sgl_rightshiftby4(opnd3);
 219                 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3);
 220                 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2);
 221                 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1);
 222                 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2);
 223                 Sgl_rightshiftby4(opnd1);
 224         }
 225         /* make sure result is left-justified */
 226         if (Sgl_iszero_sign(opnd3)) {
 227                 Sgl_leftshiftby1(opnd3);
 228         }
 229         else {
 230                 /* result mantissa >= 2. */
 231                 dest_exponent++;
 232         }
 233         /* check for denormalized result */
 234         while (Sgl_iszero_sign(opnd3)) {
 235                 Sgl_leftshiftby1(opnd3);
 236                 dest_exponent--;
 237         }
 238         /*
 239          * check for guard, sticky and inexact bits
 240          */
 241         stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1);
 242         guardbit = Sbit24(opnd3);
 243         inexact = guardbit | stickybit;
 244
 245         /* re-align mantissa */
 246         Sgl_rightshiftby8(opnd3);
 247
 248         /*
 249          * round result
 250          */
 251         if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
 252                 Sgl_clear_signexponent(opnd3);
 253                 switch (Rounding_mode()) {
 254                         case ROUNDPLUS:
 255                                 if (Sgl_iszero_sign(result))
 256                                         Sgl_increment(opnd3);
 257                                 break;
 258                         case ROUNDMINUS:
 259                                 if (Sgl_isone_sign(result))
 260                                         Sgl_increment(opnd3);
 261                                 break;
 262                         case ROUNDNEAREST:
 263                                 if (guardbit &&
 264                                     (stickybit || Sgl_isone_lowmantissa(opnd3)))
 265                                         Sgl_increment(opnd3);
 266                                 break;
 267                 }
 268                 if (Sgl_isone_hidden(opnd3)) dest_exponent++;
 269         }
 270         Sgl_set_mantissa(result,opnd3);
 271
 272         /*
 273          * Test for overflow
 274          */
 275         if (dest_exponent >= SGL_INFINITY_EXPONENT) {
 276                 /* trap if OVERFLOWTRAP enabled */
 277                 if (Is_overflowtrap_enabled()) {
 278                         /*
 279                          * Adjust bias of result
 280                          */
 281                         Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
 282                         *dstptr = result;
 283                         if (inexact) {
 284                             if (Is_inexacttrap_enabled())
 285                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
 286                             else Set_inexactflag();
 287                         }
 288                         return(OVERFLOWEXCEPTION);
 289                 }
 290                 inexact = true;
 291                 Set_overflowflag();
 292                 /* set result to infinity or largest number */
 293                 Sgl_setoverflow(result);
 294         }
 295         /*
 296          * Test for underflow
 297          */
 298         else if (dest_exponent <= 0) {
 299                 /* trap if UNDERFLOWTRAP enabled */
 300                 if (Is_underflowtrap_enabled()) {
 301                         /*
 302                          * Adjust bias of result
 303                          */
 304                         Sgl_setwrapped_exponent(result,dest_exponent,unfl);
 305                         *dstptr = result;
 306                         if (inexact) {
 307                             if (Is_inexacttrap_enabled())
 308                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
 309                             else Set_inexactflag();
 310                         }
 311                         return(UNDERFLOWEXCEPTION);
 312                 }
 313
 314                 /* Determine if should set underflow flag */
 315                 is_tiny = true;
 316                 if (dest_exponent == 0 && inexact) {
 317                         switch (Rounding_mode()) {
 318                         case ROUNDPLUS:
 319                                 if (Sgl_iszero_sign(result)) {
 320                                         Sgl_increment(opnd3);
 321                                         if (Sgl_isone_hiddenoverflow(opnd3))
 322                                                 is_tiny = false;
 323                                         Sgl_decrement(opnd3);
 324                                 }
 325                                 break;
 326                         case ROUNDMINUS:
 327                                 if (Sgl_isone_sign(result)) {
 328                                         Sgl_increment(opnd3);
 329                                         if (Sgl_isone_hiddenoverflow(opnd3))
 330                                                 is_tiny = false;
 331                                         Sgl_decrement(opnd3);
 332                                 }
 333                                 break;
 334                         case ROUNDNEAREST:
 335                                 if (guardbit && (stickybit ||
 336                                     Sgl_isone_lowmantissa(opnd3))) {
 337                                         Sgl_increment(opnd3);
 338                                         if (Sgl_isone_hiddenoverflow(opnd3))
 339                                                 is_tiny = false;
 340                                         Sgl_decrement(opnd3);
 341                                 }
 342                                 break;
 343                         }
 344                 }
 345
 346                 /*
 347                  * denormalize result or set to signed zero
 348                  */
 349                 stickybit = inexact;
 350                 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
 351
 352                 /* return zero or smallest number */
 353                 if (inexact) {
 354                         switch (Rounding_mode()) {
 355                         case ROUNDPLUS:
 356                                 if (Sgl_iszero_sign(result))
 357                                         Sgl_increment(opnd3);
 358                                 break;
 359                         case ROUNDMINUS:
 360                                 if (Sgl_isone_sign(result))
 361                                         Sgl_increment(opnd3);
 362                                 break;
 363                         case ROUNDNEAREST:
 364                                 if (guardbit && (stickybit ||
 365                                     Sgl_isone_lowmantissa(opnd3)))
 366                                         Sgl_increment(opnd3);
 367                                 break;
 368                         }
 369                 if (is_tiny) Set_underflowflag();
 370                 }
 371                 Sgl_set_exponentmantissa(result,opnd3);
 372         }
 373         else Sgl_set_exponent(result,dest_exponent);
 374         *dstptr = result;
 375
 376         /* check for inexact */
 377         if (inexact) {
 378                 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
 379                 else Set_inexactflag();
 380         }
 381         return(NOEXCEPTION);
 382 }