newlib/libm/common/s_round.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11 /*
  12 FUNCTION
  13 <<round>>, <<roundf>>--round to integer, to nearest
  14 INDEX
  15         round
  16 INDEX
  17         roundf
  18
  19 ANSI_SYNOPSIS
  20         #include <math.h>
  21         double round(double <[x]>);
  22         float roundf(float <[x]>);
  23
  24 DESCRIPTION
  25         The <<round>> functions round their argument to the nearest integer
  26         value in floating-point format, rounding halfway cases away from zero,
  27         regardless of the current rounding direction.  (While the "inexact"
  28         floating-point exception behavior is unspecified by the C standard, the
  29         <<round>> functions are written so that "inexact" is not raised if the
  30         result does not equal the argument, which behavior is as recommended by
  31         IEEE 754 for its related functions.)
  32
  33 RETURNS
  34 <[x]> rounded to an integral value.
  35
  36 PORTABILITY
  37 ANSI C, POSIX
  38
  39 SEEALSO
  40 <<nearbyint>>, <<rint>>
  41
  42 */
  43
  44 #include "fdlibm.h"
  45
  46 #ifndef _DOUBLE_IS_32BITS
  47
  48 #ifdef __STDC__
  49         double round(double x)
  50 #else
  51         double round(x)
  52         double x;
  53 #endif
  54 {
  55   /* Most significant word, least significant word. */
  56   __int32_t msw, exponent_less_1023;
  57   __uint32_t lsw;
  58
  59   EXTRACT_WORDS(msw, lsw, x);
  60
  61   /* Extract exponent field. */
  62   exponent_less_1023 = ((msw & 0x7ff00000) >> 20) - 1023;
  63
  64   if (exponent_less_1023 < 20)
  65     {
  66       if (exponent_less_1023 < 0)
  67         {
  68           msw &= 0x80000000;
  69           if (exponent_less_1023 == -1)
  70             /* Result is +1.0 or -1.0. */
  71             msw |= (1023 << 20);
  72           lsw = 0;
  73         }
  74       else
  75         {
  76           __uint32_t exponent_mask = 0x000fffff >> exponent_less_1023;
  77           if ((msw & exponent_mask) == 0 && lsw == 0)
  78             /* x in an integral value. */
  79             return x;
  80
  81           msw += 0x00080000 >> exponent_less_1023;
  82           msw &= ~exponent_mask;
  83           lsw = 0;
  84         }
  85     }
  86   else if (exponent_less_1023 > 51)
  87     {
  88       if (exponent_less_1023 == 1024)
  89         /* x is NaN or infinite. */
  90         return x + x;
  91       else
  92         return x;
  93     }
  94   else
  95     {
  96       __uint32_t exponent_mask = 0xffffffff >> (exponent_less_1023 - 20);
  97       __uint32_t tmp;
  98
  99       if ((lsw & exponent_mask) == 0)
 100         /* x is an integral value. */
 101         return x;
 102
 103       tmp = lsw + (1 << (51 - exponent_less_1023));
 104       if (tmp < lsw)
 105         msw += 1;
 106       lsw = tmp;
 107
 108       lsw &= ~exponent_mask;
 109     }
 110   INSERT_WORDS(x, msw, lsw);
 111
 112   return x;
 113 }
 114
 115 #endif /* _DOUBLE_IS_32BITS */