blob: 4ad50a1bce5a17acc750c92ee926fff1966d9cd6 [file] [log] [blame]
// Translated from C to Rust. The original C code can be found at
// https://github.com/ulfjack/ryu and carries the following license:
//
// Copyright 2018 Ulf Adams
//
// The contents of this file may be used under the terms of the Apache License,
// Version 2.0.
//
// (See accompanying file LICENSE-Apache or copy at
// http://www.apache.org/licenses/LICENSE-2.0)
//
// Alternatively, the contents of this file may be used under the terms of
// the Boost Software License, Version 1.0.
// (See accompanying file LICENSE-Boost or copy at
// https://www.boost.org/LICENSE_1_0.txt)
//
// Unless required by applicable law or agreed to in writing, this software
// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied.
use common::*;
pub const FLOAT_MANTISSA_BITS: u32 = 23;
pub const FLOAT_EXPONENT_BITS: u32 = 8;
const FLOAT_BIAS: i32 = 127;
const FLOAT_POW5_INV_BITCOUNT: i32 = 59;
const FLOAT_POW5_BITCOUNT: i32 = 61;
// This table is generated by PrintFloatLookupTable.
static FLOAT_POW5_INV_SPLIT: [u64; 32] = [
576460752303423489,
461168601842738791,
368934881474191033,
295147905179352826,
472236648286964522,
377789318629571618,
302231454903657294,
483570327845851670,
386856262276681336,
309485009821345069,
495176015714152110,
396140812571321688,
316912650057057351,
507060240091291761,
405648192073033409,
324518553658426727,
519229685853482763,
415383748682786211,
332306998946228969,
531691198313966350,
425352958651173080,
340282366920938464,
544451787073501542,
435561429658801234,
348449143727040987,
557518629963265579,
446014903970612463,
356811923176489971,
570899077082383953,
456719261665907162,
365375409332725730,
1 << 63,
];
static FLOAT_POW5_SPLIT: [u64; 47] = [
1152921504606846976,
1441151880758558720,
1801439850948198400,
2251799813685248000,
1407374883553280000,
1759218604441600000,
2199023255552000000,
1374389534720000000,
1717986918400000000,
2147483648000000000,
1342177280000000000,
1677721600000000000,
2097152000000000000,
1310720000000000000,
1638400000000000000,
2048000000000000000,
1280000000000000000,
1600000000000000000,
2000000000000000000,
1250000000000000000,
1562500000000000000,
1953125000000000000,
1220703125000000000,
1525878906250000000,
1907348632812500000,
1192092895507812500,
1490116119384765625,
1862645149230957031,
1164153218269348144,
1455191522836685180,
1818989403545856475,
2273736754432320594,
1421085471520200371,
1776356839400250464,
2220446049250313080,
1387778780781445675,
1734723475976807094,
2168404344971008868,
1355252715606880542,
1694065894508600678,
2117582368135750847,
1323488980084844279,
1654361225106055349,
2067951531382569187,
1292469707114105741,
1615587133892632177,
2019483917365790221,
];
#[cfg_attr(feature = "no-panic", inline)]
fn pow5_factor(mut value: u32) -> u32 {
let mut count = 0u32;
loop {
debug_assert!(value != 0);
let q = value / 5;
let r = value % 5;
if r != 0 {
break;
}
value = q;
count += 1;
}
count
}
// Returns true if value is divisible by 5^p.
#[cfg_attr(feature = "no-panic", inline)]
fn multiple_of_power_of_5(value: u32, p: u32) -> bool {
pow5_factor(value) >= p
}
// Returns true if value is divisible by 2^p.
#[cfg_attr(feature = "no-panic", inline)]
fn multiple_of_power_of_2(value: u32, p: u32) -> bool {
// return __builtin_ctz(value) >= p;
(value & ((1u32 << p) - 1)) == 0
}
// It seems to be slightly faster to avoid uint128_t here, although the
// generated code for uint128_t looks slightly nicer.
#[cfg_attr(feature = "no-panic", inline)]
fn mul_shift(m: u32, factor: u64, shift: i32) -> u32 {
debug_assert!(shift > 32);
// The casts here help MSVC to avoid calls to the __allmul library
// function.
let factor_lo = factor as u32;
let factor_hi = (factor >> 32) as u32;
let bits0 = m as u64 * factor_lo as u64;
let bits1 = m as u64 * factor_hi as u64;
let sum = (bits0 >> 32) + bits1;
let shifted_sum = sum >> (shift - 32);
debug_assert!(shifted_sum <= u32::max_value() as u64);
shifted_sum as u32
}
#[cfg_attr(feature = "no-panic", inline)]
fn mul_pow5_inv_div_pow2(m: u32, q: u32, j: i32) -> u32 {
debug_assert!(q < FLOAT_POW5_INV_SPLIT.len() as u32);
unsafe { mul_shift(m, *FLOAT_POW5_INV_SPLIT.get_unchecked(q as usize), j) }
}
#[cfg_attr(feature = "no-panic", inline)]
fn mul_pow5_div_pow2(m: u32, i: u32, j: i32) -> u32 {
debug_assert!(i < FLOAT_POW5_SPLIT.len() as u32);
unsafe { mul_shift(m, *FLOAT_POW5_SPLIT.get_unchecked(i as usize), j) }
}
// A floating decimal representing m * 10^e.
pub struct FloatingDecimal32 {
pub mantissa: u32,
// Decimal exponent's range is -45 to 38
// inclusive, and can fit in i16 if needed.
pub exponent: i32,
}
#[cfg_attr(feature = "no-panic", inline)]
pub fn f2d(ieee_mantissa: u32, ieee_exponent: u32) -> FloatingDecimal32 {
let (e2, m2) = if ieee_exponent == 0 {
(
// We subtract 2 so that the bounds computation has 2 additional bits.
1 - FLOAT_BIAS - FLOAT_MANTISSA_BITS as i32 - 2,
ieee_mantissa,
)
} else {
(
ieee_exponent as i32 - FLOAT_BIAS - FLOAT_MANTISSA_BITS as i32 - 2,
(1u32 << FLOAT_MANTISSA_BITS) | ieee_mantissa,
)
};
let even = (m2 & 1) == 0;
let accept_bounds = even;
// Step 2: Determine the interval of valid decimal representations.
let mv = 4 * m2;
let mp = 4 * m2 + 2;
// Implicit bool -> int conversion. True is 1, false is 0.
let mm_shift = (ieee_mantissa != 0 || ieee_exponent <= 1) as u32;
let mm = 4 * m2 - 1 - mm_shift;
// Step 3: Convert to a decimal power base using 64-bit arithmetic.
let mut vr: u32;
let mut vp: u32;
let mut vm: u32;
let e10: i32;
let mut vm_is_trailing_zeros = false;
let mut vr_is_trailing_zeros = false;
let mut last_removed_digit = 0u8;
if e2 >= 0 {
let q = log10_pow2(e2);
e10 = q as i32;
let k = FLOAT_POW5_INV_BITCOUNT + pow5bits(q as i32) - 1;
let i = -e2 + q as i32 + k;
vr = mul_pow5_inv_div_pow2(mv, q, i);
vp = mul_pow5_inv_div_pow2(mp, q, i);
vm = mul_pow5_inv_div_pow2(mm, q, i);
if q != 0 && (vp - 1) / 10 <= vm / 10 {
// We need to know one removed digit even if we are not going to loop below. We could use
// q = X - 1 above, except that would require 33 bits for the result, and we've found that
// 32-bit arithmetic is faster even on 64-bit machines.
let l = FLOAT_POW5_INV_BITCOUNT + pow5bits(q as i32 - 1) - 1;
last_removed_digit =
(mul_pow5_inv_div_pow2(mv, q - 1, -e2 + q as i32 - 1 + l) % 10) as u8;
}
if q <= 9 {
// The largest power of 5 that fits in 24 bits is 5^10, but q <= 9 seems to be safe as well.
// Only one of mp, mv, and mm can be a multiple of 5, if any.
if mv % 5 == 0 {
vr_is_trailing_zeros = multiple_of_power_of_5(mv, q);
} else if accept_bounds {
vm_is_trailing_zeros = multiple_of_power_of_5(mm, q);
} else {
vp -= multiple_of_power_of_5(mp, q) as u32;
}
}
} else {
let q = log10_pow5(-e2);
e10 = q as i32 + e2;
let i = -e2 - q as i32;
let k = pow5bits(i) - FLOAT_POW5_BITCOUNT;
let mut j = q as i32 - k;
vr = mul_pow5_div_pow2(mv, i as u32, j);
vp = mul_pow5_div_pow2(mp, i as u32, j);
vm = mul_pow5_div_pow2(mm, i as u32, j);
if q != 0 && (vp - 1) / 10 <= vm / 10 {
j = q as i32 - 1 - (pow5bits(i + 1) - FLOAT_POW5_BITCOUNT);
last_removed_digit = (mul_pow5_div_pow2(mv, (i + 1) as u32, j) % 10) as u8;
}
if q <= 1 {
// {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits.
// mv = 4 * m2, so it always has at least two trailing 0 bits.
vr_is_trailing_zeros = true;
if accept_bounds {
// mm = mv - 1 - mm_shift, so it has 1 trailing 0 bit iff mm_shift == 1.
vm_is_trailing_zeros = mm_shift == 1;
} else {
// mp = mv + 2, so it always has at least one trailing 0 bit.
vp -= 1;
}
} else if q < 31 {
// TODO(ulfjack): Use a tighter bound here.
vr_is_trailing_zeros = multiple_of_power_of_2(mv, q - 1);
}
}
// Step 4: Find the shortest decimal representation in the interval of valid representations.
let mut removed = 0i32;
let output = if vm_is_trailing_zeros || vr_is_trailing_zeros {
// General case, which happens rarely (~4.0%).
while vp / 10 > vm / 10 {
vm_is_trailing_zeros &= vm - (vm / 10) * 10 == 0;
vr_is_trailing_zeros &= last_removed_digit == 0;
last_removed_digit = (vr % 10) as u8;
vr /= 10;
vp /= 10;
vm /= 10;
removed += 1;
}
if vm_is_trailing_zeros {
while vm % 10 == 0 {
vr_is_trailing_zeros &= last_removed_digit == 0;
last_removed_digit = (vr % 10) as u8;
vr /= 10;
vp /= 10;
vm /= 10;
removed += 1;
}
}
if vr_is_trailing_zeros && last_removed_digit == 5 && vr % 2 == 0 {
// Round even if the exact number is .....50..0.
last_removed_digit = 4;
}
// We need to take vr + 1 if vr is outside bounds or we need to round up.
vr + ((vr == vm && (!accept_bounds || !vm_is_trailing_zeros)) || last_removed_digit >= 5)
as u32
} else {
// Specialized for the common case (~96.0%). Percentages below are relative to this.
// Loop iterations below (approximately):
// 0: 13.6%, 1: 70.7%, 2: 14.1%, 3: 1.39%, 4: 0.14%, 5+: 0.01%
while vp / 10 > vm / 10 {
last_removed_digit = (vr % 10) as u8;
vr /= 10;
vp /= 10;
vm /= 10;
removed += 1;
}
// We need to take vr + 1 if vr is outside bounds or we need to round up.
vr + (vr == vm || last_removed_digit >= 5) as u32
};
let exp = e10 + removed;
FloatingDecimal32 {
exponent: exp,
mantissa: output,
}
}