#[cfg(zmij_no_select_unpredictable)] mod hint; #[cfg(all(target_arch = "x86_64", target_feature = "sse2", not(miri)))] mod stdarch_x86; #[cfg(test)] mod tests; mod traits;
#[cfg(all(any(target_arch = "aarch64", target_arch = "x86_64"), not(miri)))] use core::arch::asm; #[cfg(not(zmij_no_select_unpredictable))] use core::hint; use core::mem::{self, MaybeUninit}; use core::ptr; use core::slice; use core::str; #[cfg(feature = "no-panic")] use no_panic::no_panic;
#[cfg_attr(feature = "no-panic", no_panic)] fn umul192_hi128(x_hi: u64, x_lo: u64, y: u64) -> uint128 { let p = umul128(x_hi, y); let lo = (p as u64).wrapping_add((umul128(x_lo, y) >> 64) as u64);
uint128 {
hi: (p >> 64) as u64 + u64::from(lo < p as u64),
lo,
}
}
// Computes high 64 bits of multiplication of x and y, discards the least // significant bit and rounds to odd, where x = uint128_t(x_hi << 64) | x_lo. #[cfg_attr(feature = "no-panic", no_panic)] fn umulhi_inexact_to_odd<UInt>(x_hi: u64, x_lo: u64, y: UInt) -> UInt where
UInt: traits::UInt,
{ let num_bits = mem::size_of::<UInt>() * 8; if num_bits == 64 { let p = umul192_hi128(x_hi, x_lo, y.into());
UInt::truncate(p.hi | u64::from((p.lo >> 1) != 0))
} else { let p = (umul128(x_hi, y.into()) >> 32) as u64;
UInt::enlarge((p >> 32) as u32 | u32::from((p as u32 >> 1) != 0))
}
}
// 128-bit significands of powers of 10 rounded down. // Generation with 192-bit arithmetic and compression by Dougall Johnson. static POW10_SIGNIFICANDS: Pow10SignificandsTable = { letmut data = [0; if Pow10SignificandsTable::COMPRESS { 0
} else {
Pow10SignificandsTable::NUM_POW10 * 2
}];
struct uint192 {
w0: u64, // least significant
w1: u64,
w2: u64, // most significant
}
// First element, rounded up to cancel out rounding down in the // multiplication, and minimize significant bits. letmut current = uint192 {
w0: 0xe000000000000000,
w1: 0x25e8e89c13bb0f7a,
w2: 0xff77b1fcbebcdc4f,
}; let ten = 0xa000000000000000; letmut i = 0; while i < Pow10SignificandsTable::NUM_POW10 && !Pow10SignificandsTable::COMPRESS { if Pow10SignificandsTable::SPLIT_TABLES {
data[Pow10SignificandsTable::NUM_POW10 - i - 1] = current.w2;
data[Pow10SignificandsTable::NUM_POW10 * 2 - i - 1] = current.w1;
} else {
data[i * 2] = current.w2;
data[i * 2 + 1] = current.w1;
}
let h0: u64 = umul128_hi64(current.w0, ten); let h1: u64 = umul128_hi64(current.w1, ten);
let c0: u64 = h0.wrapping_add(current.w1.wrapping_mul(ten)); let c1: u64 = ((c0 < h0) as u64 + h1).wrapping_add(current.w2.wrapping_mul(ten)); let c2: u64 = (c1 < h1) as u64 + umul128_hi64(current.w2, ten); // dodgy carry
static EXP_SHIFTS: ExpShiftTable = { letmut data = [0u8; if ExpShiftTable::ENABLE {
f64::EXP_MASK as usize + 1
} else { 1
}];
letmut raw_exp = 0; while raw_exp < data.len() && ExpShiftTable::ENABLE { letmut bin_exp = raw_exp as i32 - f64::EXP_OFFSET; if raw_exp == 0 {
bin_exp += 1;
} let dec_exp = compute_dec_exp(bin_exp, true);
data[raw_exp] = do_compute_exp_shift(bin_exp, dec_exp) as u8;
raw_exp += 1;
}
ExpShiftTable { data }
};
// Computes a shift so that, after scaling by a power of 10, the intermediate // result always has a fixed 128-bit fractional part (for double). // // Different binary exponents can map to the same decimal exponent, but place // the decimal point at different bit positions. The shift compensates for this. // // For example, both 3 * 2**59 and 3 * 2**60 have dec_exp = 2, but dividing by // 10^dec_exp puts the decimal point in different bit positions: // 3 * 2**59 / 100 = 1.72...e+16 (needs shift = 1 + 1) // 3 * 2**60 / 100 = 3.45...e+16 (needs shift = 2 + 1) #[inline] unsafefn compute_exp_shift<UInt, const ONLY_REGULAR: bool>(bin_exp: i32, dec_exp: i32) -> u8 where
UInt: traits::UInt,
{ let num_bits = mem::size_of::<UInt>() * 8; if num_bits == 64 && ExpShiftTable::ENABLE && ONLY_REGULAR { unsafe {
*EXP_SHIFTS
.data
.as_ptr()
.add((bin_exp + f64::EXP_OFFSET) as usize)
}
} else {
do_compute_exp_shift(bin_exp, dec_exp)
}
}
#[cfg_attr(feature = "no-panic", no_panic)] fn count_trailing_nonzeros(x: u64) -> usize { // We count the number of bytes until there are only zeros left. // The code is equivalent to // 8 - x.leading_zeros() / 8 // but if the BSR instruction is emitted (as gcc on x64 does with default // settings), subtracting the constant before dividing allows the compiler // to combine it with the subtraction which it inserts due to BSR counting // in the opposite direction. // // Additionally, the BSR instruction requires a zero check. Since the high // bit is unused we can avoid the zero check by shifting the datum left by // one and inserting a sentinel bit at the end. This can be faster than the // automatically inserted range check.
(70 - ((x.to_le() << 1) | 1).leading_zeros() as usize) / 8
}
// Align data since unaligned access may be slower when crossing a // hardware-specific boundary. #[repr(C, align(2))] struct Digits2([u8; 200]);
// Converts value in the range [0, 100) to a string. GCC generates a bit better // code when value is pointer-size (https://www.godbolt.org/z/5fEPMT1cc). #[cfg_attr(feature = "no-panic", no_panic)] unsafefn digits2(value: usize) -> &'static u16 {
debug_assert!(value < 100);
// Writes a significand and removes trailing zeros. value has up to 17 decimal // digits (16-17 for normals) for double (num_bits == 64) and up to 9 digits // (8-9 for normals) for float. The significant digits start from buffer[1]. // buffer[0] may contain '0' after this function if the leading digit is zero. #[cfg_attr(feature = "no-panic", no_panic)] #[inline] unsafefn write_significand<Float>(mut buffer: *mut u8, value: u64, extra_digit: bool) -> *mutu8 where
Float: FloatTraits,
{ if Float::NUM_BITS == 32 {
buffer = unsafe { write_if(buffer, (value / 100_000_000) as u32, extra_digit) }; let bcd = to_bcd8(value % 100_000_000); unsafe {
write8(buffer, bcd + ZEROS); return buffer.add(count_trailing_nonzeros(bcd));
}
}
#[cfg(not(any(
all(target_arch = "aarch64", target_feature = "neon", not(miri)),
all(target_arch = "x86_64", target_feature = "sse2", not(miri)),
)))]
{ // Digits/pairs of digits are denoted by letters: value = abbccddeeffgghhii. let abbccddee = (value / 100_000_000) as u32; let ffgghhii = (value % 100_000_000) as u32;
buffer = unsafe { write_if(buffer, abbccddee / 100_000_000, extra_digit) }; let bcd = to_bcd8(u64::from(abbccddee % 100_000_000)); unsafe {
write8(buffer, bcd + ZEROS);
} if ffgghhii == 0 { returnunsafe { buffer.add(count_trailing_nonzeros(bcd)) };
} let bcd = to_bcd8(u64::from(ffgghhii)); unsafe {
write8(buffer.add(8), bcd + ZEROS);
buffer.add(8).add(count_trailing_nonzeros(bcd))
}
}
#[cfg(all(target_arch = "aarch64", target_feature = "neon", not(miri)))]
{ // An optimized version for NEON by Dougall Johnson.
// Compiler barrier, or clang doesn't load from memory and generates 15 // more instructions. let c = unsafe {
asm!("/*{0}*/", inout(reg) c);
&*c
};
letmut hundred_million = c.hundred_million;
// Compiler barrier, or clang narrows the load to 32-bit and unpairs it. unsafe {
asm!("/*{0}*/", inout(reg) hundred_million);
}
// Equivalent to abbccddee = value / 100000000, ffgghhii = value % 100000000. let abbccddee = (umul128(value, c.mul_const) >> 90) as u64; let ffgghhii = value - abbccddee * hundred_million;
// We could probably make this bit faster, but we're preferring to // reuse the constants for now. let a = (umul128(abbccddee, c.mul_const) >> 90) as u64; let bbccddee = abbccddee - a * hundred_million;
buffer = unsafe { write_if(buffer, a as u32, extra_digit) };
unsafe { let ffgghhii_bbccddee_64: uint64x1_t =
mem::transmute::<u64, uint64x1_t>((ffgghhii << 32) | bbccddee); let bbccddee_ffgghhii: int32x2_t = vreinterpret_s32_u64(ffgghhii_bbccddee_64);
let abbccddee = (value / 100_000_000) as u32; let ffgghhii = (value % 100_000_000) as u32; let a = abbccddee / 100_000_000; let bbccddee = abbccddee % 100_000_000;
buffer = unsafe { write_if(buffer, a, extra_digit) };
letmut c = ptr::addr_of!(CONSTS); // Load constants from memory. unsafe {
asm!("/*{0}*/", inout(reg) c);
}
let div10k = unsafe { _mm_load_si128(ptr::addr_of!((*c).div10k).cast::<__m128i>()) }; let neg10k = unsafe { _mm_load_si128(ptr::addr_of!((*c).neg10k).cast::<__m128i>()) }; let div100 = unsafe { _mm_load_si128(ptr::addr_of!((*c).div100).cast::<__m128i>()) }; let div10 = unsafe { _mm_load_si128(ptr::addr_of!((*c).div10).cast::<__m128i>()) }; #[cfg(target_feature = "sse4.1")] let neg100 = unsafe { _mm_load_si128(ptr::addr_of!((*c).neg100).cast::<__m128i>()) }; #[cfg(target_feature = "sse4.1")] let neg10 = unsafe { _mm_load_si128(ptr::addr_of!((*c).neg10).cast::<__m128i>()) }; #[cfg(target_feature = "sse4.1")] let bswap = unsafe { _mm_load_si128(ptr::addr_of!((*c).bswap).cast::<__m128i>()) }; #[cfg(not(target_feature = "sse4.1"))] let hundred = unsafe { _mm_load_si128(ptr::addr_of!((*c).hundred).cast::<__m128i>()) }; #[cfg(not(target_feature = "sse4.1"))] let moddiv10 = unsafe { _mm_load_si128(ptr::addr_of!((*c).moddiv10).cast::<__m128i>()) }; let zeros = unsafe { _mm_load_si128(ptr::addr_of!((*c).zeros).cast::<__m128i>()) };
// The BCD sequences are based on ones provided by Xiang JunBo. unsafe { let x: __m128i = _mm_set_epi64x(i64::from(bbccddee), i64::from(ffgghhii)); let y: __m128i = _mm_add_epi64(
x,
_mm_mul_epu32(neg10k, _mm_srli_epi64(_mm_mul_epu32(x, div10k), DIV10K_EXP)),
);
#[cfg(target_feature = "sse4.1")] let bcd: __m128i = { // _mm_mullo_epi32 is SSE 4.1 let z: __m128i = _mm_add_epi64(
y,
_mm_mullo_epi32(neg100, _mm_srli_epi32(_mm_mulhi_epu16(y, div100), 3)),
); let big_endian_bcd: __m128i =
_mm_add_epi64(z, _mm_mullo_epi16(neg10, _mm_mulhi_epu16(z, div10))); // SSSE3
_mm_shuffle_epi8(big_endian_bcd, bswap)
};
#[cfg(not(target_feature = "sse4.1"))] let bcd: __m128i = { let y_div_100: __m128i = _mm_srli_epi16(_mm_mulhi_epu16(y, div100), 3); let y_mod_100: __m128i = _mm_sub_epi16(y, _mm_mullo_epi16(y_div_100, hundred)); let z: __m128i = _mm_or_si128(_mm_slli_epi32(y_mod_100, 16), y_div_100); let bcd_shuffled: __m128i = _mm_sub_epi16(
_mm_slli_epi16(z, 8),
_mm_mullo_epi16(moddiv10, _mm_mulhi_epu16(z, div10)),
);
_mm_shuffle_epi32(bcd_shuffled, _MM_SHUFFLE(0, 1, 2, 3))
};
let digits = _mm_or_si128(bcd, zeros);
// Count leading zeros. let mask128: __m128i = _mm_cmpgt_epi8(bcd, _mm_setzero_si128()); let mask = _mm_movemask_epi8(mask128) as u32; let len = 32 - mask.leading_zeros() as usize;
#[cfg_attr(feature = "no-panic", no_panic)] #[inline] fn to_decimal_schubfach<UInt>(bin_sig: UInt, bin_exp: i64, regular: bool) -> ToDecimalResult where
UInt: traits::UInt,
{ let num_bits = mem::size_of::<UInt>() as i32 * 8; let dec_exp = compute_dec_exp(bin_exp as i32, regular); let exp_shift = unsafe { compute_exp_shift::<UInt, false>(bin_exp as i32, dec_exp) }; letmut pow10 = unsafe { POW10_SIGNIFICANDS.get_unchecked(-dec_exp) };
// Fallback to Schubfach to guarantee correctness in boundary cases. This // requires switching to strict overestimates of powers of 10. if num_bits == 64 {
pow10.lo += 1;
} else {
pow10.hi += 1;
}
// Shift the significand so that boundaries are integer. const BOUND_SHIFT: u32 = 2; let bin_sig_shifted = bin_sig << BOUND_SHIFT;
// Compute the estimates of lower and upper bounds of the rounding interval // by multiplying them by the power of 10 and applying modified rounding. let lsb = bin_sig & UInt::from(1); let lower = (bin_sig_shifted - (UInt::from(regular) + UInt::from(1))) << exp_shift; let lower = umulhi_inexact_to_odd(pow10.hi, pow10.lo, lower) + lsb; let upper = (bin_sig_shifted + UInt::from(2)) << exp_shift; let upper = umulhi_inexact_to_odd(pow10.hi, pow10.lo, upper) - lsb;
// The idea of using a single shorter candidate is by Cassio Neri. // It is less or equal to the upper bound by construction. let shorter = (upper >> BOUND_SHIFT) / UInt::from(10) * UInt::from(10); if (shorter << BOUND_SHIFT) >= lower { return ToDecimalResult {
sig: shorter.into() as i64,
exp: dec_exp,
};
}
let scaled_sig = umulhi_inexact_to_odd(pow10.hi, pow10.lo, bin_sig_shifted << exp_shift); let longer_below = scaled_sig >> BOUND_SHIFT; let longer_above = longer_below + UInt::from(1);
// Pick the closest of longer_below and longer_above and check if it's in // the rounding interval. let cmp = scaled_sig
.wrapping_sub((longer_below + longer_above) << 1)
.to_signed(); let below_closer = cmp < UInt::from(0).to_signed()
|| (cmp == UInt::from(0).to_signed() && (longer_below & UInt::from(1)) == UInt::from(0)); let below_in = (longer_below << BOUND_SHIFT) >= lower; let dec_sig = if below_closer & below_in {
longer_below
} else {
longer_above
};
ToDecimalResult {
sig: dec_sig.into() as i64,
exp: dec_exp,
}
}
// Here be s. // Converts a binary FP number bin_sig * 2**bin_exp to the shortest decimal // representation, where bin_exp = raw_exp - exp_offset. #[cfg_attr(feature = "no-panic", no_panic)] #[inline] fn to_decimal_fast<Float, UInt>(bin_sig: UInt, raw_exp: i64, regular: bool) -> ToDecimalResult where
Float: FloatTraits,
UInt: traits::UInt,
{ let bin_exp = raw_exp - i64::from(Float::EXP_OFFSET); let num_bits = mem::size_of::<UInt>() as i32 * 8; // An optimization from yy by Yaoyuan Guo: while regular { let dec_exp = if USE_UMUL128_HI64 {
umul128_hi64(bin_exp as u64, 0x4d10500000000000) as i32
} else {
compute_dec_exp(bin_exp as i32, true)
}; let exp_shift = unsafe { compute_exp_shift::<UInt, true>(bin_exp as i32, dec_exp) }; let pow10 = unsafe { POW10_SIGNIFICANDS.get_unchecked(-dec_exp) };
let integral; // integral part of bin_sig * pow10 let fractional; // fractional part of bin_sig * pow10 if num_bits == 64 { let p = umul192_hi128(pow10.hi, pow10.lo, (bin_sig << exp_shift).into());
integral = UInt::truncate(p.hi);
fractional = p.lo;
} else { let p = umul128(pow10.hi, (bin_sig << exp_shift).into());
integral = UInt::truncate((p >> 64) as u64);
fractional = p as u64;
} const HALF_ULP: u64 = 1 << 63;
// Exact half-ulp tie when rounding to nearest integer. let cmp = fractional.wrapping_sub(HALF_ULP) as i64; if cmp == 0 { break;
}
// An optimization of integral % 10 by Dougall Johnson. Relies on range // calculation: (max_bin_sig << max_exp_shift) * max_u128. // (1 << 63) / 5 == (1 << 64) / 10 without an intermediate int128. const DIV10_SIG64: u64 = (1 << 63) / 5 + 1; let div10 = umul128_hi64(integral.into(), DIV10_SIG64); #[allow(unused_mut)] letmut digit = integral.into() - div10 * 10; // or it narrows to 32-bit and doesn't use madd/msub #[cfg(all(any(target_arch = "aarch64", target_arch = "x86_64"), not(miri)))] unsafe {
asm!("/*{0}*/", inout(reg) digit);
}
// Switch to a fixed-point representation with the least significant // integral digit in the upper bits and fractional digits in the lower // bits. let num_integral_bits = if num_bits == 64 { 4 } else { 32 }; let num_fractional_bits = 64 - num_integral_bits; let ten = 10u64 << num_fractional_bits; // Fixed-point remainder of the scaled significand modulo 10. let scaled_sig_mod10 = (digit << num_fractional_bits) | (fractional >> num_integral_bits);
// scaled_half_ulp = 0.5 * pow10 in the fixed-point format. // dec_exp is chosen so that 10**dec_exp <= 2**bin_exp < 10**(dec_exp + 1). // Since 1ulp == 2**bin_exp it will be in the range [1, 10) after scaling // by 10**dec_exp. Add 1 to combine the shift with division by two. let scaled_half_ulp = pow10.hi >> (num_integral_bits - exp_shift + 1); let upper = scaled_sig_mod10 + scaled_half_ulp;
// value = 5.0507837461e-27 // next = 5.0507837461000010e-27 // // c = integral.fractional' = 50507837461000003.153987... (value) // 50507837461000010.328635... (next) // scaled_half_ulp = 3.587324... // // fractional' = fractional / 2**64, fractional = 2840565642863009226 // // 50507837461000000 c upper 50507837461000010 // s l| L | S // ───┬────┬────┼────┬────┬────┼*-──┼────┬────┬───*┬────┬────┬────┼-*--┬─── // 8 9 0 1 2 3 4 5 6 7 8 9 0 | 1 // └─────────────────┼─────────────────┘ next // 1ulp // // s - shorter underestimate, S - shorter overestimate // l - longer underestimate, L - longer overestimate
// Check for boundary case when rounding down to nearest 10 and // near-boundary case when rounding up to nearest 10. // Case where upper == ten is insufficient: 1.342178e+08f. if ten.wrapping_sub(upper) <= 1// upper == ten || upper == ten - 1
|| scaled_sig_mod10 == scaled_half_ulp
{ break;
}
let shorter = (integral.into() - digit) as i64; let longer = (integral.into() + u64::from(cmp >= 0)) as i64; let dec_sig = select_if_less(scaled_sig_mod10, scaled_half_ulp, shorter, longer); return ToDecimalResult {
sig: select_if_less(ten, upper, shorter + 10, dec_sig),
exp: dec_exp,
};
}
to_decimal_schubfach(bin_sig, bin_exp, regular)
}
/// Writes the shortest correctly rounded decimal representation of `value` to /// `buffer`. `buffer` should point to a buffer of size `buffer_size` or larger. #[cfg_attr(feature = "no-panic", no_panic)] unsafefn write<Float>(value: Float, mut buffer: *mut u8) -> *mut u8 where
Float: FloatTraits,
{ let bits = value.to_bits(); // It is beneficial to extract exponent and significand early. let bin_exp = Float::get_exp(bits); // binary exponent let bin_sig = Float::get_sig(bits); // binary significand
/// Safe API for formatting floating point numbers to text. /// /// ## Example /// /// ``` /// let mut buffer = zmij::Buffer::new(); /// let printed = buffer.format_finite(1.234); /// assert_eq!(printed, "1.234"); /// ``` pubstruct Buffer {
bytes: [MaybeUninit<u8>; BUFFER_SIZE],
}
impl Buffer { /// This is a cheap operation; you don't need to worry about reusing buffers /// for efficiency. #[inline] #[cfg_attr(feature = "no-panic", no_panic)] pubfn new() -> Self { let bytes = [MaybeUninit::<u8>::uninit(); BUFFER_SIZE];
Buffer { bytes }
}
/// Print a floating point number into this buffer and return a reference to /// its string representation within the buffer. /// /// # Special cases /// /// This function formats NaN as the string "NaN", positive infinity as /// "inf", and negative infinity as "-inf" to match std::fmt. /// /// If your input is known to be finite, you may get better performance by /// calling the `format_finite` method instead of `format` to avoid the /// checks for special cases. #[cfg_attr(feature = "no-panic", no_panic)] pubfn format<F: Float>(&mutself, f: F) -> &str { if f.is_nonfinite() {
f.format_nonfinite()
} else { self.format_finite(f)
}
}
/// Print a floating point number into this buffer and return a reference to /// its string representation within the buffer. /// /// # Special cases /// /// This function **does not** check for NaN or infinity. If the input /// number is not a finite float, the printed representation will be some /// correctly formatted but unspecified numerical value. /// /// Please check [`is_finite`] yourself before calling this function, or /// check [`is_nan`] and [`is_infinite`] and handle those cases yourself. /// /// [`is_finite`]: f64::is_finite /// [`is_nan`]: f64::is_nan /// [`is_infinite`]: f64::is_infinite #[cfg_attr(feature = "no-panic", no_panic)] pubfn format_finite<F: Float>(&mutself, f: F) -> &str { unsafe { let end = f.write_to_zmij_buffer(self.bytes.as_mut_ptr().cast::<u8>()); let len = end.offset_from(self.bytes.as_ptr().cast::<u8>()) as usize; let slice = slice::from_raw_parts(self.bytes.as_ptr().cast::<u8>(), len);
str::from_utf8_unchecked(slice)
}
}
}
/// A floating point number, f32 or f64, that can be written into a /// [`zmij::Buffer`][Buffer]. /// /// This trait is sealed and cannot be implemented for types outside of the /// `zmij` crate. #[allow(unknown_lints)] // rustc older than 1.74 #[allow(private_bounds)] pubtrait Float: private::Sealed {} impl Float for f32 {} impl Float for f64 {}
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.