source: trunk/sys/libm/s_rint.c @ 128

Last change on this file since 128 was 1, checked in by alain, 8 years ago

First import

File size: 2.0 KB
Line 
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/*
15 * rint(x)
16 * Return x rounded to integral value according to the prevailing
17 * rounding mode.
18 * Method:
19 *      Using floating addition.
20 * Exception:
21 *      Inexact flag raised if x not equal to rint(x).
22 */
23
24#include <libm/fdlibm.h>
25
26#ifdef __STDC__
27static const double
28#else
29static double 
30#endif
31one = 1.0,
32TWO52[2]={
33  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
34 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
35};
36
37#ifdef __STDC__
38        double rint(double x)
39#else
40        double rint(x)
41        double x;
42#endif
43{
44        int i0,n0,j0,sx;
45        unsigned i,i1;
46        double w,t;
47        n0 = (*((int *)&one)>>29)^1;
48        i0 =  *(n0+(int*)&x);
49        sx = (i0>>31)&1;
50        i1 =  *(1-n0+(int*)&x);
51        j0 = ((i0>>20)&0x7ff)-0x3ff;
52        if(j0<20) {
53            if(j0<0) { 
54                if(((i0&0x7fffffff)|i1)==0) return x;
55                i1 |= (i0&0x0fffff);
56                i0 &= 0xfffe0000;
57                i0 |= ((i1|-i1)>>12)&0x80000;
58                *(n0+(int*)&x)=i0;
59                w = TWO52[sx]+x;
60                t =  w-TWO52[sx];
61                i0 = *(n0+(int*)&t);
62                *(n0+(int*)&t) = (i0&0x7fffffff)|(sx<<31);
63                return t;
64            } else {
65                i = (0x000fffff)>>j0;
66                if(((i0&i)|i1)==0) return x; /* x is integral */
67                i>>=1;
68                if(((i0&i)|i1)!=0) {
69                    if(j0==19) i1 = 0x40000000; else
70                    i0 = (i0&(~i))|((0x20000)>>j0);
71                }
72            }
73        } else if (j0>51) {
74            if(j0==0x400) return x+x;   /* inf or NaN */
75            else return x;              /* x is integral */
76        } else {
77            i = ((unsigned)(0xffffffff))>>(j0-20);
78            if((i1&i)==0) return x;     /* x is integral */
79            i>>=1;
80            if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
81        }
82        *(n0+(int*)&x) = i0;
83        *(1-n0+(int*)&x) = i1;
84        w = TWO52[sx]+x;
85        return w-TWO52[sx];
86}
Note: See TracBrowser for help on using the repository browser.