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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
#ifndef _ROUND_INTERNAL_H
/*
* round_internal.h
*
* $Id$
*
* Provides a generic implementation of the numerical rounding
* algorithm, which is shared by all functions in the `round()',
* `lround()' and `llround()' families.
*
* Written by Keith Marshall <keithmarshall@users.sourceforge.net>
*
* This is free software. You may redistribute and/or modify it as you
* see fit, without restriction of copyright.
*
* This software is provided "as is", in the hope that it may be useful,
* but WITHOUT WARRANTY OF ANY KIND, not even any implied warranty of
* MERCHANTABILITY, nor of FITNESS FOR ANY PARTICULAR PURPOSE. At no
* time will the author accept any form of liability for any damages,
* however caused, resulting from the use of this software.
*
*/
#define _ROUND_INTERNAL_H
#include <math.h>
#include <fenv.h>
#define TYPE_PASTE( NAME, TYPE ) NAME##TYPE
#define INPUT_TYPE INPUT_TYPEDEF( FUNCTION )
#define INPUT_TYPEDEF( FUNCTION ) TYPE_PASTE( FUNCTION, _input_type )
/*
* The types for the formal parameter, to each of the derived functions.
*/
#define round_input_type double
#define roundf_input_type float
#define roundl_input_type long double
#define lround_input_type double
#define lroundf_input_type float
#define lroundl_input_type long double
#define llround_input_type double
#define llroundf_input_type float
#define llroundl_input_type long double
#define RETURN_TYPE RETURN_TYPEDEF( FUNCTION )
#define RETURN_TYPEDEF( FUNCTION ) TYPE_PASTE( FUNCTION, _return_type )
/*
* The types for the return value, from each of the derived functions.
*/
#define round_return_type double
#define roundf_return_type float
#define roundl_return_type long double
#define lround_return_type long
#define lroundf_return_type long
#define lroundl_return_type long
#define llround_return_type long long
#define llroundf_return_type long long
#define llroundl_return_type long long
#define MAX_RETURN_VALUE RETURN_MAX( FUNCTION )
#define RETURN_MAX( FUNCTION ) TYPE_PASTE( FUNCTION, _return_max )
/*
* The maximum values which may be returned by each of the derived functions
* in the `lround' or the `llround' families.
*/
#define lround_return_max LONG_MAX
#define lroundf_return_max LONG_MAX
#define lroundl_return_max LONG_MAX
#define llround_return_max LLONG_MAX
#define llroundf_return_max LLONG_MAX
#define llroundl_return_max LLONG_MAX
#define MIN_RETURN_VALUE RETURN_MIN( FUNCTION )
#define RETURN_MIN( FUNCTION ) TYPE_PASTE( FUNCTION, _return_min )
/*
* The minimum values which may be returned by each of the derived functions
* in the `lround' or the `llround' families.
*/
#define lround_return_min LONG_MIN
#define lroundf_return_min LONG_MIN
#define lroundl_return_min LONG_MIN
#define llround_return_min LLONG_MIN
#define llroundf_return_min LLONG_MIN
#define llroundl_return_min LLONG_MIN
#define REF_VALUE( VALUE ) REF_TYPE( FUNCTION, VALUE )
#define REF_TYPE( FUNC, VAL ) TYPE_PASTE( FUNC, _ref )( VAL )
/*
* Macros for expressing constant values of the appropriate data type,
* for use in each of the derived functions.
*/
#define round_ref( VALUE ) VALUE
#define lround_ref( VALUE ) VALUE
#define llround_ref( VALUE ) VALUE
#define roundl_ref( VALUE ) TYPE_PASTE( VALUE, L )
#define lroundl_ref( VALUE ) TYPE_PASTE( VALUE, L )
#define llroundl_ref( VALUE ) TYPE_PASTE( VALUE, L )
#define roundf_ref( VALUE ) TYPE_PASTE( VALUE, F )
#define lroundf_ref( VALUE ) TYPE_PASTE( VALUE, F )
#define llroundf_ref( VALUE ) TYPE_PASTE( VALUE, F )
static __inline__
INPUT_TYPE __attribute__(( always_inline )) round_internal( INPUT_TYPE x )
#define ROUND_MODES ( FE_TONEAREST | FE_UPWARD | FE_DOWNWARD | FE_TOWARDZERO )
{
/* Generic helper function, for rounding of the input parameter value to
* the nearest integer value.
*/
INPUT_TYPE z;
unsigned short saved_CW, tmp_required_CW;
/* Rounding method suggested by Danny Smith <dannysmith@users.sf.net>
*
* Save the FPU control word state, set rounding mode TONEAREST, round the
* input value, then restore the original FPU control word state.
*/
__asm__( "fnstcw %0;" : "=m"( saved_CW ));
tmp_required_CW = ( saved_CW & ~ROUND_MODES ) | FE_TONEAREST;
__asm__( "fldcw %0;" :: "m"( tmp_required_CW ));
__asm__( "frndint;" : "=t"( z ) : "0"( x ));
__asm__( "fldcw %0;" :: "m"( saved_CW ));
/* We now have a possible rounded value; unfortunately the FPU uses the
* `round-to-even' rule for exact mid-way cases, where both C99 and POSIX
* require us to always round away from zero, so we need to adjust those
* mid-way cases which the FPU rounded in the wrong direction.
*
* Correction method suggested by Greg Chicares <gchicares@sbcglobal.net>
*/
return x < REF_VALUE( 0.0 )
? /*
* For negative input values, an incorrectly rounded value will be
* exactly 0.5 greater than the original value; when we find such an
* exact rounding offset, we must subtract an additional 1.0 from the
* rounded result, otherwise we return the rounded result unchanged.
*/
z - x == REF_VALUE( 0.5 ) ? z - REF_VALUE( 1.0 ) : z
: /* For positive input values, an incorrectly rounded value will be
* exactly 0.5 less than the original value; when we find such an exact
* rounding offset, we must add an additional 1.0 to the rounded result,
* otherwise we return the rounded result unchanged.
*/
x - z == REF_VALUE( 0.5 ) ? z + REF_VALUE( 1.0 ) : z;
}
#endif /* !defined _ROUND_INTERNAL_H: $RCSfile$: end of file */
|
