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
|
/* adapted from Tango version 0.99.9, BSD Licensed
*/
/* 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 BIG_ENDIAN == 1
# 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 BIG_ENDIAN == 1
# 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
|
