sysdeps/ieee754/dbl-64/s_rint.c

   1 /* @(#)s_rint.c 5.1 93/09/24 */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12
  13 /*
  14  * rint(x)
  15  * Return x rounded to integral value according to the prevailing
  16  * rounding mode.
  17  * Method:
  18  *      Using floating addition.
  19  * Exception:
  20  *      Inexact flag raised if x not equal to rint(x).
  21  */
  22
  23 #include <math.h>
  24 #include <math_private.h>
  25
  26 static const double
  27   TWO52[2] = {
  28   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  29  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  30 };
  31
  32 double
  33 __rint (double x)
  34 {
  35   int32_t i0, j0, sx;
  36   double w, t;
  37   GET_HIGH_WORD (i0, x);
  38   sx = (i0 >> 31) & 1;
  39   j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
  40   if (j0 < 52)
  41     {
  42       if (j0 < 0)
  43         {
  44           w = TWO52[sx] + x;
  45           t = w - TWO52[sx];
  46           GET_HIGH_WORD (i0, t);
  47           SET_HIGH_WORD (t, (i0 & 0x7fffffff) | (sx << 31));
  48           return t;
  49         }
  50     }
  51   else
  52     {
  53       if (j0 == 0x400)
  54         return x + x;                   /* inf or NaN */
  55       else
  56         return x;                       /* x is integral */
  57     }
  58   w = TWO52[sx] + x;
  59   return w - TWO52[sx];
  60 }
  61 #ifndef __rint
  62 weak_alias (__rint, rint)
  63 # ifdef NO_LONG_DOUBLE
  64 strong_alias (__rint, __rintl)
  65 weak_alias (__rint, rintl)
  66 # endif
  67 #endif