sysdeps/ieee754/dbl-64/s_nearbyint.c

   1 /* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>.  */
   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 #if defined(LIBM_SCCS) && !defined(lint)
  14 static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
  15 #endif
  16
  17 /*
  18  * rint(x)
  19  * Return x rounded to integral value according to the prevailing
  20  * rounding mode.
  21  * Method:
  22  *      Using floating addition.
  23  * Exception:
  24  *      Inexact flag raised if x not equal to rint(x).
  25  */
  26
  27 #include <fenv.h>
  28 #include <math.h>
  29 #include <math_private.h>
  30
  31 static const double
  32   TWO52[2] = {
  33   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  34  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  35 };
  36
  37 double
  38 __nearbyint (double x)
  39 {
  40   fenv_t env;
  41   int32_t i0, j0, sx;
  42   double w, t;
  43   GET_HIGH_WORD (i0, x);
  44   sx = (i0 >> 31) & 1;
  45   j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
  46   if (j0 < 52)
  47     {
  48       if (j0 < 0)
  49         {
  50           libc_feholdexcept (&env);
  51           w = TWO52[sx] + x;
  52           t = w - TWO52[sx];
  53           math_force_eval (t);
  54           libc_fesetenv (&env);
  55           GET_HIGH_WORD (i0, t);
  56           SET_HIGH_WORD (t, (i0 & 0x7fffffff) | (sx << 31));
  57           return t;
  58         }
  59     }
  60   else
  61     {
  62       if (j0 == 0x400)
  63         return x + x;                   /* inf or NaN */
  64       else
  65         return x;                       /* x is integral */
  66     }
  67   libc_feholdexcept (&env);
  68   w = TWO52[sx] + x;
  69   t = w - TWO52[sx];
  70   math_force_eval (t);
  71   libc_fesetenv (&env);
  72   return t;
  73 }
  74 weak_alias (__nearbyint, nearbyint)
  75 #ifdef NO_LONG_DOUBLE
  76 strong_alias (__nearbyint, __nearbyintl)
  77 weak_alias (__nearbyint, nearbyintl)
  78 #endif