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

use super::sort::{sort_by};
use num::{Float,zero};
use unreachable::unreachable;

/// Sorts floating point numbers efficiently.
/// The ordering used is
/// | -inf | < 0 | -0 | +0 | > 0 | +inf | NaN |
pub fn sort_floats<T: Float>(v: &mut [T]) {
    /*
     * We don't have hardware support for a total order on floats. NaN is not
     * smaller or greater than any number. We want NaNs to be last, so we could
     * just use is_nan() in the comparison function. It turns out that this hurts
     * performance a lot, and in most cases we probably don't have any NaNs anyway.
     * 
     * The solution is to first move all NaNs to the end of the array, and then
     * sort the remainder with efficient comparisons. After sorting, the zeros
     * might be in the wrong order since -0 and 0 compare equal. We want
     * -0 to be sorted before 0. We binary search to find the interval containing
     * all zeros and perform a counting sort on it.
     */

    if v.len() <= 1 {
        return;
    }

    // First we move all NaNs to the end
    let mut rnan = v.len() - 1;
    // Skip NaNs already in place
    while rnan > 0 && v[rnan].is_nan() {
        rnan -= 1;
    }
    let mut p = rnan;
    while p > 0 {
        p -= 1;
        if v[p].is_nan() {
            v.swap(p, rnan);
            rnan -= 1;
        }
    }

    // Sort the non-NaN part with efficient comparisons
    sort_by(&mut v[..rnan + 1], &|x: &T, y: &T|
        match x.partial_cmp(y) {
            Some(ord) => ord,
            None      => unsafe { unreachable() }
        });


    let left = find_first_zero(&v[..rnan + 1]);

    // Count zeros of each type and then fill them in in the right order
    let mut zeros = 0;
    let mut neg_zeros = 0;
    for x in v[left..].iter() {
        if *x != zero() {
            break;
        }
        if x.is_sign_negative() {
            neg_zeros += 1;
        } else {
            zeros += 1;
        }
    }
    for x in v[left..].iter_mut() {
        if neg_zeros > 0 {
            *x = Float::neg_zero();
            neg_zeros -= 1;
        } else if zeros > 0 {
             *x = zero();
             zeros -= 1;
        } else {
            break;
        }
    }
}

/// Find the first zero in `v`.
/// If there is no zero, it return v.len() - 1
fn find_first_zero<T: Float>(v: &[T]) -> usize {
    if v.len() == 0 { return 0; }
    let mut hi = v.len() - 1;
    let mut left = 0;

    while left < hi {
        let mid = ((hi - left) / 2) + left;
        if v[mid] < zero() {
            left = mid + 1;
        } else {
            hi = mid;
        }
    }

    while left < v.len() && v[left] < zero() {
        left += 1;
    }
    return left;
}