Subversion Repositories shark

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 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
#ifndef lint
14
static char rcsid[] = "$\Id: s_rint.c,v 1.3.2.1 1997/02/23 11:03:20 joerg Exp $";
15
#endif
16
 
17
/*
18
 * rint(x)
19
 * Return x rounded to integral value according to the prevailing
20
 * rounding mode.
21
 * Method:
22
 *      Using floating addition.
23
 * Exception:
24
 *      Inexact flag raised if x not equal to rint(x).
25
 */
26
 
27
#include "math.h"
28
#include "math_private.h"
29
 
30
/*
31
 * TWO23 is long double instead of double to avoid a bug in gcc.  Without
32
 * this, gcc thinks that TWO23[sx]+x and w-TWO23[sx] already have double
33
 * precision and doesn't clip them to double precision when they are
34
 * assigned and returned.  Use long double even in the !__STDC__ case in
35
 * case this is compiled with gcc -traditional.
36
 */
37
#ifdef __STDC__
38
static const long double
39
#else
40
static long double
41
#endif
42
TWO52[2]={
43
  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
44
 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
45
};
46
 
47
#ifdef __STDC__
48
        double __generic_rint(double x)
49
#else
50
        double __generic_rint(x)
51
        double x;
52
#endif
53
{
54
        int32_t i0,j0,sx;
55
        u_int32_t i,i1;
56
        double w,t;
57
        EXTRACT_WORDS(i0,i1,x);
58
        sx = (i0>>31)&1;
59
        j0 = ((i0>>20)&0x7ff)-0x3ff;
60
        if(j0<20) {
61
            if(j0<0) {
62
                if(((i0&0x7fffffff)|i1)==0) return x;
63
                i1 |= (i0&0x0fffff);
64
                i0 &= 0xfffe0000;
65
                i0 |= ((i1|-i1)>>12)&0x80000;
66
                SET_HIGH_WORD(x,i0);
67
                w = TWO52[sx]+x;
68
                t =  w-TWO52[sx];
69
                GET_HIGH_WORD(i0,t);
70
                SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
71
                return t;
72
            } else {
73
                i = (0x000fffff)>>j0;
74
                if(((i0&i)|i1)==0) return x; /* x is integral */
75
                i>>=1;
76
                if(((i0&i)|i1)!=0) {
77
                    if(j0==19) i1 = 0x40000000; else
78
                    i0 = (i0&(~i))|((0x20000)>>j0);
79
                }
80
            }
81
        } else if (j0>51) {
82
            if(j0==0x400) return x+x;   /* inf or NaN */
83
            else return x;              /* x is integral */
84
        } else {
85
            i = ((u_int32_t)(0xffffffff))>>(j0-20);
86
            if((i1&i)==0) return x;     /* x is integral */
87
            i>>=1;
88
            if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
89
        }
90
        INSERT_WORDS(x,i0,i1);
91
        w = TWO52[sx]+x;
92
        return w-TWO52[sx];
93
}