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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
//! Tools for ULP-based comparison of floating point numbers.
use std::mem;

/// Represents the result of an ULP-based comparison between two floating point numbers.
#[derive(Debug, Copy, Clone, PartialEq)]
pub enum UlpComparisonResult
{
    /// Signifies an exact match between two floating point numbers.
    ExactMatch,
    /// The difference in ULP between two floating point numbers.
    Difference(u64),
    /// The two floating point numbers have different signs,
    /// and cannot be compared in a meaningful way.
    IncompatibleSigns,
    /// One or both of the two floating point numbers is a NaN,
    /// in which case the ULP comparison is not meaningful.
    Nan
}

/// Floating point types for which two instances can be compared for Unit in the Last Place (ULP) difference.
///
/// Implementing this trait enables the usage of the `ulp` comparator in
/// [assert_matrix_eq!](../macro.assert_matrix_eq!.html) for the given type.
///
/// The definition here leverages the fact that for two adjacent floating point numbers,
/// their integer representations are also adjacent.
///
/// A somewhat accessible (but not exhaustive) guide on the topic is available in the popular article
/// [Comparing Floating Point Numbers, 2012 Edition]
/// (https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/).
///
/// Implementations for `f32` and `f64` are already available, and so users should not normally
/// need to implement this. In the case when a custom implementation is necessary,
/// please see the possible return values for [UlpComparisonResult](ulp/enum.UlpComparisonResult.html).
/// Otherwise, we can recommend to read the source code of the included `f32` and `f64` implementations.
pub trait Ulp {
    /// Returns the difference between two floating point numbers, measured in ULP.
    fn ulp_diff(a: &Self, b: &Self) -> UlpComparisonResult;
}

macro_rules! impl_float_ulp {
    ($ftype:ty, $itype:ty) => {
        impl Ulp for $ftype {
            fn ulp_diff(a: &Self, b: &Self) -> UlpComparisonResult {
                if a == b {
                    UlpComparisonResult::ExactMatch
                } else if a.is_nan() || b.is_nan() {
                    // ULP comparison does not make much sense for NaN
                    UlpComparisonResult::Nan
                } else if a.is_sign_positive() != b.is_sign_positive() {
                    // ULP is not meaningful when the signs of the two numbers differ
                    UlpComparisonResult::IncompatibleSigns
                } else {
                    // Otherwise, we compute the ULP diff as the difference of the signed integer representations
                    let a_int = unsafe { mem::transmute::<$ftype, $itype>(a.to_owned()) };
                    let b_int = unsafe { mem::transmute::<$ftype, $itype>(b.to_owned()) };
                    UlpComparisonResult::Difference((b_int - a_int).abs() as u64)
                }
            }
        }
    }
}

impl_float_ulp!(f32, i32);
impl_float_ulp!(f64, i64);

#[cfg(test)]
mod tests {
    use super::Ulp;
    use super::UlpComparisonResult;
    use std::mem;
    use std::{f32, f64};
    use quickcheck::TestResult;

    #[test]
    fn plus_minus_zero_is_exact_match_f32() {
        assert!(f32::ulp_diff(&0.0, &0.0) == UlpComparisonResult::ExactMatch);
        assert!(f32::ulp_diff(&-0.0, &-0.0) == UlpComparisonResult::ExactMatch);
        assert!(f32::ulp_diff(&0.0, &-0.0) == UlpComparisonResult::ExactMatch);
        assert!(f32::ulp_diff(&-0.0, &0.0) == UlpComparisonResult::ExactMatch);
    }

    #[test]
    fn plus_minus_zero_is_exact_match_f64() {
        assert!(f64::ulp_diff(&0.0, &0.0) == UlpComparisonResult::ExactMatch);
        assert!(f64::ulp_diff(&-0.0, &-0.0) == UlpComparisonResult::ExactMatch);
        assert!(f64::ulp_diff(&0.0, &-0.0) == UlpComparisonResult::ExactMatch);
        assert!(f64::ulp_diff(&-0.0, &0.0) == UlpComparisonResult::ExactMatch);
    }

    #[test]
    fn f32_double_nan() {
        assert!(f32::ulp_diff(&f32::NAN, &f32::NAN) == UlpComparisonResult::Nan);
    }

    #[test]
    fn f64_double_nan() {
        assert!(f64::ulp_diff(&f64::NAN, &f64::NAN) == UlpComparisonResult::Nan);
    }

    quickcheck! {
        fn property_exact_match_for_finite_f32_self_comparison(x: f32) -> TestResult {
            if x.is_finite() {
                TestResult::from_bool(f32::ulp_diff(&x, &x) == UlpComparisonResult::ExactMatch)
            } else {
                TestResult::discard()
            }
        }
    }

    quickcheck! {
        fn property_exact_match_for_finite_f64_self_comparison(x: f64) -> TestResult {
            if x.is_finite() {
                TestResult::from_bool(f64::ulp_diff(&x, &x) == UlpComparisonResult::ExactMatch)
            } else {
                TestResult::discard()
            }
        }
    }

    quickcheck! {
        fn property_recovers_ulp_diff_when_f32_constructed_from_i32(a: i32, b: i32) -> TestResult {
            if a == b {
                // Ignore self-comparisons, as it makes the below test have more complicated logic,
                // and moreover we test self-comparisons in another property.
                return TestResult::discard();
            }

            let x = unsafe { mem::transmute::<i32, f32>(a) };
            let y = unsafe { mem::transmute::<i32, f32>(b) };

            // Discard the input if it's non-finite or has different signs
            if x.is_finite() && y.is_finite() && x.signum() == y.signum() {
                TestResult::from_bool(f32::ulp_diff(&x, &y) == UlpComparisonResult::Difference((b - a).abs() as u64))
            } else {
                TestResult::discard()
            }
        }
    }

    quickcheck! {
        fn property_recovers_ulp_diff_when_f64_constructed_from_i64(a: i64, b: i64) -> TestResult {
            if a == b {
                // Ignore self-comparisons, as it makes the below test have more complicated logic,
                // and moreover we test self-comparisons in another property.
                return TestResult::discard();
            }

            let x = unsafe { mem::transmute::<i64, f64>(a) };
            let y = unsafe { mem::transmute::<i64, f64>(b) };

            // Discard the input if it's non-finite or has different signs
            if x.is_finite() && y.is_finite() && x.signum() == y.signum() {
                TestResult::from_bool(f64::ulp_diff(&x, &y) == UlpComparisonResult::Difference((b - a).abs() as u64))
            } else {
                TestResult::discard()
            }
        }
    }

    quickcheck! {
        fn property_f32_incompatible_signs_yield_corresponding_enum_value(x: f32, y: f32) -> TestResult {
            if x.signum() == y.signum() {
                TestResult::discard()
            } else if x.is_nan() || y.is_nan() {
                TestResult::discard()
            } else {
                TestResult::from_bool(f32::ulp_diff(&x, &y) == UlpComparisonResult::IncompatibleSigns)
            }
        }
    }

    quickcheck! {
        fn property_f64_incompatible_signs_yield_corresponding_enum_value(x: f64, y: f64) -> TestResult {
            if x.signum() == y.signum() {
                TestResult::discard()
            } else if x.is_nan() || y.is_nan() {
                TestResult::discard()
            } else {
                TestResult::from_bool(f64::ulp_diff(&x, &y) == UlpComparisonResult::IncompatibleSigns)
            }
        }
    }

    quickcheck! {
        fn property_f32_nan_gives_nan_enum_value(x: f32) -> bool {
            f32::ulp_diff(&f32::NAN, &x) == UlpComparisonResult::Nan
            && f32::ulp_diff(&x, &f32::NAN) == UlpComparisonResult::Nan
        }
    }

    quickcheck! {
        fn property_f64_nan_gives_nan_enum_value(x: f64) -> bool {
            f64::ulp_diff(&f64::NAN, &x) == UlpComparisonResult::Nan
            && f64::ulp_diff(&x, &f64::NAN) == UlpComparisonResult::Nan
        }
    }
}