sys/arch/hppa/spmath/sfdiv.c

   1 /*      $NetBSD: sfdiv.c,v 1.3 2005/12/11 12:17:40 christos Exp $       */
   2
   3 /*      $OpenBSD: sfdiv.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: sfdiv.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 Divide
  52  */
  53 int
  54 sgl_fdiv(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_isinfinity(opnd2)) {
  78                                         /*
  79                                          * invalid since both operands
  80                                          * are infinity
  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                         /*
 135                          * return zero
 136                          */
 137                         Sgl_setzero_exponentmantissa(result);
 138                         *dstptr = result;
 139                         return(NOEXCEPTION);
 140                 }
 141                 /*
 142                  * is NaN; signaling or quiet?
 143                  */
 144                 if (Sgl_isone_signaling(opnd2)) {
 145                         /* trap if INVALIDTRAP enabled */
 146                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 147                         /* make NaN quiet */
 148                         Set_invalidflag();
 149                         Sgl_set_quiet(opnd2);
 150                 }
 151                 /*
 152                  * return quiet NaN
 153                  */
 154                 *dstptr = opnd2;
 155                 return(NOEXCEPTION);
 156         }
 157         /*
 158          * check for division by zero
 159          */
 160         if (Sgl_iszero_exponentmantissa(opnd2)) {
 161                 if (Sgl_iszero_exponentmantissa(opnd1)) {
 162                         /* invalid since both operands are zero */
 163                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 164                         Set_invalidflag();
 165                         Sgl_makequietnan(result);
 166                         *dstptr = result;
 167                         return(NOEXCEPTION);
 168                 }
 169                 if (Is_divisionbyzerotrap_enabled())
 170                         return(DIVISIONBYZEROEXCEPTION);
 171                 Set_divisionbyzeroflag();
 172                 Sgl_setinfinity_exponentmantissa(result);
 173                 *dstptr = result;
 174                 return(NOEXCEPTION);
 175         }
 176         /*
 177          * Generate exponent
 178          */
 179         dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
 180
 181         /*
 182          * Generate mantissa
 183          */
 184         if (Sgl_isnotzero_exponent(opnd1)) {
 185                 /* set hidden bit */
 186                 Sgl_clear_signexponent_set_hidden(opnd1);
 187         }
 188         else {
 189                 /* check for zero */
 190                 if (Sgl_iszero_mantissa(opnd1)) {
 191                         Sgl_setzero_exponentmantissa(result);
 192                         *dstptr = result;
 193                         return(NOEXCEPTION);
 194                 }
 195                 /* is denormalized; want to normalize */
 196                 Sgl_clear_signexponent(opnd1);
 197                 Sgl_leftshiftby1(opnd1);
 198                 Sgl_normalize(opnd1,dest_exponent);
 199         }
 200         /* opnd2 needs to have hidden bit set with msb in hidden bit */
 201         if (Sgl_isnotzero_exponent(opnd2)) {
 202                 Sgl_clear_signexponent_set_hidden(opnd2);
 203         }
 204         else {
 205                 /* is denormalized; want to normalize */
 206                 Sgl_clear_signexponent(opnd2);
 207                 Sgl_leftshiftby1(opnd2);
 208                 while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
 209                         Sgl_leftshiftby8(opnd2);
 210                         dest_exponent += 8;
 211                 }
 212                 if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
 213                         Sgl_leftshiftby4(opnd2);
 214                         dest_exponent += 4;
 215                 }
 216                 while(Sgl_iszero_hidden(opnd2)) {
 217                         Sgl_leftshiftby1(opnd2);
 218                         dest_exponent += 1;
 219                 }
 220         }
 221
 222         /* Divide the source mantissas */
 223
 224         /*
 225          * A non_restoring divide algorithm is used.
 226          */
 227         Sgl_subtract(opnd1,opnd2,opnd1);
 228         Sgl_setzero(opnd3);
 229         for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
 230                 Sgl_leftshiftby1(opnd1);
 231                 Sgl_leftshiftby1(opnd3);
 232                 if (Sgl_iszero_sign(opnd1)) {
 233                         Sgl_setone_lowmantissa(opnd3);
 234                         Sgl_subtract(opnd1,opnd2,opnd1);
 235                 }
 236                 else Sgl_addition(opnd1,opnd2,opnd1);
 237         }
 238         if (count <= SGL_P) {
 239                 Sgl_leftshiftby1(opnd3);
 240                 Sgl_setone_lowmantissa(opnd3);
 241                 Sgl_leftshift(opnd3,SGL_P-count);
 242                 if (Sgl_iszero_hidden(opnd3)) {
 243                         Sgl_leftshiftby1(opnd3);
 244                         dest_exponent--;
 245                 }
 246         }
 247         else {
 248                 if (Sgl_iszero_hidden(opnd3)) {
 249                         /* need to get one more bit of result */
 250                         Sgl_leftshiftby1(opnd1);
 251                         Sgl_leftshiftby1(opnd3);
 252                         if (Sgl_iszero_sign(opnd1)) {
 253                                 Sgl_setone_lowmantissa(opnd3);
 254                                 Sgl_subtract(opnd1,opnd2,opnd1);
 255                         }
 256                         else Sgl_addition(opnd1,opnd2,opnd1);
 257                         dest_exponent--;
 258                 }
 259                 if (Sgl_iszero_sign(opnd1)) guardbit = true;
 260                 stickybit = Sgl_all(opnd1);
 261         }
 262         inexact = guardbit | stickybit;
 263
 264         /*
 265          * round result
 266          */
 267         if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
 268                 Sgl_clear_signexponent(opnd3);
 269                 switch (Rounding_mode()) {
 270                         case ROUNDPLUS:
 271                                 if (Sgl_iszero_sign(result))
 272                                         Sgl_increment_mantissa(opnd3);
 273                                 break;
 274                         case ROUNDMINUS:
 275                                 if (Sgl_isone_sign(result))
 276                                         Sgl_increment_mantissa(opnd3);
 277                                 break;
 278                         case ROUNDNEAREST:
 279                                 if (guardbit &&
 280                                     (stickybit || Sgl_isone_lowmantissa(opnd3)))
 281                                         Sgl_increment_mantissa(opnd3);
 282                 }
 283                 if (Sgl_isone_hidden(opnd3)) dest_exponent++;
 284         }
 285         Sgl_set_mantissa(result,opnd3);
 286
 287         /*
 288          * Test for overflow
 289          */
 290         if (dest_exponent >= SGL_INFINITY_EXPONENT) {
 291                 /* trap if OVERFLOWTRAP enabled */
 292                 if (Is_overflowtrap_enabled()) {
 293                         /*
 294                          * Adjust bias of result
 295                          */
 296                         Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
 297                         *dstptr = result;
 298                         if (inexact) {
 299                             if (Is_inexacttrap_enabled())
 300                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
 301                             else Set_inexactflag();
 302                         }
 303                         return(OVERFLOWEXCEPTION);
 304                 }
 305                 Set_overflowflag();
 306                 /* set result to infinity or largest number */
 307                 Sgl_setoverflow(result);
 308                 inexact = true;
 309         }
 310         /*
 311          * Test for underflow
 312          */
 313         else if (dest_exponent <= 0) {
 314                 /* trap if UNDERFLOWTRAP enabled */
 315                 if (Is_underflowtrap_enabled()) {
 316                         /*
 317                          * Adjust bias of result
 318                          */
 319                         Sgl_setwrapped_exponent(result,dest_exponent,unfl);
 320                         *dstptr = result;
 321                         if (inexact) {
 322                             if (Is_inexacttrap_enabled())
 323                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
 324                             else Set_inexactflag();
 325                         }
 326                         return(UNDERFLOWEXCEPTION);
 327                 }
 328
 329                 /* Determine if should set underflow flag */
 330                 is_tiny = true;
 331                 if (dest_exponent == 0 && inexact) {
 332                         switch (Rounding_mode()) {
 333                         case ROUNDPLUS:
 334                                 if (Sgl_iszero_sign(result)) {
 335                                         Sgl_increment(opnd3);
 336                                         if (Sgl_isone_hiddenoverflow(opnd3))
 337                                                 is_tiny = false;
 338                                         Sgl_decrement(opnd3);
 339                                 }
 340                                 break;
 341                         case ROUNDMINUS:
 342                                 if (Sgl_isone_sign(result)) {
 343                                         Sgl_increment(opnd3);
 344                                         if (Sgl_isone_hiddenoverflow(opnd3))
 345                                                 is_tiny = false;
 346                                         Sgl_decrement(opnd3);
 347                                 }
 348                                 break;
 349                         case ROUNDNEAREST:
 350                                 if (guardbit && (stickybit ||
 351                                     Sgl_isone_lowmantissa(opnd3))) {
 352                                         Sgl_increment(opnd3);
 353                                         if (Sgl_isone_hiddenoverflow(opnd3))
 354                                                 is_tiny = false;
 355                                         Sgl_decrement(opnd3);
 356                                 }
 357                                 break;
 358                         }
 359                 }
 360
 361                 /*
 362                  * denormalize result or set to signed zero
 363                  */
 364                 stickybit = inexact;
 365                 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
 366
 367                 /* return rounded number */
 368                 if (inexact) {
 369                         switch (Rounding_mode()) {
 370                         case ROUNDPLUS:
 371                                 if (Sgl_iszero_sign(result)) {
 372                                         Sgl_increment(opnd3);
 373                                 }
 374                                 break;
 375                         case ROUNDMINUS:
 376                                 if (Sgl_isone_sign(result))  {
 377                                         Sgl_increment(opnd3);
 378                                 }
 379                                 break;
 380                         case ROUNDNEAREST:
 381                                 if (guardbit && (stickybit ||
 382                                     Sgl_isone_lowmantissa(opnd3))) {
 383                                         Sgl_increment(opnd3);
 384                                 }
 385                                 break;
 386                         }
 387                         if (is_tiny) Set_underflowflag();
 388                 }
 389                 Sgl_set_exponentmantissa(result,opnd3);
 390         }
 391         else Sgl_set_exponent(result,dest_exponent);
 392         *dstptr = result;
 393         /* check for inexact */
 394         if (inexact) {
 395                 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
 396                 else  Set_inexactflag();
 397         }
 398         return(NOEXCEPTION);
 399 }