Quelle  ulps.rs

// Copyright 2014-2018 Optimal Computing (NZ) Ltd.
// Licensed under the MIT license.  See LICENSE for details.

#[cfg(feature="num_traits")]
use num_traits::NumCast;

/// A trait for floating point numbers which computes the number of representable
/// values or ULPs (Units of Least Precision) that separate the two given values.
#[cfg(feature="num_traits")]
pub trait Ulps {
    type U: Copy + NumCast;

    /// The number of representable values or ULPs (Units of Least Precision) that
    /// separate `self` and `other`.  The result `U` is an integral value, and will
    /// be zero if `self` and `other` are exactly equal.
    fn ulps(&self, other: &Self) -> <Self as Ulps>::U;

    /// The next representable number above this one
    fn next(&self) -> Self;

    /// The previous representable number below this one
    fn prev(&self) -> Self;
}

#[cfg(not(feature="num_traits"))]
pub trait Ulps {
    type U: Copy;

    /// The number of representable values or ULPs (Units of Least Precision) that
    /// separate `self` and `other`.  The result `U` is an integral value, and will
    /// be zero if `self` and `other` are exactly equal.
    fn ulps(&self, other: &Self) -> <Self as Ulps>::U;

    /// The next representable number above this one
    fn next(&self) -> Self;

    /// The previous representable number below this one
    fn prev(&self) -> Self;
}

impl Ulps for f32 {
    type U = i32;

    fn ulps(&self, other: &f32) -> i32 {

        // IEEE754 defined floating point storage representation to
        // maintain their order when their bit patterns are interpreted as
        // integers.  This is a huge boon to the task at hand, as we can
        // reinterpret them as integers to find out how many ULPs apart any
        // two floats are

        // Setup integer representations of the input
        let ai32: i32 = self.to_bits() as i32;
        let bi32: i32 = other.to_bits() as i32;

        ai32.wrapping_sub(bi32)
    }

    fn next(&self) -> Self {
        if self.is_infinite() && *self > 0.0 {
            *self
        } else if *self == -0.0 && self.is_sign_negative() {
            0.0
        } else {
            let mut u = self.to_bits();
            if *self >= 0.0 {
                u += 1;
            } else {
                u -= 1;
            }
            f32::from_bits(u)
        }
    }

    fn prev(&self) -> Self {
        if self.is_infinite() && *self < 0.0 {
            *self
        } else if *self == 0.0 && self.is_sign_positive() {
            -0.0
        } else {
            let mut u = self.to_bits();
            if *self <= -0.0 {
                u += 1;
            } else {
                u -= 1;
            }
            f32::from_bits(u)
        }
    }
}

#[test]
fn f32_ulps_test1() {
    let x: f32 = 1000000_f32;
    let y: f32 = 1000000.1_f32;
    println!("DIST IS {}",x.ulps(&y));
    assert!(x.ulps(&y) == -2);
}

#[test]
fn f32_ulps_test2() {
    let pzero: f32 = f32::from_bits(0x00000000_u32);
    let nzero: f32 = f32::from_bits(0x80000000_u32);
    println!("DIST IS {}",pzero.ulps(&nzero));
    assert!(pzero.ulps(&nzero) == -2147483648);
}
#[test]
fn f32_ulps_test3() {
    let pinf: f32 = f32::from_bits(0x7f800000_u32);
    let ninf: f32 = f32::from_bits(0xff800000_u32);
    println!("DIST IS {}",pinf.ulps(&ninf));
    assert!(pinf.ulps(&ninf) == -2147483648);
}

#[test]
fn f32_ulps_test4() {
    let x: f32 = f32::from_bits(0x63a7f026_u32);
    let y: f32 = f32::from_bits(0x63a7f023_u32);
    println!("DIST IS {}",x.ulps(&y));
    assert!(x.ulps(&y) == 3);
}

#[test]
fn f32_ulps_test5() {
    let x: f32 = 2.0;
    let ulps: i32 = x.to_bits() as i32;
    let x2: f32 = <f32>::from_bits(ulps as u32);
    assert_eq!(x, x2);
}

#[test]
fn f32_ulps_test6() {
    let negzero: f32 = -0.;
    let zero: f32 = 0.;
    assert_eq!(negzero.next(), zero);
    assert_eq!(zero.prev(), negzero);
    assert!(negzero.prev() < 0.0);
    assert!(zero.next() > 0.0);
}

impl Ulps for f64 {
    type U = i64;

    fn ulps(&self, other: &f64) -> i64 {

        // IEEE754 defined floating point storage representation to
        // maintain their order when their bit patterns are interpreted as
        // integers.  This is a huge boon to the task at hand, as we can
        // reinterpret them as integers to find out how many ULPs apart any
        // two floats are

        // Setup integer representations of the input
        let ai64: i64 = self.to_bits() as i64;
        let bi64: i64 = other.to_bits() as i64;

        ai64.wrapping_sub(bi64)
    }

    fn next(&self) -> Self {
        if self.is_infinite() && *self > 0.0 {
            *self
        } else if *self == -0.0 && self.is_sign_negative() {
            0.0
        } else {
            let mut u = self.to_bits();
            if *self >= 0.0 {
                u += 1;
            } else {
                u -= 1;
            }
            f64::from_bits(u)
        }
    }

    fn prev(&self) -> Self {
        if self.is_infinite() && *self < 0.0 {
            *self
        } else if *self == 0.0 && self.is_sign_positive() {
            -0.0
        } else {
            let mut u = self.to_bits();
            if *self <= -0.0 {
                u += 1;
            } else {
                u -= 1;
            }
            f64::from_bits(u)
        }
    }
}

#[test]
fn f64_ulps_test1() {
    let x: f64 = 1000000_f64;
    let y: f64 = 1000000.00000001_f64;
    println!("DIST IS {}",x.ulps(&y));
    assert!(x.ulps(&y) == -86);
}

#[test]
fn f64_ulps_test2() {
    let pzero: f64 = f64::from_bits(0x0000000000000000_u64);
    let nzero: f64 = f64::from_bits(0x8000000000000000_u64);
    println!("DIST IS {}",pzero.ulps(&nzero));
    assert!(pzero.ulps(&nzero) == -9223372036854775808i64);
}
#[test]
fn f64_ulps_test3() {
    let pinf: f64 = f64::from_bits(0x7f80000000000000_u64);
    let ninf: f64 = f64::from_bits(0xff80000000000000_u64);
    println!("DIST IS {}",pinf.ulps(&ninf));
    assert!(pinf.ulps(&ninf) == -9223372036854775808i64);
}

#[test]
fn f64_ulps_test4() {
    let x: f64 = f64::from_bits(0xd017f6cc63a7f026_u64);
    let y: f64 = f64::from_bits(0xd017f6cc63a7f023_u64);
    println!("DIST IS {}",x.ulps(&y));
    assert!(x.ulps(&y) == 3);
}

#[test]
fn f64_ulps_test5() {
    let x: f64 = 2.0;
    let ulps: i64 = x.to_bits() as i64;
    let x2: f64 = <f64>::from_bits(ulps as u64);
    assert_eq!(x, x2);
}

#[test]
fn f64_ulps_test6() {
    let negzero: f64 = -0.;
    let zero: f64 = 0.;
    assert_eq!(negzero.next(), zero);
    assert_eq!(zero.prev(), negzero);
    assert!(negzero.prev() < 0.0);
    assert!(zero.next() > 0.0);
}