blob: 62d5936d9b68dad99b742b0d9e87395758c51345 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
#include <limits.h>
#include <math.h>
#include <endian.h>
typedef union
{
struct {
#if (__BYTE_ORDER == __BIG_ENDIAN)
unsigned long int hi;
unsigned long int lo;
#else
unsigned long int lo;
unsigned long int hi;
#endif
} words;
double dbl;
} DblInHex;
static const unsigned long int signMask = 0x80000000ul;
static const double twoTo52 = 4503599627370496.0;
/*******************************************************************************
* *
* The function round rounds its double argument to integral value *
* according to the "add half to the magnitude and truncate" rounding of *
* Pascal's Round function and FORTRAN's ANINT function and returns the *
* result in double format. This function signals inexact if an ordered *
* return value is not equal to the operand. *
* *
*******************************************************************************/
libm_hidden_proto(round)
double round ( double x )
{
DblInHex argument, OldEnvironment;
register double y, z;
register unsigned long int xHead;
register long int target;
argument.dbl = x;
xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead < 0x3fe00000ul )
/*******************************************************************************
* Is |x| < 0.5? *
*******************************************************************************/
{
if ( ( xHead | argument.words.lo ) != 0ul )
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 0.0 );
else
return ( -0.0 );
}
/*******************************************************************************
* Is 0.5 ² |x| < 1.0? *
*******************************************************************************/
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 1.0 );
else
return ( -1.0 );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{ // positive x
y = ( x + twoTo52 ) - twoTo52; // round at binary point
if ( y == x ) // exact case
return ( x );
z = x + 0.5; // inexact case
y = ( z + twoTo52 ) - twoTo52; // round at binary point
if ( y > z )
return ( y - 1.0 );
else
return ( y );
}
/*******************************************************************************
* Is x < 0? *
*******************************************************************************/
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary point
if ( y == x )
return ( x );
z = x - 0.5;
y = ( z - twoTo52 ) + twoTo52; // round at binary point
if ( y < z )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
libm_hidden_def(round)
|
