decimal_scaled/

consts.rs

//! Mathematical constants and float-compatibility constants for every
//! decimal width.
//!
//! # Constants provided
//!
//! The [`DecimalConsts`] trait exposes `pi`, `tau`, `half_pi`,
//! `quarter_pi`, `golden`, and `e` as methods on every width. The
//! native-tier (`D38` and narrower) impls live here; the wide tier
//! (`D76` / `D153` / `D307`) impls live in `consts_wide.rs`.
//!
//! Two inherent associated constants, `EPSILON` and `MIN_POSITIVE`, are
//! provided as analogues to `f64::EPSILON` and `f64::MIN_POSITIVE` so
//! that generic code parameterised over numeric types continues to
//! compile when `T` is any of the decimal widths.
//!
//! # Precision strategy
//!
//! Constants are derived from raw integer references — no `f64`
//! anywhere. Each tier has its own reference at the tier's maximum
//! storage precision; the rescale to the caller's `SCALE` is always
//! **downward**, never upward, so half-to-even rounding always lands
//! on the **correctly-rounded** value at the target scale:
//!
//! | Tier           | Reference storage | `SCALE_REF` (= reference digits) | Source file       |
//! |----------------|-------------------|----------------------------------|-------------------|
//! | D9 / D18 / D38 | `Int256`          | 75                               | this file         |
//! | D76            | `Int256`          | 75                               | `consts_wide.rs`  |
//! | D153           | `Int512`          | 153                              | `consts_wide.rs`  |
//! | D307           | `Int1024`         | 307                              | `consts_wide.rs`  |
//!
//! The rescale from `SCALE_REF` to the caller's `SCALE` uses integer
//! division with the crate-default [`RoundingMode`] (half-to-even by
//! default; overridable via the `rounding-*` Cargo features). Going
//! through `f64` would cap precision at ~15–17 decimal digits; the
//! raw-integer path preserves the full per-tier reference width.
//!
//! **0.5 ULP at every supported scale**, on every width, with no
//! exceptions in the precision contract. The only constraint is the
//! width's *storage range*: a value that mathematically exceeds the
//! type's `Storage::MAX / 10^SCALE` cannot be represented at all. At
//! `D38<38>` the storage range is approximately ±1.7, so the four
//! larger-magnitude constants — `pi ≈ 3.14`, `tau ≈ 6.28`, `e ≈ 2.72`,
//! `golden ≈ 1.62` — overflow `i128` and the corresponding methods
//! panic with a clear "constant out of storage range" message;
//! `half_pi ≈ 1.57` and `quarter_pi ≈ 0.79` fit and remain
//! correctly-rounded to 0.5 ULP.
//!
//! [`RoundingMode`]: crate::rounding::RoundingMode
//!
//! # Sources
//!
//! Each raw constant is the half-to-even rounding of the canonical
//! decimal expansion to the tier's `SCALE_REF` fractional digits. ISO
//! 80000-2 (pi, tau, pi/2, pi/4), OEIS A001113 (e), OEIS A001622
//! (golden ratio).

use crate::core_type::D38;
use crate::d_w128_kernels::Fixed;
use crate::wide_int::Int256;

/// Reference scale for every constant in this file: the 75-digit
/// representation that fits an `Int256` (`2 · 128` bits). Every D38
/// scale (0..=38) is at most 38 digits, so we always rescale **down**
/// from 75 → SCALE, never up. The half-to-even rescale-down step is
/// performed by [`Fixed::round_to_i128`] (`Fixed` is the same 256-bit
/// guard-digit type the strict transcendentals use), giving 0.5 ULP at
/// the caller's `SCALE` for every value that fits `i128` at that
/// scale.
///
/// # Precision
///
/// N/A: constant value, no arithmetic performed.
const SCALE_REF: u32 = 75;

// Raw decimal strings at 75 fractional digits, materialised at build
// time by `build.rs` (the same hand-rolled multi-precision generator
// that emits the wide-tier constants). Sources: ISO 80000-2 (pi, tau,
// pi/2, pi/4), OEIS A001113 (e), OEIS A001622 (golden ratio).
//
// The build-time string -> Int256 parse is `const fn` (via
// `Int256::from_str_radix`, base 10 only). The 75-digit reference is
// the largest decimal expansion that always fits Int256 for the
// biggest of these constants (tau ≈ 6.28×10⁷⁵ < Int256::MAX ≈
// 5.78×10⁷⁶); a single shared SCALE_REF keeps the rescale helpers
// uniform across all six methods on the trait.

include!(concat!(env!("OUT_DIR"), "/wide_consts.rs"));

const PI_RAW: Int256 = match Int256::from_str_radix(PI_D76_S75, 10) {
    Ok(v) => v,
    Err(_) => panic!("consts: PI_D76_S75 not parseable"),
};
const TAU_RAW: Int256 = match Int256::from_str_radix(TAU_D76_S75, 10) {
    Ok(v) => v,
    Err(_) => panic!("consts: TAU_D76_S75 not parseable"),
};
const HALF_PI_RAW: Int256 = match Int256::from_str_radix(HALF_PI_D76_S75, 10) {
    Ok(v) => v,
    Err(_) => panic!("consts: HALF_PI_D76_S75 not parseable"),
};
const QUARTER_PI_RAW: Int256 = match Int256::from_str_radix(QUARTER_PI_D76_S75, 10) {
    Ok(v) => v,
    Err(_) => panic!("consts: QUARTER_PI_D76_S75 not parseable"),
};
const E_RAW: Int256 = match Int256::from_str_radix(E_D76_S75, 10) {
    Ok(v) => v,
    Err(_) => panic!("consts: E_D76_S75 not parseable"),
};
const GOLDEN_RAW: Int256 = match Int256::from_str_radix(GOLDEN_D76_S75, 10) {
    Ok(v) => v,
    Err(_) => panic!("consts: GOLDEN_D76_S75 not parseable"),
};

/// Rescale a 75-digit `Int256` reference down to the caller's `TARGET`
/// scale as an `i128`, half-to-even. Panics if the value at `TARGET`
/// does not fit `i128` (the type's storage range at that scale just
/// doesn't include this constant — e.g. `pi ≈ 3.14` at `D38<38>` would
/// need `3.14 × 10^38 ≈ 3.14e38`, which exceeds `i128::MAX ≈ 1.7e38`).
fn rescale_75_to_target<const TARGET: u32>(raw: Int256, name: &'static str) -> i128 {
    let limbs = raw.0;  // [u128; 2], little-endian
    let f = Fixed { negative: false, mag: limbs };
    match f.round_to_i128(SCALE_REF, TARGET) {
        Some(v) => v,
        None => panic!(
            "D38 constant out of storage range: {name} cannot fit i128 at SCALE = {TARGET} \
             (storage range is ±i128::MAX / 10^SCALE)",
            name = name,
            TARGET = TARGET,
        ),
    }
}

/// Well-known mathematical constants available on every decimal width
/// (`D9` / `D18` / `D38` / `D76` / `D153` / `D307`).
///
/// Import this trait to call `D38s12::pi()`, `D76::<35>::e()`, etc.
///
/// All returned values are computed from a raw integer reference at
/// the tier's maximum storage precision (75 digits for D9/D18/D38 and
/// D76; 153 for D153; 307 for D307) without passing through `f64`,
/// then rescaled down to the caller's `SCALE` with half-to-even
/// rounding. The result is **within 0.5 ULP** of the canonical
/// decimal expansion at every supported scale on every width.
///
/// The one situation where a method does not return a value is when
/// the constant's magnitude exceeds the type's storage range at the
/// caller's `SCALE` — e.g. `D38<38>::pi()` would need `3.14 × 10³⁸`,
/// which exceeds `i128::MAX ≈ 1.7×10³⁸`. The method panics with a
/// clear "constant out of storage range" message in that case.
pub trait DecimalConsts: Sized {
    /// Pi (~3.14159265...). One half-turn in radians.
    ///
    /// Source: ISO 80000-2 / OEIS A000796. Rescaled per-tier (see the
    /// module-level table) to the caller's `SCALE` via the crate-default
    /// rounding mode.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    fn pi() -> Self;

    /// Tau (~6.28318530...). One full turn in radians.
    ///
    /// Defined as `2 * pi`. Rescaled per-tier (see the module-level table) to the caller's `SCALE` via the crate-default rounding mode.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    fn tau() -> Self;

    /// Half-pi (~1.57079632...). One quarter-turn in radians.
    ///
    /// Defined as `pi / 2`. Rescaled per-tier (see the module-level table) to the caller's `SCALE` via the crate-default rounding mode.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    fn half_pi() -> Self;

    /// Quarter-pi (~0.78539816...). One eighth-turn in radians.
    ///
    /// Defined as `pi / 4`. Rescaled per-tier (see the module-level table) to the caller's `SCALE` via the crate-default rounding mode.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    fn quarter_pi() -> Self;

    /// The golden ratio (~1.61803398...). Dimensionless.
    ///
    /// Defined as `(1 + sqrt(5)) / 2`. Source: OEIS A001622. Rescaled
    /// per-tier (see the module-level table) to the caller's `SCALE`
    /// via the crate-default rounding mode.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    fn golden() -> Self;

    /// Euler's number (~2.71828182...). Dimensionless.
    ///
    /// Source: OEIS A001113. Rescaled per-tier (see the module-level table) to the caller's `SCALE` via the crate-default rounding mode.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    fn e() -> Self;
}

// Public-to-crate helpers that return each constant's rescaled bits at
// the caller's target SCALE. Used by the `decl_decimal_consts!` macro
// to provide DecimalConsts for narrower widths (D9, D18) without
// duplicating the rescale logic.

pub(crate) fn pi_at_target<const TARGET: u32>() -> i128 {
    rescale_75_to_target::<TARGET>(PI_RAW, "pi")
}
pub(crate) fn tau_at_target<const TARGET: u32>() -> i128 {
    rescale_75_to_target::<TARGET>(TAU_RAW, "tau")
}
pub(crate) fn half_pi_at_target<const TARGET: u32>() -> i128 {
    rescale_75_to_target::<TARGET>(HALF_PI_RAW, "half_pi")
}
pub(crate) fn quarter_pi_at_target<const TARGET: u32>() -> i128 {
    rescale_75_to_target::<TARGET>(QUARTER_PI_RAW, "quarter_pi")
}
pub(crate) fn golden_at_target<const TARGET: u32>() -> i128 {
    rescale_75_to_target::<TARGET>(GOLDEN_RAW, "golden")
}
pub(crate) fn e_at_target<const TARGET: u32>() -> i128 {
    rescale_75_to_target::<TARGET>(E_RAW, "e")
}

// The `DecimalConsts` impl for `D38<SCALE>` is emitted by the
// `decl_decimal_consts!` macro — the same macro D9 / D18 / D76+ use.
// It expands to `Self(pi_at_target::<SCALE>())` etc.; each
// `*_at_target` helper above rescales the 75-digit Int256 reference
// down to the caller's `SCALE` via half-to-even and narrows to i128
// (or panics with a clear message if the constant's magnitude
// exceeds the storage range at that scale).
crate::macros::consts::decl_decimal_consts!(D38, i128);

// Inherent associated constants: EPSILON / MIN_POSITIVE.
//
// These mirror `f64::EPSILON` and `f64::MIN_POSITIVE` so that generic
// numeric code that calls `T::EPSILON` or `T::MIN_POSITIVE` compiles
// when `T = D38<SCALE>`. For D38 both equal `D38(1)` -- the smallest
// representable positive value (1 LSB = 10^-SCALE). There are no subnormals.

impl<const SCALE: u32> D38<SCALE> {
    /// Smallest representable positive value: 1 LSB = `10^-SCALE`.
    ///
    /// Provided as an analogue to `f64::EPSILON` for generic numeric code.
    /// Note that this differs from the f64 definition ("difference between
    /// 1.0 and the next-larger f64"): for `D38` the LSB is uniform across
    /// the entire representable range.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    pub const EPSILON: Self = Self(1);

    /// Smallest positive value (equal to [`Self::EPSILON`]).
    ///
    /// Provided as an analogue to `f64::MIN_POSITIVE` for generic numeric
    /// code. Unlike `f64`, `D38` has no subnormals, so `MIN_POSITIVE`
    /// and `EPSILON` are the same value.
    ///
    /// # Precision
    ///
    /// N/A: constant value, no arithmetic performed.
    pub const MIN_POSITIVE: Self = Self(1);
}

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

    // Bit-exact assertions at SCALE = 12.
    //
    // At SCALE = 12 each constant is the 37-digit raw integer divided by
    // 10^23, rounded half-to-even.

    /// pi at SCALE=12: raw / 10^23.
    /// Truncated 13 digits: 3_141_592_653_589.
    /// 14th digit is 7 (from position 14 of the raw) -> round up.
    /// Expected: 3_141_592_653_590.
    #[test]
    fn pi_is_bit_exact_at_scale_12() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        assert_eq!(D38s12::pi().to_bits(), 3_141_592_653_590_i128);
    }

    /// tau at SCALE=12: raw / 10^23.
    /// Truncated 13 digits: 6_283_185_307_179.
    /// 14th digit is 5 -> round up. Expected: 6_283_185_307_180.
    #[test]
    fn tau_is_bit_exact_at_scale_12() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        assert_eq!(D38s12::tau().to_bits(), 6_283_185_307_180_i128);
    }

    /// half_pi at SCALE=12: raw / 10^23.
    /// Truncated 13 digits: 1_570_796_326_794.
    /// 14th digit is 8 -> round up. Expected: 1_570_796_326_795.
    #[test]
    fn half_pi_is_bit_exact_at_scale_12() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        assert_eq!(D38s12::half_pi().to_bits(), 1_570_796_326_795_i128);
    }

    /// quarter_pi at SCALE=12: raw / 10^23.
    /// Truncated 12 digits: 785_398_163_397.
    /// 13th digit is 4 -> no round-up. Expected: 785_398_163_397.
    #[test]
    fn quarter_pi_is_bit_exact_at_scale_12() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        assert_eq!(D38s12::quarter_pi().to_bits(), 785_398_163_397_i128);
    }

    /// e at SCALE=12: raw / 10^23.
    /// Truncated 13 digits: 2_718_281_828_459.
    /// 14th digit is 0 -> no round-up. Expected: 2_718_281_828_459.
    #[test]
    fn e_is_bit_exact_at_scale_12() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        assert_eq!(D38s12::e().to_bits(), 2_718_281_828_459_i128);
    }

    /// golden at SCALE=12: raw / 10^23.
    /// Truncated 13 digits: 1_618_033_988_749.
    /// 14th digit is 8 -> round up. Expected: 1_618_033_988_750.
    #[test]
    fn golden_is_bit_exact_at_scale_12() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        assert_eq!(D38s12::golden().to_bits(), 1_618_033_988_750_i128);
    }

    // Closeness checks against core::f64::consts.
    // These verify that the correct reference digits were selected; the
    // bit-exact tests above are the primary acceptance criteria.

    /// pi() converted to f64 is within 1e-11 of `core::f64::consts::PI`.
    /// At SCALE=12, 1 LSB = 1e-12, so 1e-11 covers rescale rounding plus
    /// the f64 conversion step.
    #[test]
    fn pi_close_to_f64_pi() {
        let diff = (D38s12::pi().to_f64() - core::f64::consts::PI).abs();
        assert!(diff < 1e-11, "pi diverges from f64 PI by {diff}");
    }

    #[test]
    fn tau_close_to_f64_tau() {
        let diff = (D38s12::tau().to_f64() - core::f64::consts::TAU).abs();
        assert!(diff < 1e-11, "tau diverges from f64 TAU by {diff}");
    }

    #[test]
    fn half_pi_close_to_f64_frac_pi_2() {
        let diff =
            (D38s12::half_pi().to_f64() - core::f64::consts::FRAC_PI_2).abs();
        assert!(diff < 1e-11, "half_pi diverges from f64 FRAC_PI_2 by {diff}");
    }

    #[test]
    fn quarter_pi_close_to_f64_frac_pi_4() {
        let diff =
            (D38s12::quarter_pi().to_f64() - core::f64::consts::FRAC_PI_4).abs();
        assert!(
            diff < 1e-11,
            "quarter_pi diverges from f64 FRAC_PI_4 by {diff}"
        );
    }

    #[test]
    fn e_close_to_f64_e() {
        let diff = (D38s12::e().to_f64() - core::f64::consts::E).abs();
        assert!(diff < 1e-11, "e diverges from f64 E by {diff}");
    }

    /// golden() converted to f64 is within 1e-11 of the closed form
    /// `(1 + sqrt(5)) / 2`. Requires std for `f64::sqrt`.
    #[cfg(feature = "std")]
    #[test]
    fn golden_close_to_closed_form() {
        let expected = (1.0_f64 + 5.0_f64.sqrt()) / 2.0;
        let diff = (D38s12::golden().to_f64() - expected).abs();
        assert!(diff < 1e-11, "golden diverges from closed-form by {diff}");
    }

    // EPSILON / MIN_POSITIVE

    #[test]
    fn epsilon_is_one_ulp() {
        assert_eq!(D38s12::EPSILON.to_bits(), 1_i128);
        assert!(D38s12::EPSILON > D38s12::ZERO);
    }

    #[test]
    fn min_positive_is_one_ulp() {
        assert_eq!(D38s12::MIN_POSITIVE.to_bits(), 1_i128);
        assert_eq!(D38s12::MIN_POSITIVE, D38s12::EPSILON);
    }

    /// At SCALE = 6 the LSB is 10^-6; EPSILON is still raw 1.
    #[test]
    fn epsilon_at_scale_6_is_one_ulp() {
        type D6 = D38<6>;
        assert_eq!(D6::EPSILON.to_bits(), 1_i128);
        assert_eq!(D6::MIN_POSITIVE.to_bits(), 1_i128);
    }

    // Cross-scale exercises

    /// At SCALE = 6, pi() should equal 3.141593 (rounded half-to-even from
    /// 3.1415926535...). Expected raw bits: 3_141_593.
    #[test]
    fn pi_at_scale_6_is_bit_exact() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        type D6 = D38<6>;
        assert_eq!(D6::pi().to_bits(), 3_141_593_i128);
    }

    /// At SCALE = 0, pi() rounds to 3 (first fractional digit is 1, no
    /// round-up).
    #[test]
    fn pi_at_scale_0_is_three() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        type D0 = D38<0>;
        assert_eq!(D0::pi().to_bits(), 3_i128);
    }

    /// `D38<37>::pi()` is the canonical pi rounded half-to-even to 37
    /// fractional digits. The 75-digit Int256 reference is rescaled
    /// down to 37 digits; the result is bit-identical to the
    /// hand-tabulated constant.
    #[test]
    fn pi_at_scale_37_matches_canonical_37_digit_rounding() {
        type D37 = D38<37>;
        // pi to 38 digits: 3.14159265358979323846264338327950288420
        //                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        //                   keep 37 frac digits; the 38th digit is 0
        //                   so half-to-even rounds down — no bump.
        let expected: i128 = 31_415_926_535_897_932_384_626_433_832_795_028_842;
        assert_eq!(D37::pi().to_bits(), expected);
    }

    /// `D38<38>` cannot represent pi (storage range is ~±1.7 at
    /// SCALE=38, pi is 3.14). The method panics with a clear
    /// "constant out of storage range" message rather than silently
    /// returning a fudged value.
    #[test]
    #[should_panic(expected = "out of storage range")]
    fn pi_at_scale_38_panics_storage_range() {
        let _ = D38::<38>::pi();
    }

    /// `D38<38>` storage range covers ±1.7, which DOES include
    /// `half_pi ≈ 1.57` and `quarter_pi ≈ 0.79`. They must be
    /// correctly rounded to 0.5 ULP (= 1 LSB).
    #[test]
    fn half_pi_and_quarter_pi_at_scale_38_are_correctly_rounded() {
        // half_pi to 38 digits: 1.57079632679489661923132169163975144210
        let expected_half_pi: i128 = 157_079_632_679_489_661_923_132_169_163_975_144_210;
        let got = D38::<38>::half_pi().to_bits();
        let diff = (got - expected_half_pi).abs();
        assert!(diff <= 1, "half_pi: got {got}, expected {expected_half_pi}, diff {diff} > 1 LSB");

        // quarter_pi to 38 digits: 0.78539816339744830961566084581987572105
        let expected_quarter_pi: i128 = 78_539_816_339_744_830_961_566_084_581_987_572_105;
        let got = D38::<38>::quarter_pi().to_bits();
        let diff = (got - expected_quarter_pi).abs();
        assert!(diff <= 1, "quarter_pi: got {got}, expected {expected_quarter_pi}, diff {diff} > 1 LSB");
    }

    /// Negative-side rounding: negating pi gives the expected raw bits.
    #[test]
    fn neg_pi_round_trip() {
        if !crate::rounding::DEFAULT_IS_HALF_TO_EVEN { return; }
        let pi = D38s12::pi();
        let neg_pi = -pi;
        assert_eq!(neg_pi.to_bits(), -3_141_592_653_590_i128);
    }

    // (`rescale_from_ref` boundary tests removed: the rounding logic now
    // lives in `D38::rescale` / `src/rounding.rs::apply_rounding` and is
    // covered by the tests in those modules.)
}