gwseq-io 0.2.0

//! Formatting for the text exporters.
//!
//! Two functions, and between them they answer one question: what is the
//! shortest decimal string that reads back as this float?
//!
//! Rust answers the hard half for free. `{}` and `{:e}` both give the shortest
//! digits that round-trip — provably shortest, not "shortest of the precisions
//! we tried" — so nothing here computes a digit. What is left is the choice of
//! *form*, and Rust has no opinion: `{}` never uses an exponent, so `1e30`
//! comes out as thirty-one characters of zeroes, and `{:e}` always does, so
//! `1` comes out as `1e0`.
//!
//! The rule here is C's, minus its digit-counting: **an exponent when the
//! decimal exponent is below -4 or at least 15, and a plain decimal
//! otherwise.** That is the convention every reader of a bedGraph has seen for
//! thirty years — `0.0001` plain and `1e-5` not — and it keeps whole numbers
//! looking like whole numbers, which "whichever is shorter" would not: a
//! coverage of 10000 is not improved by being written `1e4`.
//!
//! The upper cut is 15 rather than the precision, because past it a plain
//! decimal is mostly zeroes carrying no information. Below it, a long plain
//! number is still a number a person can read.
//!
//! Not C's `%.*g`, which walks fixed precisions from `digits10` upward and
//! stops at the first that round-trips — spending six digits on `0.3` where
//! one will do.

use std::fmt::Write as _;

/// Below this decimal exponent, and at or above [`EXPONENT_HIGH`], the
/// exponential form is used. C's `%g` uses -4 for the same reason: `0.0001`
/// reads as a number and `0.00001` reads as a typo.
const EXPONENT_LOW: i32 = -4;

/// Past this, a plain decimal is mostly zeroes.
const EXPONENT_HIGH: i32 = 15;

/// Append `value` as the shortest string that reads back as the same `f32`.
pub fn push_float(out: &mut String, value: f32) {
    push(
        out,
        value as f64,
        |text| text.parse::<f32>() == Ok(value),
        {
            let mut s = String::new();
            let _ = write!(s, "{value}");
            s
        },
        {
            let mut s = String::new();
            let _ = write!(s, "{value:e}");
            s
        },
    );
}

/// The same for an `f64`.
///
/// Separate because the shortest digits differ by width: `0.1f32` and `0.1f64`
/// are different numbers and need different strings, and formatting a double
/// through the single-precision path would round it first.
pub fn push_double(out: &mut String, value: f64) {
    push(
        out,
        value,
        |text| text.parse::<f64>() == Ok(value),
        {
            let mut s = String::new();
            let _ = write!(s, "{value}");
            s
        },
        {
            let mut s = String::new();
            let _ = write!(s, "{value:e}");
            s
        },
    );
}

/// Pick the form and write it.
///
/// `plain` and `exponential` are the two shortest-round-trip renderings; this
/// only decides which. The `round_trips` check is a belt-and-braces assertion
/// in debug builds — Rust guarantees both forms read back, and a change to
/// that guarantee should fail loudly here rather than quietly in a file.
fn push(
    out: &mut String,
    value: f64,
    round_trips: impl Fn(&str) -> bool,
    plain: String,
    exponential: String,
) {
    if value.is_nan() {
        out.push_str("nan");
        return;
    }
    if value.is_infinite() {
        out.push_str(if value < 0.0 { "-inf" } else { "inf" });
        return;
    }
    // The decimal exponent, read off the form that spells it out. `{:e}` is
    // always `<mantissa>e<exp>`, so the tail after the `e` is the exponent and
    // no logarithm is involved — `log10` would round the wrong way at exact
    // powers of ten, which are precisely the boundary cases.
    let exponent: i32 = exponential
        .rsplit_once('e')
        .and_then(|(_, e)| e.parse().ok())
        .unwrap_or(0);
    let chosen = if value == 0.0 || (EXPONENT_LOW..EXPONENT_HIGH).contains(&exponent) {
        &plain
    } else {
        &exponential
    };
    debug_assert!(
        round_trips(chosen),
        "{chosen} does not read back as {value:e}"
    );
    out.push_str(chosen);
}

#[cfg(test)]
mod tests {
    use super::*;

    fn g(value: f32) -> String {
        let mut s = String::new();
        push_float(&mut s, value);
        s
    }

    fn d(value: f64) -> String {
        let mut s = String::new();
        push_double(&mut s, value);
        s
    }

    #[test]
    fn ordinary_values_read_the_way_they_were_written() {
        for (value, expected) in [
            (1.0f32, "1"),
            (0.0, "0"),
            (-2.0, "-2"),
            (12345.0, "12345"),
            (0.1, "0.1"),
            (0.5, "0.5"),
            (12345.678, "12345.678"),
            (999_999.0, "999999"),
        ] {
            assert_eq!(g(value), expected, "for {value:e}");
        }
    }

    /// The whole reason there are two forms: at these magnitudes the plain one
    /// is a wall of zeroes and the exponential one is three characters.
    #[test]
    fn extreme_magnitudes_use_the_exponent() {
        for (value, expected) in [
            (3e-7f32, "3e-7"),
            (1.5e-10, "1.5e-10"),
            (-2.5e-8, "-2.5e-8"),
            (1e20, "1e20"),
            (1e-5, "1e-5"),
        ] {
            assert_eq!(g(value), expected, "for {value:e}");
        }
    }

    /// The form follows the exponent, not the character count. A whole number
    /// stays a whole number even when `1e4` would be shorter.
    #[test]
    fn the_form_follows_the_exponent() {
        assert_eq!(g(10_000.0), "10000"); // exp 4 — inside the plain range
        assert_eq!(g(1e7), "10000000"); // exp 7 — still plain, still readable
        assert_eq!(g(1e-4), "0.0001"); // exp -4 — the boundary, plain
        assert_eq!(g(0.001), "0.001"); // exp -3 — plain
        assert_eq!(g(1e-5), "1e-5"); // exp -5 — past it
        assert_eq!(g(1e15), "1e15"); // exp 15 — past the upper cut
        assert_eq!(g(1e14), "100000000000000"); // exp 14 — just inside
        assert_eq!(g(1.0), "1");
        assert_eq!(g(0.0), "0");
        assert_eq!(g(-0.0), "-0");
    }

    /// Shorter than `%g` where `%g` had to spend a whole precision step, which
    /// is most non-terminating decimals.
    #[test]
    fn shortest_means_shortest_not_a_precision_that_happened_to_work() {
        // `%.6g` of 0.3f32 is "0.3", but of 0.1f32 + 0.2f32 it is "0.300000".
        assert_eq!(g(0.3), "0.3");
        assert_eq!(g(1.0 / 3.0), "0.33333334");
        // An f64 third needs its own digits, not an f32's.
        assert_eq!(d(1.0f64 / 3.0), "0.3333333333333333");
    }

    #[test]
    fn every_value_round_trips() {
        let mut state = 0x9E37_79B9u32;
        for _ in 0..50_000 {
            state ^= state << 13;
            state ^= state >> 17;
            state ^= state << 5;
            let value = f32::from_bits(state);
            if !value.is_finite() {
                continue;
            }
            let text = g(value);
            assert_eq!(
                text.parse::<f32>().unwrap().to_bits(),
                value.to_bits(),
                "{value:e} formatted as {text}"
            );
        }
    }

    #[test]
    fn every_double_round_trips() {
        let mut state = 0x9E37_79B9_7F4A_7C15u64;
        for _ in 0..50_000 {
            state ^= state << 13;
            state ^= state >> 7;
            state ^= state << 17;
            let value = f64::from_bits(state);
            if !value.is_finite() {
                continue;
            }
            let text = d(value);
            assert_eq!(
                text.parse::<f64>().unwrap().to_bits(),
                value.to_bits(),
                "{value:e} formatted as {text}"
            );
        }
    }

    #[test]
    fn non_finite_values_use_the_c_spellings() {
        assert_eq!(g(f32::NAN), "nan");
        assert_eq!(g(f32::INFINITY), "inf");
        assert_eq!(g(f32::NEG_INFINITY), "-inf");
        assert_eq!(d(f64::NAN), "nan");
    }
}