Line data    Source code

       1             : /* @(#)s_rint.c 5.1 93/09/24 */
       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             : //#include <sys/cdefs.h>
      14             : //__FBSDID("$FreeBSD$");
      15             : 
      16             : /*
      17             :  * rint(x)
      18             :  * Return x rounded to integral value according to the prevailing
      19             :  * rounding mode.
      20             :  * Method:
      21             :  *      Using floating addition.
      22             :  * Exception:
      23             :  *      Inexact flag raised if x not equal to rint(x).
      24             :  */
      25             : 
      26             : #include <float.h>
      27             : 
      28             : #include "math_private.h"
      29             : 
      30             : static const double
      31             : TWO52[2]={
      32             :   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
      33             :  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
      34             : };
      35             : 
      36             : double
      37           0 : rint(double x)
      38             : {
      39             :         int32_t i0,j0,sx;
      40             :         u_int32_t i,i1;
      41             :         double w,t;
      42           0 :         EXTRACT_WORDS(i0,i1,x);
      43           0 :         sx = (i0>>31)&1;
      44           0 :         j0 = ((i0>>20)&0x7ff)-0x3ff;
      45           0 :         if(j0<20) {
      46           0 :             if(j0<0) {
      47           0 :                 if(((i0&0x7fffffff)|i1)==0) return x;
      48           0 :                 i1 |= (i0&0x0fffff);
      49           0 :                 i0 &= 0xfffe0000;
      50           0 :                 i0 |= ((i1|-i1)>>12)&0x80000;
      51           0 :                 SET_HIGH_WORD(x,i0);
      52           0 :                 STRICT_ASSIGN(double,w,TWO52[sx]+x);
      53           0 :                 t =  w-TWO52[sx];
      54           0 :                 GET_HIGH_WORD(i0,t);
      55           0 :                 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
      56           0 :                 return t;
      57             :             } else {
      58           0 :                 i = (0x000fffff)>>j0;
      59           0 :                 if(((i0&i)|i1)==0) return x; /* x is integral */
      60           0 :                 i>>=1;
      61           0 :                 if(((i0&i)|i1)!=0) {
      62             :                     /*
      63             :                      * Some bit is set after the 0.5 bit.  To avoid the
      64             :                      * possibility of errors from double rounding in
      65             :                      * w = TWO52[sx]+x, adjust the 0.25 bit to a lower
      66             :                      * guard bit.  We do this for all j0<=51.  The
      67             :                      * adjustment is trickiest for j0==18 and j0==19
      68             :                      * since then it spans the word boundary.
      69             :                      */
      70           0 :                     if(j0==19) i1 = 0x40000000; else
      71           0 :                     if(j0==18) i1 = 0x80000000; else
      72           0 :                     i0 = (i0&(~i))|((0x20000)>>j0);
      73             :                 }
      74             :             }
      75           0 :         } else if (j0>51) {
      76           0 :             if(j0==0x400) return x+x;   /* inf or NaN */
      77           0 :             else return x;              /* x is integral */
      78             :         } else {
      79           0 :             i = ((u_int32_t)(0xffffffff))>>(j0-20);
      80           0 :             if((i1&i)==0) return x; /* x is integral */
      81           0 :             i>>=1;
      82           0 :             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
      83             :         }
      84           0 :         INSERT_WORDS(x,i0,i1);
      85           0 :         STRICT_ASSIGN(double,w,TWO52[sx]+x);
      86           0 :         return w-TWO52[sx];
      87             : }