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
|
#include <float.h>
#include <math.h>
#include <stdint.h>
#define REAL double
#define REAL_ABS fabs
#define REAL_MIN_NORMAL DBL_MIN
#define REAL_EPSILON DBL_EPSILON
#define REAL_MAX DBL_MAX
#define REAL_MANT_DIG DBL_MANT_DIG
#define FEQREL feqrel
#include "feqrel_source.c"
union double_t {
double d;
uint64_t w;
};
int
identical (double x, double y)
{
union double_t ux = { x };
union double_t uy = { y };
return ux.w == uy.w;
}
double
copysign (double x, double y)
{
union double_t ux = { x };
union double_t uy = { y };
union double_t uz;
uint64_t val = ux.w & 0x7FFFFFFFFFFFFFFFULL;
uint64_t sign = uy.w & 0x8000000000000000ULL;
uz.w = sign | val;
return uz.d;
}
/* ported from tango/math/IEEE.d */
double
hsi7_nextup (double x)
{
union double_t ps = { x };
if ((ps.w & 0x7FF0000000000000ULL) == 0x7FF0000000000000ULL) {
/* First, deal with NANs and infinity */
if (x == -INFINITY) return -REAL_MAX;
return x; // +INF and NAN are unchanged.
}
if (ps.w & 0x8000000000000000ULL) { /* Negative number */
if (ps.w == 0x8000000000000000ULL) { /* it was negative zero */
ps.w = 0x0000000000000001ULL; /* change to smallest subnormal */
return ps.d;
}
--ps.w;
} else { /* Positive number */
++ps.w;
}
return ps.d;
}
/* ported from tango/math/IEEE.d */
double
hsi7_nextdown (double x)
{
return -hsi7_nextup(-x);
}
/* ported from tango/math/IEEE.d */
double
ieeemean (double x, double y)
{
if (!((x>=0 && y>=0) || (x<=0 && y<=0))) return NAN;
union double_t ul;
union double_t xl = { x };
union double_t yl = { y };
uint64_t m = ( (xl.w & 0x7FFFFFFFFFFFFFFFULL)
+ (yl.w & 0x7FFFFFFFFFFFFFFFULL) ) >> 1;
m |= (xl.w & 0x8000000000000000ULL);
ul.w = m;
return ul.d;
}
double
mknan (uint64_t payload)
{
union double_t x = { NAN };
/* get sign, exponent, and quiet bit from NAN */
x.w &= 0xFFF8000000000000ULL;
/* ignore sign, exponent, and quiet bit in payload */
payload &= 0x0007FFFFFFFFFFFFULL;
x.w |= payload;
return x.d;
}
uint64_t
getnan (double x)
{
union double_t payload = { x };
/* clear sign, exponent, and quiet bit */
payload.w &= 0x0007FFFFFFFFFFFFULL;
return payload.w;
}
|
