Safety of nbody power optimization (f**-1.5) == 1/(sqrt(f) * f) [Always within 1 ULP?]

opt.rs

extern crate rand; // 0.7.3

use rand::prelude::*;
use rand::distributions::Uniform;

fn clear(f: f64) -> f64 {
    f.powf(-1.5)
}

fn fast(f: f64) -> f64 {
    let s = f.sqrt();
    1.0 / (s * f)
}

fn cmp(f: f64) {
    let a = clear(f);
    let b = fast(f);
    let u = if a.is_nan() || b.is_nan() {
        if a.is_nan() == b.is_nan() {
            0
        } else {
            127
        }
    } else {
        diff_ulps(a, b)
    };
    println!("Diff for {}: {} and {} -> {} ulps ({})", f, a, b, u, (a - b).abs());
}

fn diff_ulps(a: f64, b: f64) -> u64 {
    (raw_ulps(a) as i64).wrapping_sub(raw_ulps(b) as i64).abs() as u64
}
fn raw_ulps(f: f64) -> u64 {
    assert!(!f.is_nan());
    let b: u64 = f.to_bits();
    const TOP: u64 = 0x8000_0000_0000_0000;
    let b = if (b as i64) < 0 { TOP.wrapping_sub(b) } else { b };
    assert!((b as i64) >= 0);
    b
}



fn main() {
    for &f in &[-1.0, 0.0, 0.5, 1.0, 2.37, 42.8, f64::INFINITY, 1027.65] {
        cmp(f)
    }
    println!("Testing random:");
    let dist = Uniform::new_inclusive(0.0f64, f64::MAX);
    let mut rng = thread_rng();
    const RUNS: u32 = 25;
    for i in 1..=RUNS {
        print!("{}/{}", i, RUNS);
        let f = rng.sample(dist);
        cmp(f);
    }
}