Add some tests for pow

These probably aren't comprehensive but they cover all the edge cases identified in the original musl source.

Author: rrbutani

libm/src/math/pow.rs

@@ -408,3 +408,192 @@ pub fn pow(x: f64, y: f64) -> f64 {
 
     return s * z;
 }
+
+/// Special cases:
+
+///      20. (anything) ** 1 is (anything)
+///      21. (anything) ** -1 is 1/(anything)
+///      22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+///      23. (-anything except 0 and inf) ** (non-integer) is NAN
+
+#[cfg(test)]
+mod tests {
+    // #[macro_use]
+    extern crate std;
+
+    use self::std::f64::consts::{E, PI};
+    use self::std::f64::{EPSILON, INFINITY, MAX, MIN, MIN_POSITIVE, NAN, NEG_INFINITY};
+    use super::pow;
+
+    // const TESTCASES: &[f64] = &[1.0, 0.0, PI, -PI, E, -E, MIN, MAX, MIN_POSITIVE, NAN, INFINITY, NEG_INFINITY];
+
+    const POS_ZERO: &[f64] = &[0.0];
+    const NEG_ZERO: &[f64] = &[-0.0];
+    const POS_ONE: &[f64] = &[1.0];
+    const NEG_ONE: &[f64] = &[-1.0];
+    const POS_FLOATS: &[f64] = &[E, PI, MAX];
+    const NEG_FLOATS: &[f64] = &[-E, -PI, MIN];
+    const POS_SMALL_FLOATS: &[f64] = &[(1.0 / 2.0), MIN_POSITIVE, EPSILON];
+    const NEG_SMALL_FLOATS: &[f64] = &[-(1.0 / 2.0), -MIN_POSITIVE, -EPSILON];
+    const POS_EVENS: &[f64] = &[2.0, 6.0, 8.0, 10.0, 22.0, 100.0];
+    const NEG_EVENS: &[f64] = &[-8.0, -2.0];
+    const POS_ODDS: &[f64] = &[3.0, 7.0];
+    const NEG_ODDS: &[f64] = &[-7.0, -3.0];
+    const NANS: &[f64] = &[NAN];
+    // const EDGES: &[f64] = &[MIN, MAX, MIN_POSITIVE, EPSILON];
+    const POS_INF: &[f64] = &[INFINITY];
+    const NEG_INF: &[f64] = &[NEG_INFINITY];
+
+    const ALL: &[&[f64]] = &[
+        POS_ZERO, NEG_ZERO, NANS, NEG_SMALL_FLOATS, POS_SMALL_FLOATS, NEG_FLOATS, POS_FLOATS, NEG_EVENS, POS_EVENS, NEG_ODDS, POS_ODDS,
+        NEG_INF, POS_INF, NEG_ONE, POS_ONE,
+    ];
+    const POS: &[&[f64]] = &[POS_ZERO, POS_ODDS, POS_ONE, POS_FLOATS, POS_EVENS, POS_INF];
+    const NEG: &[&[f64]] = &[NEG_ZERO, NEG_ODDS, NEG_ONE, NEG_FLOATS, NEG_EVENS, NEG_INF];
+
+    fn pow_test(base: f64, exponent: f64, expected: f64) {
+        let res = pow(base, exponent);
+        assert!(if expected.is_nan() {res.is_nan()} else {pow(base, exponent) == expected},
+                "{} ** {} was {} instead of {}", base, exponent, res, expected);
+    }
+
+    fn test_sets_as_base(sets: &[&[f64]], exponent: f64, expected: f64) {
+        sets.iter()
+            .for_each(|s| s.iter().for_each(|val| pow_test(*val, exponent, expected)));
+    }
+
+    fn test_sets_as_exponent(base: f64, sets: &[&[f64]], expected: f64) {
+        sets.iter()
+            .for_each(|s| s.iter().for_each(|val| pow_test(base, *val, expected)));
+    }
+
+    fn test_sets(sets: &[&[f64]], computed: &Fn(f64) -> f64, expected: &Fn(f64) -> f64) {
+        sets.iter()
+            .for_each(|s| s.iter().for_each(|val| {
+                let exp = expected(*val);
+                let res = computed(*val);
+
+                assert!(if exp.is_nan() {res.is_nan()} else {exp == res},
+                        "test for {} was {} instead of {}", val, res, exp);
+            }));
+    }
+
+    ///      1.  (anything) ** 0  is 1
+    #[test]
+    fn zero_as_exponent() {
+        test_sets_as_base(ALL, 0.0, 1.0);
+        test_sets_as_base(ALL, -0.0, 1.0);
+    }
+
+    ///      2.  1 ** (anything)  is 1
+    #[test]
+    fn one_as_base() {
+        test_sets_as_exponent(1.0, ALL, 1.0);
+    }
+
+    ///      3.  (anything except 1) ** NAN is NAN
+    ///      4.  NAN ** (anything except 0) is NAN
+    #[test]
+    fn nan_inputs() {
+        // NAN as the base:
+        // (NAN ^ anything *but 0* should be NAN)
+        test_sets_as_exponent(NAN, &ALL[2..], NAN);
+
+        // NAN as the exponent:
+        // (anything *but 1* ^ NAN should be NAN)
+        test_sets_as_base(&ALL[..(ALL.len() - 2)], NAN, NAN);
+    }
+
+    ///      16. +INF ** (+anything except 0,NAN) is +INF
+    ///      17. +INF ** (-anything except 0,NAN) is +0
+    ///      18. -INF ** (+odd integer) is -INF
+    ///      19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer)
+    #[test]
+    fn infinity_as_base() {
+        // Positive Infinity as the base:
+        // (+Infinity ^ positive anything but 0 and NAN should be +Infinity)
+        test_sets_as_exponent(INFINITY, &POS[1..], INFINITY);
+
+        // (+Infinity ^ negative anything except 0 and NAN should be 0.0)
+        test_sets_as_exponent(INFINITY, &NEG[1..], 0.0);
+
+        // Negative Infinity as the base:
+        // (-Infinity ^ positive odd ints should be -Infinity)
+        test_sets_as_exponent(NEG_INFINITY, &[POS_ODDS], NEG_INFINITY);
+
+        // (-Infinity ^ anything but odd ints should be == -0 ^ (-anything))
+        // We can lump in pos/neg odd ints here because they don't seem to
+        // cause panics (div by zero) in release mode (I think).
+        test_sets(ALL, &|v: f64| pow(NEG_INFINITY, v), &|v: f64| pow(-0.0, -v));
+    }
+
+    ///      5.  +-(|x| > 1) **  +INF is +INF
+    ///      6.  +-(|x| > 1) **  -INF is +0
+    ///      7.  +-(|x| < 1) **  +INF is +0
+    ///      8.  +-(|x| < 1) **  -INF is +INF
+    ///      9.  -1          ** +-INF is 1
+    #[test]
+    fn infinity_as_exponent() {
+        // Positive/Negative base greater than 1:
+        // (pos/neg > 1 ^ Infinity should be Infinity - note this excludes NAN as the base)
+        test_sets_as_base(&ALL[5..(ALL.len() - 2)], INFINITY, INFINITY);
+
+        // (pos/neg > 1 ^ -Infinity should be 0.0)
+        test_sets_as_base(&ALL[5..(ALL.len() - 2)], NEG_INFINITY, 0.0);
+
+        // Positive/Negative base less than 1:
+        let base_below_one = &[POS_ZERO, NEG_ZERO, NEG_SMALL_FLOATS, POS_SMALL_FLOATS];
+
+        // (pos/neg < 1 ^ Infinity should be 0.0 - this also excludes NAN as the base)
+        test_sets_as_base(base_below_one, INFINITY, 0.0);
+
+        // (pos/neg < 1 ^ -Infinity should be Infinity)
+        test_sets_as_base(base_below_one, NEG_INFINITY, INFINITY);
+
+        // Positive/Negative 1 as the base:
+        // (pos/neg 1 ^ Infinity should be 1)
+        test_sets_as_base(&[NEG_ONE, POS_ONE], INFINITY, 1.0);
+
+        // (pos/neg 1 ^ -Infinity should be 1)
+        test_sets_as_base(&[NEG_ONE, POS_ONE], NEG_INFINITY, 1.0);
+    }
+
+    ///      10. +0 ** (+anything except 0, NAN)               is +0
+    ///      11. -0 ** (+anything except 0, NAN, odd integer)  is +0
+    ///      12. +0 ** (-anything except 0, NAN)               is +INF, raise divbyzero
+    ///      13. -0 ** (-anything except 0, NAN, odd integer)  is +INF, raise divbyzero
+    ///      14. -0 ** (+odd integer) is -0
+    ///      15. -0 ** (-odd integer) is -INF, raise divbyzero
+    #[test]
+    fn zero_as_base() {
+        // Positive Zero as the base:
+        // (+0 ^ anything positive but 0 and NAN should be +0)
+        test_sets_as_exponent(0.0, &POS[1..], 0.0);
+
+        // (+0 ^ anything negative but 0 and NAN should be Infinity)
+        // (this should panic because we're dividing by zero but won't because release mode, I think)
+        test_sets_as_exponent(0.0, &NEG[1..], INFINITY);
+
+        // Negative Zero as the base:
+        // (-0 ^ anything positive but 0, NAN, and odd ints should be +0)
+        test_sets_as_exponent(-0.0,  &POS[3..], 0.0);
+
+        // (-0 ^ anything negative but 0, NAN, and odd ints should be Infinity)
+        // (should panic because of divide by zero)
+        test_sets_as_exponent(-0.0, &NEG[3..], INFINITY);
+
+        // (-0 ^ positive odd ints should be -0)
+        test_sets_as_exponent(-0.0, &[POS_ODDS], -0.0);
+
+        // (-0 ^ negative odd ints should be -Infinity)
+        // (should panic because of divide by zero)
+        test_sets_as_exponent(-0.0, &[NEG_ODDS], NEG_INFINITY);
+    }
+
+    #[test]
+    fn normal_cases() {
+        assert_eq!(pow(2.0, 20.0), (1 << 20) as f64);
+        assert_eq!(pow(-1.0, 9.0), -1.0);
+        assert!(pow(-1.0, 2.2).is_nan());
+    }
+}