source: trunk/libs/libmath/s_rint.c @ 567

Last change on this file since 567 was 469, checked in by alain, 6 years ago

1) Introduce the libsemaphore library.
2) Introduce a small libmath library, required by the "fft" application..
3) Introduce the multithreaded "fft" application.
4) Fix a bad synchronisation bug in the Copy-On-Write mechanism.

File size: 2.3 KB
Line 
1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/*
13 * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner)
14 */
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 "math.h"
27#include "math_private.h"
28
29static const double
30TWO52[2]={
31  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
32 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
33};
34
35double rint(double x)
36{
37        int32_t i0, j0, sx;
38        uint32_t i,i1;
39        double t;
40        /* We use w = x + 2^52; t = w - 2^52; trick to round x to integer.
41         * This trick requires that compiler does not optimize it
42         * by keeping intermediate result w in a register wider than double.
43         * Declaring w volatile assures that value gets truncated to double
44         * (unfortunately, it also forces store+load):
45         */
46        volatile double w;
47
48        EXTRACT_WORDS(i0,i1,x);
49        /* Unbiased exponent */
50        j0 = ((((uint32_t)i0) >> 20)&0x7ff)-0x3ff;
51
52        if (j0 > 51) {
53                //Why bother? Just returning x works too
54                //if (j0 == 0x400)  /* inf or NaN */
55                //      return x+x;
56                return x;  /* x is integral */
57        }
58
59        /* Sign */
60        sx = ((uint32_t)i0) >> 31;
61
62        if (j0<20) {
63            if (j0<0) { /* |x| < 1 */
64                if (((i0&0x7fffffff)|i1)==0) return x;
65                i1 |= (i0&0x0fffff);
66                i0 &= 0xfffe0000;
67                i0 |= ((i1|-i1)>>12)&0x80000;
68                SET_HIGH_WORD(x,i0);
69                w = TWO52[sx]+x;
70                t = w-TWO52[sx];
71                GET_HIGH_WORD(i0,t);
72                SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
73                return t;
74            } else {
75                i = (0x000fffff)>>j0;
76                if (((i0&i)|i1)==0) return x; /* x is integral */
77                i>>=1;
78                if (((i0&i)|i1)!=0) {
79                    if (j0==19) i1 = 0x40000000;
80                    else i0 = (i0&(~i))|((0x20000)>>j0);
81                }
82            }
83        } else {
84            i = ((uint32_t)(0xffffffff))>>(j0-20);
85            if ((i1&i)==0) return x;    /* x is integral */
86            i>>=1;
87            if ((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
88        }
89        INSERT_WORDS(x,i0,i1);
90        w = TWO52[sx]+x;
91        return w-TWO52[sx];
92}
Note: See TracBrowser for help on using the repository browser.