newlib/libm/common/s_lrint.c

   1
   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 FUNCTION
  15 <<lrint>>, <<lrintf>>, <<llrint>>, <<llrintf>>---round to integer
  16 INDEX
  17         lrint
  18 INDEX
  19         lrintf
  20 INDEX
  21         llrint
  22 INDEX
  23         llrintf
  24
  25 SYNOPSIS
  26         #include <math.h>
  27         long int lrint(double <[x]>);
  28         long int lrintf(float <[x]>);
  29         long long int llrint(double <[x]>);
  30         long long int llrintf(float <[x]>);
  31
  32 DESCRIPTION
  33 The <<lrint>> and <<llrint>> functions round their argument to the nearest
  34 integer value, using the current rounding direction.  If the rounded value is
  35 outside the range of the return type, the numeric result is unspecified.  A
  36 range error may occur if the magnitude of <[x]> is too large.
  37 The "inexact" floating-point exception is raised in implementations that
  38 support it when the result differs in value from the argument (i.e., when
  39 a fraction actually has been truncated).
  40
  41 RETURNS
  42 <[x]> rounded to an integral value, using the current rounding direction.
  43
  44 SEEALSO
  45 <<lround>>
  46
  47 PORTABILITY
  48 ANSI C, POSIX
  49
  50 */
  51
  52 /*
  53  * lrint(x)
  54  * Return x rounded to integral value according to the prevailing
  55  * rounding mode.
  56  * Method:
  57  *      Using floating addition.
  58  * Exception:
  59  *      Inexact flag raised if x not equal to lrint(x).
  60  */
  61
  62 #include "fdlibm.h"
  63
  64 #ifndef _DOUBLE_IS_32BITS
  65
  66 #ifdef __STDC__
  67 static const double
  68 #else
  69 static double
  70 #endif
  71
  72 /* Adding a double, x, to 2^52 will cause the result to be rounded based on
  73    the fractional part of x, according to the implementation's current rounding
  74    mode.  2^52 is the smallest double that can be represented using all 52 significant
  75    digits. */
  76 TWO52[2]={
  77   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  78  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  79 };
  80
  81 #ifdef __STDC__
  82         long int lrint(double x)
  83 #else
  84         long int lrint(x)
  85         double x;
  86 #endif
  87 {
  88   __int32_t i0,j0,sx;
  89   __uint32_t i1;
  90   double t;
  91   volatile double w;
  92   long int result;
  93
  94   EXTRACT_WORDS(i0,i1,x);
  95
  96   /* Extract sign bit. */
  97   sx = (i0>>31)&1;
  98
  99   /* Extract exponent field. */
 100   j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
 101   /* j0 in [-1023,1024] */
 102
 103   if(j0 < 20)
 104     {
 105       /* j0 in [-1023,19] */
 106       w = TWO52[sx] + x;
 107       t = w - TWO52[sx];
 108       GET_HIGH_WORD(i0, t);
 109       /* Detect the all-zeros representation of plus and
 110          minus zero, which fails the calculation below. */
 111       if ((i0 & ~(1L << 31)) == 0)
 112           return 0;
 113       j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
 114       i0 &= 0x000fffff;
 115       i0 |= 0x00100000;
 116       result = (j0 < 0 ? 0 : i0 >> (20 - j0));
 117     }
 118   else if (j0 < (int)(8 * sizeof (long int)) - 1)
 119     {
 120       /* 32bit return: j0 in [20,30] */
 121       /* 64bit return: j0 in [20,62] */
 122       if (j0 >= 52)
 123         /* 64bit return: j0 in [52,62] */
 124         /* 64bit return: left shift amt in [32,42] */
 125         result = ((long int) ((i0 & 0x000fffff) | 0x00100000) << (j0 - 20)) |
 126                 /* 64bit return: right shift amt in [0,10] */
 127                    ((long int) i1 << (j0 - 52));
 128       else
 129         {
 130           /* 32bit return: j0 in [20,30] */
 131           /* 64bit return: j0 in [20,51] */
 132           w = TWO52[sx] + x;
 133           t = w - TWO52[sx];
 134           EXTRACT_WORDS (i0, i1, t);
 135           j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
 136           i0 &= 0x000fffff;
 137           i0 |= 0x00100000;
 138           /* After round:
 139            * 32bit return: j0 in [20,31];
 140            * 64bit return: j0 in [20,52] */
 141           /* 32bit return: left shift amt in [0,11] */
 142           /* 64bit return: left shift amt in [0,32] */
 143           /* ***32bit return: right shift amt in [32,21] */
 144           /* ***64bit return: right shift amt in [32,0] */
 145           result = ((long int) i0 << (j0 - 20))
 146                         | SAFE_RIGHT_SHIFT (i1, (52 - j0));
 147         }
 148     }
 149   else
 150     {
 151       return (long int) x;
 152     }
 153
 154   return sx ? -result : result;
 155 }
 156
 157 #endif /* _DOUBLE_IS_32BITS */