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

use traits::FloatEFT;
use twosum::fasttwosum;
use split::{split, safesplit_straight};
#[cfg(any(feature = "use-fma", feature = "doc"))]
use fma::{fma, Fma};

#[inline]
pub fn twoproduct<T: FloatEFT>(a: T, b: T) -> (T, T) {
    let prod = a.clone() * b.clone();
    let ((a1, a2), (b1, b2)) = (split(a), split(b));
    (prod.clone(),
     a2.clone() * b2.clone() - (((prod - a1.clone() * b1.clone()) - a1 * b2) - a2 * b1))
}

#[inline]
pub fn safetwoproduct_branch<T: FloatEFT>(a: T, b: T) -> (T, T) {
    let prod = a.clone() * b.clone();
    let ((a1, a2), (b1, b2)) = if a.clone().abs() >=
                                  T::one() / (T::min_positive() / T::epsilon()) {
        (split(a * T::epsilon()), split(b / T::epsilon()))
    } else if b.clone().abs() >= (T::min_positive() / T::epsilon()) {
        (split(a / T::epsilon()), split(b * T::epsilon()))
    } else {
        (split(a), split(b))
    };
    let tmp = if prod.clone().abs() > T::radix() / T::min_positive() {
        ((prod.clone() / T::radix()) - (a1.clone() / T::radix()) * b1.clone()) * T::radix()
    } else {
        prod.clone() - a1.clone() * b1.clone()
    };
    (prod, a2.clone() * b2.clone() - ((tmp - a1 * b2) - a2 * b1))
}

#[inline]
pub fn safetwoproduct_straight<T: FloatEFT>(a: T, b: T) -> (T, T) {
    let prod = a.clone() * b.clone();
    let ((a1, a2, a3), (b1, b2, b3)) = (safesplit_straight(a.clone()),
                                        safesplit_straight(b.clone()));
    let two_a1b1 = T::radix() * (a1.clone() * b1.clone());
    let mid = fasttwosum(prod.clone(), -two_a1b1.clone());
    (prod,
     ((T::radix() * T::radix() * a2.clone()) * b2.clone() -
      ((((mid.0 - two_a1b1) + mid.1) - (T::radix() * T::radix() * b2) * a1) -
       (T::radix() * T::radix() * a2) * b1)) + a * b3 + b * a3)
}

#[cfg(any(feature = "use-fma", feature = "doc"))]
#[inline]
pub fn safetwoproduct_fma<T: FloatEFT + Fma>(a: T, b: T) -> (T, T) {
    let prod = a.clone() * b.clone();
    (prod.clone(), fma(a, b, -prod))
}

#[cfg(test)]
mod tests {
    extern crate rand;

    use std::f64;
    use self::rand::Rng;

    use super::*;

    #[test]
    fn normal() {
        let mut rng = rand::thread_rng();
        for _ in 0..10000000 {
            let (l, r) = (rng.gen_range::<f64>(2f64.powi(-510), 2f64.powi(510)) *
                          rng.choose(&[1., -1.]).unwrap(),
                          rng.gen_range::<f64>(2f64.powi(-510), 2f64.powi(510)) *
                          rng.choose(&[1., -1.]).unwrap());
            let (a1, b1) = safetwoproduct_branch(l, r);
            let (a2, b2) = safetwoproduct_straight(l, r);
            assert!((a1 == a2) && (b1 == b2));

            #[cfg(feature = "use-fma")]
            {
                let (a2, b2) = safetwoproduct_fma(l, r);
                assert!((a1 == a2) && (b1 == b2));
            }
        }
    }

    #[test]
    fn extreme() {
        let mut rng = rand::thread_rng();
        for _ in 0..10000000 {
            let (l, r) = (rng.gen_range::<f64>(2f64.powi(510), f64::MAX) *
                          rng.choose(&[1., -1.]).unwrap(),
                          rng.gen_range::<f64>(f64::MIN_POSITIVE, 2f64.powi(-510)) *
                          rng.choose(&[1., -1.]).unwrap());
            let (a1, b1) = safetwoproduct_branch(l, r);
            let (a2, b2) = safetwoproduct_straight(l, r);
            assert!((a1 == a2) && (b1 == b2));

            #[cfg(feature = "use-fma")]
            {
                let (a2, b2) = safetwoproduct_fma(l, r);
                assert!((a1 == a2) && (b1 == b2));
            }
        }
    }

    #[test]
    fn corner_case() {
        let (l, r) = (6.929001713869936e+236, 2.5944475251952003e+71);
        let (a1, b1) = safetwoproduct_branch(l, r);
        let (a2, b2) = safetwoproduct_straight(l, r);
        assert!((a1 == a2) && (b1 == b2));

        #[cfg(feature = "use-fma")]
        {
            let (a2, b2) = safetwoproduct_fma(l, r);
            assert!((a1 == a2) && (b1 == b2));
        }
    }

    #[test]
    fn subnormal_without_underflow() {
        let mut rng = rand::thread_rng();
        for _ in 0..10000000 {
            let (l, r) = ((rng.gen_range::<i64>(-0x1F_FFFF_FFFF_FFFF, 0x20_0000_0000_0000) as
                           f64) * 2f64.powi(-1022) * 2f64.powi(-52),
                          rng.gen_range::<f64>(2f64.powi(53), f64::MAX) *
                          rng.choose(&[1., -1.]).unwrap());
            let (a1, b1) = safetwoproduct_branch(l, r);
            let (a2, b2) = safetwoproduct_straight(l, r);
            assert!((a1 == a2) && (b1 == b2));

            #[cfg(feature = "use-fma")]
            {
                let (a2, b2) = safetwoproduct_fma(l, r);
                assert!((a1 == a2) && (b1 == b2));
            }
        }
    }
}