Subversion Repositories Kolibri OS

Rev

Rev 4872 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  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. ANSI_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.       if(j0 < -1)
  107.         return 0;
  108.       else
  109.         {
  110.           /* j0 in [0,19] */
  111.           /* shift amt in [0,19] */
  112.           w = TWO52[sx] + x;
  113.           t = w - TWO52[sx];
  114.           GET_HIGH_WORD(i0, t);
  115.           /* Detect the all-zeros representation of plus and
  116.              minus zero, which fails the calculation below. */
  117.           if ((i0 & ~(1L << 31)) == 0)
  118.               return 0;
  119.           /* After round:  j0 in [0,20] */
  120.           j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
  121.           i0 &= 0x000fffff;
  122.           i0 |= 0x00100000;
  123.           /* shift amt in [20,0] */
  124.           result = i0 >> (20 - j0);
  125.         }
  126.     }
  127.   else if (j0 < (int)(8 * sizeof (long int)) - 1)
  128.     {
  129.       /* 32bit return: j0 in [20,30] */
  130.       /* 64bit return: j0 in [20,62] */
  131.       if (j0 >= 52)
  132.         /* 64bit return: j0 in [52,62] */
  133.         /* 64bit return: left shift amt in [32,42] */
  134.         result = ((long int) ((i0 & 0x000fffff) | 0x0010000) << (j0 - 20)) |
  135.                 /* 64bit return: right shift amt in [0,10] */
  136.                    (i1 << (j0 - 52));
  137.       else
  138.         {
  139.           /* 32bit return: j0 in [20,30] */
  140.           /* 64bit return: j0 in [20,51] */
  141.           w = TWO52[sx] + x;
  142.           t = w - TWO52[sx];
  143.           EXTRACT_WORDS (i0, i1, t);
  144.           j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
  145.           i0 &= 0x000fffff;
  146.           i0 |= 0x00100000;
  147.           /* After round:
  148.            * 32bit return: j0 in [20,31];
  149.            * 64bit return: j0 in [20,52] */
  150.           /* 32bit return: left shift amt in [0,11] */
  151.           /* 64bit return: left shift amt in [0,32] */
  152.           /* ***32bit return: right shift amt in [32,21] */
  153.           /* ***64bit return: right shift amt in [32,0] */
  154.           result = ((long int) i0 << (j0 - 20))
  155.                         | SAFE_RIGHT_SHIFT (i1, (52 - j0));
  156.         }
  157.     }
  158.   else
  159.     {
  160.       return (long int) x;
  161.     }
  162.  
  163.   return sx ? -result : result;
  164. }
  165.  
  166. #endif /* _DOUBLE_IS_32BITS */
  167.