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