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
156
157
158
159
160
161
162
163
164
|
/* adapted from Tango version 0.99.9, BSD Licensed
*/
/* Define WORDS_BIGENDIAN on big-endian platforms; otherwise, assume
* little-endian.
*
* Endianness detection modified from http://esr.ibiblio.org/?p=5095 and
* https://gist.github.com/panzi/6856583 .
*
* In the subsequent code we rely on endianness being the same for integers
* and floats; this is true for modern systems but false in a few historical
* machines and some old ARM processors; see
* http://en.wikipedia.org/wiki/Endianness#Floating-point_and_endianness .
*/
#if defined(__BYTE_ORDER__)
/* (1) Try to use the GCC/Clang __BYTE_ORDER__ defines */
# if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__
# define WORDS_BIGENDIAN 1
# elif __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
# undef WORDS_BIGENDIAN
# else
# error byte order not supported
# endif
/* (2) Otherwise, try the __BIG_ENDIAN__ and __LITTLE_ENDIAN__ defines */
#elif defined(__BIG_ENDIAN__)
# define WORDS_BIGENDIAN 1
#elif defined(__LITTLE_ENDIAN__)
# undef WORDS_BIGENDIAN
#else
/* (3) As a last resort, try platform-specific fallbacks: */
# if defined(_WIN16) || defined(_WIN32) || defined(_WIN64)
/* (a) Assume Windows is little endian (http://stackoverflow.com/a/6449581 ) */
# undef WORDS_BIGENDIAN
# else
/* (b) Otherwise, use endian.h */
# if (defined(__DragonFly__) || defined(__FreeBSD__) || defined(__NetBSD__) \
|| defined(__OpenBSD__))
# include <sys/endian.h>
# else
# include <endian.h>
# endif
# if __BYTE_ORDER == __BIG_ENDIAN
# define WORDS_BIGENDIAN 1
# elif __BYTE_ORDER == __LITTLE_ENDIAN
# undef WORDS_BIGENDIAN
# else
# error platform not supported: unable to determine endianess
# endif
# endif
#endif
/* REAL_EXPMASK is a ushort mask to select the exponent portion (without sign)
* REAL_SIGNMASK is a ushort mask to select the sign bit.
* REAL_EXPPOS_SHORT is the index of the exponent when represented as a uint16_t array.
* REAL_SIGNPOS_BYTE is the index of the sign when represented as a uint8_t array.
* REAL_RECIP_EPSILON is the value such that
* (smallest_denormal) * REAL_RECIP_EPSILON == REAL_MIN_NORMAL
*/
#define REAL_RECIP_EPSILON (1 / REAL_EPSILON)
#if REAL_MANT_DIG == 24
# define REAL_EXPMASK ((uint16_t) 0x7F80)
# define REAL_SIGNMASK ((uint16_t) 0x8000)
# define REAL_EXPBIAS ((uint16_t) 0x3F00)
# define REAL_EXPBIAS_INT32 ((uint32_t) 0x7F800000)
# define REAL_MANTISSAMASK_INT32 ((uint32_t) 0x007FFFFF)
# if WORDS_BIGENDIAN
# define REAL_EXPPOS_INT16 0
# else
# define REAL_EXPPOS_INT16 1
# endif
#
#elif REAL_MANT_DIG == 53 /* double */
# define REAL_EXPMASK ((uint16_t) 0x7FF0)
# define REAL_SIGNMASK ((uint16_t) 0x8000)
# define REAL_EXPBIAS ((uint16_t) 0x3FE0)
# define REAL_EXPBIAS_INT32 ((uint32_t) 0x7FF00000)
# define REAL_MANTISSAMASK_INT32 ((uint32_t) 0x000FFFFF); /* for the MSB only */
# if WORDS_BIGENDIAN
# define REAL_EXPPOS_INT16 0
# define REAL_SIGNPOS_BYTE 0
# else
# define REAL_EXPPOS_INT16 3
# define REAL_SIGNPOS_BYTE 7
# endif
#endif
int
FEQREL (REAL x, REAL y)
{
/* Public Domain. Original Author: Don Clugston, 18 Aug 2005.
* Ported to C by Patrick Perry, 26 Feb 2010.
*/
if (x == y) return REAL_MANT_DIG; /* ensure diff!= 0, cope with INF. */
REAL diff = REAL_ABS(x - y);
union { REAL r; uint16_t w[sizeof(REAL)/2]; } pa = { x };
union { REAL r; uint16_t w[sizeof(REAL)/2]; } pb = { y };
union { REAL r; uint16_t w[sizeof(REAL)/2]; } pd = { diff };
/* The difference in abs(exponent) between x or y and abs(x-y)
* is equal to the number of significand bits of x which are
* equal to y. If negative, x and y have different exponents.
* If positive, x and y are equal to 'bitsdiff' bits.
* AND with 0x7FFF to form the absolute value.
* To avoid out-by-1 errors, we subtract 1 so it rounds down
* if the exponents were different. This means 'bitsdiff' is
* always 1 lower than we want, except that if bitsdiff==0,
* they could have 0 or 1 bits in common.
*/
#if REAL_MANT_DIG == 53 /* double */
int bitsdiff = (( ((pa.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)
+ (pb.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)
- ((uint16_t) 0x8000 - REAL_EXPMASK)) >> 1)
- (pd.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)) >> 4;
#elif REAL_MANT_DIG == 24 /* float */
int bitsdiff = (( ((pa.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)
+ (pb.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)
- ((uint16_t) 0x8000 - REAL_EXPMASK)) >> 1)
- (pd.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)) >> 7;
#else
# error unsuported floating-point mantissa size
#endif
if ((pd.w[REAL_EXPPOS_INT16] & REAL_EXPMASK) == 0) {
/* Difference is denormal
* For denormals, we need to add the number of zeros that
* lie at the start of diff's significand.
* We do this by multiplying by 2^REAL_MANT_DIG
*/
pd.r *= REAL_RECIP_EPSILON;
#if REAL_MANT_DIG == 53 /* double */
return (bitsdiff + REAL_MANT_DIG
- (pd.w[REAL_EXPPOS_INT16] >> 4));
#elif REAL_MANT_DIG == 24 /* float */
return (bitsdiff + REAL_MANT_DIG
- (pd.w[REAL_EXPPOS_INT16] >> 7));
#else
# error unsuported floating-point mantissa size
#endif
}
if (bitsdiff > 0)
return bitsdiff + 1; /* add the 1 we subtracted before */
/* Avoid out-by-1 errors when factor is almost 2. */
return (bitsdiff == 0
&& !((pa.w[REAL_EXPPOS_INT16]
^ pb.w[REAL_EXPPOS_INT16]) & REAL_EXPMASK)) ? 1 : 0;
}
#undef REAL_RECIP_EPSILON
#undef REAL_EXPMASK
#undef REAL_SIGNMASK
#undef REAL_EXPBIAS
#undef REAL_EXPBIAS_INT32
#undef REAL_MANTISSAMASK_INT32
#undef REAL_EXPPOS_INT16
#undef REAL_SIGNPOS_BYTE
|
