newlib/libm/common/s_llrint.c

   1 /* lrint adapted to be llrint for Newlib, 2009 by Craig Howland.  */
   2 /* @(#)s_lrint.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 /*
  15  * llrint(x)
  16  * Return x rounded to integral value according to the prevailing
  17  * rounding mode.
  18  * Method:
  19  *      Using floating addition.
  20  * Exception:
  21  *      Inexact flag raised if x not equal to llrint(x).
  22  */
  23
  24 #include "fdlibm.h"
  25
  26 #ifndef _DOUBLE_IS_32BITS
  27
  28 #ifdef __STDC__
  29 static const double
  30 #else
  31 static double
  32 #endif
  33
  34 /* Adding a double, x, to 2^52 will cause the result to be rounded based on
  35    the fractional part of x, according to the implementation's current rounding
  36    mode.  2^52 is the smallest double that can be represented using all 52 significant
  37    digits. */
  38 TWO52[2]={
  39   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  40  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  41 };
  42
  43 long long int
  44 #ifdef __STDC__
  45         llrint(double x)
  46 #else
  47         llrint(x)
  48         double x;
  49 #endif
  50 {
  51   __int32_t i0,j0,sx;
  52   __uint32_t i1;
  53   double t;
  54   volatile double w;
  55   long long int result;
  56
  57   EXTRACT_WORDS(i0,i1,x);
  58
  59   /* Extract sign bit. */
  60   sx = (i0>>31)&1;
  61
  62   /* Extract exponent field. */
  63   j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
  64   /* j0 in [-1023,1024] */
  65
  66   if(j0 < 20)
  67     {
  68       /* j0 in [-1023,19] */
  69       w = TWO52[sx] + x;
  70       t = w - TWO52[sx];
  71       GET_HIGH_WORD(i0, t);
  72       /* Detect the all-zeros representation of plus and
  73          minus zero, which fails the calculation below. */
  74       if ((i0 & ~((__int32_t)1 << 31)) == 0)
  75           return 0;
  76       j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
  77       i0 &= 0x000fffff;
  78       i0 |= 0x00100000;
  79       result = (j0 < 0 ? 0 : i0 >> (20 - j0));
  80     }
  81   else if (j0 < (int)(8 * sizeof (long long int)) - 1)
  82     {
  83       /* 64bit return: j0 in [20,62] */
  84       if (j0 >= 52)
  85         /* 64bit return: j0 in [52,62] */
  86         /* 64bit return: left shift amt in [32,42] */
  87         result = ((long long int) ((i0 & 0x000fffff) | 0x00100000) << (j0 - 20)) |
  88                 /* 64bit return: right shift amt in [0,10] */
  89                    ((long long int) i1 << (j0 - 52));
  90       else
  91         {
  92           /* 64bit return: j0 in [20,51] */
  93           w = TWO52[sx] + x;
  94           t = w - TWO52[sx];
  95           EXTRACT_WORDS (i0, i1, t);
  96           j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
  97           i0 &= 0x000fffff;
  98           i0 |= 0x00100000;
  99           /* After round:
 100            * 64bit return: j0 in [20,52] */
 101           /* 64bit return: left shift amt in [0,32] */
 102           /* ***64bit return: right shift amt in [32,0] */
 103           result = ((long long int) i0 << (j0 - 20))
 104                         | SAFE_RIGHT_SHIFT (i1, (52 - j0));
 105         }
 106     }
 107   else
 108     {
 109       return (long long int) x;
 110     }
 111
 112   return sx ? -result : result;
 113 }
 114
 115 #endif /* _DOUBLE_IS_32BITS */