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
|
/*
* Functions for adding two double precision numbers with rounding to
* infinity or rounding to negative infinity without using <fenv.h>.
*/
#ifndef ROUND_H
#define ROUND_H
#if defined(__cplusplus)
#include <cmath>
using std::isnan;
#endif
#include <math.h>
/* Computes fl(a+b) and err(a+b). */
static inline double two_sum(double a, double b, double *err)
{
volatile double s = a + b;
volatile double c = s - a;
volatile double d = b - c;
volatile double e = s - c;
*err = (a - e) + d;
return s;
}
double add_round_up(double a, double b)
{
double s, err;
if (isnan(a) || isnan(b)) {
return NAN;
}
s = two_sum(a, b, &err);
if (err > 0) {
/* fl(a + b) rounded down */
return nextafter(s, INFINITY);
}
else {
/* fl(a + b) rounded up or didn't round */
return s;
}
}
double add_round_down(double a, double b)
{
double s, err;
if (isnan(a) || isnan(b)) {
return NAN;
}
s = two_sum(a, b, &err);
if (err < 0) {
return nextafter(s, -INFINITY);
}
else {
return s;
}
}
/* Helper code for testing _round.h. */
#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L) || defined(__cplusplus)
/* We have C99, or C++11 or higher; both have fenv.h */
#include <fenv.h>
#else
int fesetround(int round)
{
return -1;
}
int fegetround()
{
return -1;
}
#define FE_UPWARD -1
#define FE_DOWNWARD -1
#endif
/* SH4 is not C99 compliant, see https://github.com/scipy/scipy/issues/15584 */
#ifndef FE_UPWARD
#define FE_UPWARD -1
#endif
#ifndef FE_DOWNWARD
#define FE_DOWNWARD -1
#endif
#endif /* _round.h */
|
