Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1.  
  2. /* @(#)s_rint.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. <<rint>>, <<rintf>>--round to integer
  16. INDEX
  17.         rint
  18. INDEX
  19.         rintf
  20.  
  21. ANSI_SYNOPSIS
  22.         #include <math.h>
  23.         double rint(double <[x]>);
  24.         float rintf(float <[x]>);
  25.  
  26. DESCRIPTION
  27.         The <<rint>> functions round their argument to an integer value in
  28.         floating-point format, using the current rounding direction.  They
  29.         raise the "inexact" floating-point exception if the result differs
  30.         in value from the argument.  See the <<nearbyint>> functions for the
  31.         same function with the "inexact" floating-point exception never being
  32.         raised.  Newlib does not directly support floating-point exceptions.
  33.         The <<rint>> functions are written so that the "inexact" exception is
  34.         raised in hardware implementations that support it, even though Newlib
  35.         does not provide access.
  36.  
  37. RETURNS
  38. <[x]> rounded to an integral value, using the current rounding direction.
  39.  
  40. PORTABILITY
  41. ANSI C, POSIX
  42.  
  43. SEEALSO
  44. <<nearbyint>>, <<round>>
  45.  
  46. */
  47.  
  48. /*
  49.  * rint(x)
  50.  * Return x rounded to integral value according to the prevailing
  51.  * rounding mode.
  52.  * Method:
  53.  *      Using floating addition.
  54.  *      Whenever a fraction is present, if the second or any following bit after
  55.  *      the radix point is set, limit to the second radix point to avoid
  56.  *      possible double rounding in the TWO52 +- steps (in case guard bits are
  57.  *      used).  Specifically, if have any, chop off bits past the 2nd place and
  58.  *      set the second place.
  59.  *      (e.g.   2.0625=0b10.0001 => 0b10.01=2.25;
  60.  *              2.3125=0b10.011  => 0b10.01=2.25;
  61.  *              1.5625= 0b1.1001 =>  0b1.11=1.75;
  62.  *              1.9375= 0b1.1111 =>  0b1.11=1.75.
  63.  *      Pseudo-code:  if(x.frac & ~0b0.10) x.frac = (x.frac & 0b0.11) | 0b0.01;).
  64.  * Exception:
  65.  *      Inexact flag raised if x not equal to rint(x).
  66.  */
  67.  
  68. #include "fdlibm.h"
  69.  
  70. #ifndef _DOUBLE_IS_32BITS
  71.  
  72. #ifdef __STDC__
  73. static const double
  74. #else
  75. static double
  76. #endif
  77. TWO52[2]={
  78.   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  79.  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  80. };
  81.  
  82. #ifdef __STDC__
  83.         double rint(double x)
  84. #else
  85.         double rint(x)
  86.         double x;
  87. #endif
  88. {
  89.         __int32_t i0,j0,sx;
  90.         __uint32_t i,i1;
  91.         double t;
  92.         volatile double w;
  93.         EXTRACT_WORDS(i0,i1,x);
  94.         sx = (i0>>31)&1;                /* sign */
  95.         j0 = ((i0>>20)&0x7ff)-0x3ff;    /* exponent */
  96.         if(j0<20) {                     /* no integral bits in LS part */
  97.             if(j0<0) {                  /* x is fractional or 0 */
  98.                 if(((i0&0x7fffffff)|i1)==0) return x;   /* x == 0 */
  99.                 i1 |= (i0&0x0fffff);
  100.                 i0 &= 0xfffe0000;
  101.                 i0 |= ((i1|-i1)>>12)&0x80000;
  102.                 SET_HIGH_WORD(x,i0);
  103.                 w = TWO52[sx]+x;
  104.                 t =  w-TWO52[sx];
  105.                 GET_HIGH_WORD(i0,t);
  106.                 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
  107.                 return t;
  108.             } else {                    /* x has integer and maybe fraction */
  109.                 i = (0x000fffff)>>j0;
  110.                 if(((i0&i)|i1)==0) return x; /* x is integral */
  111.                 i>>=1;
  112.                 if(((i0&i)|i1)!=0) {
  113.                     /* 2nd or any later bit after radix is set */
  114.                     if(j0==19) i1 = 0x80000000; else i1 = 0;
  115.                     i0 = (i0&(~i))|((0x40000)>>j0);
  116.                 }
  117.             }
  118.         } else if (j0>51) {
  119.             if(j0==0x400) return x+x;   /* inf or NaN */
  120.             else return x;              /* x is integral */
  121.         } else {
  122.             i = ((__uint32_t)(0xffffffff))>>(j0-20);
  123.             if((i1&i)==0) return x;     /* x is integral */
  124.             i>>=1;
  125.             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
  126.         }
  127.         INSERT_WORDS(x,i0,i1);
  128.         w = TWO52[sx]+x;
  129.         return w-TWO52[sx];
  130. }
  131.  
  132. #endif /* _DOUBLE_IS_32BITS */
  133.