safeeft/
split.rs

1use traits::FloatEFT;
2
3#[inline]
4pub fn split<T: FloatEFT>(a: T) -> (T, T) {
5    let tmp = a.clone() * T::split_coef();
6    let x = tmp.clone() - (tmp - a.clone());
7    (x.clone(), a - x)
8}
9
10#[inline]
11pub fn safesplit_branch<T: FloatEFT>(a: T) -> (T, T) {
12    // unsafe when usp(a) >= 2^997 (a >= 0x1.FFFFFF8000000p+1022) <-- ?????
13    // theoritically split(DOUBLE_MAX) == (succ(DOUBLE_MAX)==2^1023, some negative float)
14    if a > T::one() {
15        let t = split(a * T::epsilon());
16        (t.0 / T::epsilon(), t.1 / T::epsilon())
17    } else {
18        split(a)
19    }
20}
21
22#[inline]
23pub fn safesplit_straight<T: FloatEFT>(a: T) -> (T, T, T) {
24    // Returns a_high, a_low, a_err which satisfy a == 2 * a_high + 2 * a_low + a_err.
25    // 2 * a_high may overflow, so to get a, you should write a_high + (a_high + (2.*a_low + a_err))
26    let aa = a.clone() / T::radix();
27    let err = a - aa.clone() * T::radix(); // if usp(a) == 2^-1074, err == 2^-1074, else 0.
28
29    let step = (((aa.clone() + T::min_positive()) - aa.clone()) / T::min_positive()) /
30               (T::epsilon() * T::epsilon() * T::epsilon()) + T::epsilon();
31    let split_shift = split(aa * step.clone());
32
33    (split_shift.0 / step.clone(), split_shift.1 / step, err)
34}
35
36#[cfg(test)]
37mod tests {
38    extern crate num_traits;
39    extern crate rand;
40
41    use std::f64;
42    use self::rand::Rng;
43    use self::num_traits::Float;
44
45    use super::*;
46
47    #[allow(dead_code)]
48    fn is_split(a: f64, b: f64) -> bool {
49        let (a_mant, a_expo, _) = a.integer_decode();
50        let (b_mant, b_expo, _) = b.integer_decode();
51
52        ((a_expo - b_expo >= 27) && (((a_mant | b_mant) & 0x7FFFFFFu64) == 0)) |
53        ((a_expo - b_expo == 26) && ((a_mant & 0x7FFFFFFu64) == 0) && (a_mant % 2 == 0) &&
54         (b_mant == 0x10_0000_0000_0000)) |
55        ((a_expo == b_expo) && (b_expo == -1022) &&
56         (b_mant <= 2u64.pow(a_mant.trailing_zeros() - 1)))
57    }
58
59    #[test]
60    fn normal() {
61        let mut rng = rand::thread_rng();
62        for _ in 0..10000000 {
63            let fl = rng.gen_range(2f64.powi(-1022 + 27), 1e+308) *
64                     (1. - 2. * (rng.gen_range(0, 1) as f64));
65
66            let s1 = split(fl);
67            if s1.0.is_nan() {
68                continue;
69            }
70            assert!(s1.0 + s1.1 == fl);
71
72            let s2 = safesplit_straight(fl);
73
74            let (s2h_prop, s2l_prop) = (s2.0.integer_decode(), s2.1.integer_decode());
75            assert!((s2h_prop.0 & 0x7FFFFFF) == 0);
76            assert!((s2l_prop.0 & 0x7FFFFFF) == 0);
77            assert!(s2.2.abs() <= f64::MIN_POSITIVE * 2f64.powi(-52));
78            assert!(s2.0.abs() * 2f64.powi(-26) >= s2.1.abs());
79
80            assert_eq!(s1.0, s2.0 * 2.);
81            assert_eq!(s1.1, s2.1 * 2.);
82            //assert_eq!(s1.0 + s1.1, s2.0 * 2. + s2.1 * 2. + s2.2);
83        }
84    }
85
86    #[test]
87    fn large() {
88        let mut rng = rand::thread_rng();
89        for _ in 0..10000000 {
90            let fl = (rng.gen_range::<i64>(-0x1F_FFFF_FFFF_FFFF, 0x20_0000_0000_0000) as f64) *
91                     2f64.powi(1022 - 52);
92
93            let s = safesplit_straight(fl);
94
95            if s.0.is_nan() {
96                unreachable!()
97            }
98
99            let (sh_prop, sl_prop) = (s.0.integer_decode(), s.1.integer_decode());
100            println!("{:e}, {:b}, {:b}", s.0, sh_prop.0, sh_prop.1);
101            assert!((sh_prop.0 & 0x7FFFFFF) == 0);
102            assert!((sl_prop.0 & 0x7FFFFFF) == 0);
103
104            assert_eq!(((s.2 + s.1 * 2.) + s.0) + s.0, fl);
105            assert!(s.0.abs() * 2f64.powi(-26) >= s.1.abs());
106        }
107
108        let s = safesplit_straight(f64::MAX);
109
110        let (sh_prop, sl_prop) = (s.0.integer_decode(), s.1.integer_decode());
111        assert!((sh_prop.0 & 0x7FFFFFF) == 0);
112        assert!((sl_prop.0 & 0x7FFFFFF) == 0);
113
114        assert_eq!(((s.2 + s.1 * 2.) + s.0) + s.0, f64::MAX);
115        assert!(s.0.abs() * 2f64.powi(-26) >= s.1.abs());
116    }
117
118    #[test]
119    fn subnormal() {
120        let mut rng = rand::thread_rng();
121        for _ in 0..10000000 {
122            let fl = ((rng.gen_range::<i64>(-0xFFFF_FFFF_FFFF, 0x1_0000_0000_0000) as f64) *
123                      2f64.powi(-1022)) * 2f64.powi(-52);
124            let s = safesplit_straight(fl);
125            assert_eq!(((s.2 + s.1 * 2.) + s.0) + s.0, fl);
126            assert!(s.0.abs() * 2f64.powi(-26) >= s.1.abs());
127        }
128    }
129}
safeeft/split.rs

safeeft/
split.rs
