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
|
/*
* 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
#include <numpy/npy_math.h>
#include "_c99compat.h"
#include "cephes/dd_idefs.h"
double add_round_up(double a, double b)
{
double s, err;
if (sc_isnan(a) || sc_isnan(b)) {
return NPY_NAN;
}
s = two_sum(a, b, &err);
if (err > 0) {
/* fl(a + b) rounded down */
return npy_nextafter(s, NPY_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 (sc_isnan(a) || sc_isnan(b)) {
return NPY_NAN;
}
s = two_sum(a, b, &err);
if (err < 0) {
return npy_nextafter(s, -NPY_INFINITY);
}
else {
return s;
}
}
/* Helper code for testing _round.h. */
#if __STDC_VERSION__ >= 199901L
/* We have C99 */
#include <fenv.h>
#else
int fesetround(int round)
{
return -1;
}
int fegetround()
{
return -1;
}
#define FE_UPWARD -1
#define FE_DOWNWARD -1
#endif
#endif /* _round.h */
|
