Branch data     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 "cdefs-compat.h"
      14                 :             : //__FBSDID("$FreeBSD: src/lib/msun/src/s_rint.c,v 1.16 2008/02/22 02:30:35 das Exp $");
      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                 :             : #include <openlibm_math.h>
      28                 :             : 
      29                 :             : #include "math_private.h"
      30                 :             : 
      31                 :             : static const double
      32                 :             : TWO52[2]={
      33                 :             :   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
      34                 :             :  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
      35                 :             : };
      36                 :             : 
      37                 :             : OLM_DLLEXPORT double
      38                 :          23 : rint(double x)
      39                 :             : {
      40                 :             :         int32_t i0,j0,sx;
      41                 :             :         u_int32_t i,i1;
      42                 :             :         double w,t;
      43                 :          23 :         EXTRACT_WORDS(i0,i1,x);
      44                 :          23 :         sx = (i0>>31)&1;
      45                 :          23 :         j0 = ((i0>>20)&0x7ff)-0x3ff;
      46         [ +  + ]:          23 :         if(j0<20) {
      47         [ +  + ]:          18 :             if(j0<0) {
      48         [ +  + ]:           8 :                 if(((i0&0x7fffffff)|i1)==0) return x;
      49                 :           4 :                 i1 |= (i0&0x0fffff);
      50                 :           4 :                 i0 &= 0xfffe0000;
      51                 :           4 :                 i0 |= ((i1|-i1)>>12)&0x80000;
      52                 :           4 :                 SET_HIGH_WORD(x,i0);
      53                 :           4 :                 STRICT_ASSIGN(double,w,TWO52[sx]+x);
      54                 :           4 :                 t =  w-TWO52[sx];
      55                 :           4 :                 GET_HIGH_WORD(i0,t);
      56                 :           4 :                 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
      57                 :           4 :                 return t;
      58                 :             :             } else {
      59                 :          10 :                 i = (0x000fffff)>>j0;
      60         [ -  + ]:          10 :                 if(((i0&i)|i1)==0) return x; /* x is integral */
      61                 :          10 :                 i>>=1;
      62         [ -  + ]:          10 :                 if(((i0&i)|i1)!=0) {
      63                 :             :                     /*
      64                 :             :                      * Some bit is set after the 0.5 bit.  To avoid the
      65                 :             :                      * possibility of errors from double rounding in
      66                 :             :                      * w = TWO52[sx]+x, adjust the 0.25 bit to a lower
      67                 :             :                      * guard bit.  We do this for all j0<=51.  The
      68                 :             :                      * adjustment is trickiest for j0==18 and j0==19
      69                 :             :                      * since then it spans the word boundary.
      70                 :             :                      */
      71         [ #  # ]:           0 :                     if(j0==19) i1 = 0x40000000; else
      72         [ #  # ]:           0 :                     if(j0==18) i1 = 0x80000000; else
      73                 :           0 :                     i0 = (i0&(~i))|((0x20000)>>j0);
      74                 :             :                 }
      75                 :             :             }
      76         [ +  - ]:           5 :         } else if (j0>51) {
      77         [ +  - ]:           5 :             if(j0==0x400) return x+x;   /* inf or NaN */
      78                 :           0 :             else return x;              /* x is integral */
      79                 :             :         } else {
      80                 :           0 :             i = ((u_int32_t)(0xffffffff))>>(j0-20);
      81         [ #  # ]:           0 :             if((i1&i)==0) return x; /* x is integral */
      82                 :           0 :             i>>=1;
      83         [ #  # ]:           0 :             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
      84                 :             :         }
      85                 :          10 :         INSERT_WORDS(x,i0,i1);
      86                 :          10 :         STRICT_ASSIGN(double,w,TWO52[sx]+x);
      87                 :          10 :         return w-TWO52[sx];
      88                 :             : }
      89                 :             : 
      90                 :             : #if (LDBL_MANT_DIG == 53)
      91                 :             : openlibm_weak_reference(rint, rintl);
      92                 :             : #endif