ryu/
f2s.rs

Help
1// Translated from C to Rust. The original C code can be found at
2// https://github.com/ulfjack/ryu and carries the following license:
3//
4// Copyright 2018 Ulf Adams
5//
6// The contents of this file may be used under the terms of the Apache License,
7// Version 2.0.
8//
9//    (See accompanying file LICENSE-Apache or copy at
10//     http://www.apache.org/licenses/LICENSE-2.0)
11//
12// Alternatively, the contents of this file may be used under the terms of
13// the Boost Software License, Version 1.0.
14//    (See accompanying file LICENSE-Boost or copy at
15//     https://www.boost.org/LICENSE_1_0.txt)
16//
17// Unless required by applicable law or agreed to in writing, this software
18// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
19// KIND, either express or implied.
20
21use common::*;
22
23pub const FLOAT_MANTISSA_BITS: u32 = 23;
24pub const FLOAT_EXPONENT_BITS: u32 = 8;
25
26const FLOAT_BIAS: i32 = 127;
27const FLOAT_POW5_INV_BITCOUNT: i32 = 59;
28const FLOAT_POW5_BITCOUNT: i32 = 61;
29
30// This table is generated by PrintFloatLookupTable.
31static FLOAT_POW5_INV_SPLIT: [u64; 32] = [
32    576460752303423489,
33    461168601842738791,
34    368934881474191033,
35    295147905179352826,
36    472236648286964522,
37    377789318629571618,
38    302231454903657294,
39    483570327845851670,
40    386856262276681336,
41    309485009821345069,
42    495176015714152110,
43    396140812571321688,
44    316912650057057351,
45    507060240091291761,
46    405648192073033409,
47    324518553658426727,
48    519229685853482763,
49    415383748682786211,
50    332306998946228969,
51    531691198313966350,
52    425352958651173080,
53    340282366920938464,
54    544451787073501542,
55    435561429658801234,
56    348449143727040987,
57    557518629963265579,
58    446014903970612463,
59    356811923176489971,
60    570899077082383953,
61    456719261665907162,
62    365375409332725730,
63    1 << 63,
64];
65
66static FLOAT_POW5_SPLIT: [u64; 47] = [
67    1152921504606846976,
68    1441151880758558720,
69    1801439850948198400,
70    2251799813685248000,
71    1407374883553280000,
72    1759218604441600000,
73    2199023255552000000,
74    1374389534720000000,
75    1717986918400000000,
76    2147483648000000000,
77    1342177280000000000,
78    1677721600000000000,
79    2097152000000000000,
80    1310720000000000000,
81    1638400000000000000,
82    2048000000000000000,
83    1280000000000000000,
84    1600000000000000000,
85    2000000000000000000,
86    1250000000000000000,
87    1562500000000000000,
88    1953125000000000000,
89    1220703125000000000,
90    1525878906250000000,
91    1907348632812500000,
92    1192092895507812500,
93    1490116119384765625,
94    1862645149230957031,
95    1164153218269348144,
96    1455191522836685180,
97    1818989403545856475,
98    2273736754432320594,
99    1421085471520200371,
100    1776356839400250464,
101    2220446049250313080,
102    1387778780781445675,
103    1734723475976807094,
104    2168404344971008868,
105    1355252715606880542,
106    1694065894508600678,
107    2117582368135750847,
108    1323488980084844279,
109    1654361225106055349,
110    2067951531382569187,
111    1292469707114105741,
112    1615587133892632177,
113    2019483917365790221,
114];
115
116#[cfg_attr(feature = "no-panic", inline)]
117fn pow5_factor(mut value: u32) -> u32 {
118    let mut count = 0u32;
119    loop {
120        debug_assert!(value != 0);
121        let q = value / 5;
122        let r = value % 5;
123        if r != 0 {
124            break;
125        }
126        value = q;
127        count += 1;
128    }
129    count
130}
131
132// Returns true if value is divisible by 5^p.
133#[cfg_attr(feature = "no-panic", inline)]
134fn multiple_of_power_of_5(value: u32, p: u32) -> bool {
135    pow5_factor(value) >= p
136}
137
138// Returns true if value is divisible by 2^p.
139#[cfg_attr(feature = "no-panic", inline)]
140fn multiple_of_power_of_2(value: u32, p: u32) -> bool {
141    // return __builtin_ctz(value) >= p;
142    (value & ((1u32 << p) - 1)) == 0
143}
144
145// It seems to be slightly faster to avoid uint128_t here, although the
146// generated code for uint128_t looks slightly nicer.
147#[cfg_attr(feature = "no-panic", inline)]
148fn mul_shift(m: u32, factor: u64, shift: i32) -> u32 {
149    debug_assert!(shift > 32);
150
151    // The casts here help MSVC to avoid calls to the __allmul library
152    // function.
153    let factor_lo = factor as u32;
154    let factor_hi = (factor >> 32) as u32;
155    let bits0 = m as u64 * factor_lo as u64;
156    let bits1 = m as u64 * factor_hi as u64;
157
158    let sum = (bits0 >> 32) + bits1;
159    let shifted_sum = sum >> (shift - 32);
160    debug_assert!(shifted_sum <= u32::max_value() as u64);
161    shifted_sum as u32
162}
163
164#[cfg_attr(feature = "no-panic", inline)]
165fn mul_pow5_inv_div_pow2(m: u32, q: u32, j: i32) -> u32 {
166    debug_assert!(q < FLOAT_POW5_INV_SPLIT.len() as u32);
167    unsafe { mul_shift(m, *FLOAT_POW5_INV_SPLIT.get_unchecked(q as usize), j) }
168}
169
170#[cfg_attr(feature = "no-panic", inline)]
171fn mul_pow5_div_pow2(m: u32, i: u32, j: i32) -> u32 {
172    debug_assert!(i < FLOAT_POW5_SPLIT.len() as u32);
173    unsafe { mul_shift(m, *FLOAT_POW5_SPLIT.get_unchecked(i as usize), j) }
174}
175
176// A floating decimal representing m * 10^e.
177pub struct FloatingDecimal32 {
178    pub mantissa: u32,
179    // Decimal exponent's range is -45 to 38
180    // inclusive, and can fit in i16 if needed.
181    pub exponent: i32,
182}
183
184#[cfg_attr(feature = "no-panic", inline)]
185pub fn f2d(ieee_mantissa: u32, ieee_exponent: u32) -> FloatingDecimal32 {
186    let (e2, m2) = if ieee_exponent == 0 {
187        (
188            // We subtract 2 so that the bounds computation has 2 additional bits.
189            1 - FLOAT_BIAS - FLOAT_MANTISSA_BITS as i32 - 2,
190            ieee_mantissa,
191        )
192    } else {
193        (
194            ieee_exponent as i32 - FLOAT_BIAS - FLOAT_MANTISSA_BITS as i32 - 2,
195            (1u32 << FLOAT_MANTISSA_BITS) | ieee_mantissa,
196        )
197    };
198    let even = (m2 & 1) == 0;
199    let accept_bounds = even;
200
201    // Step 2: Determine the interval of valid decimal representations.
202    let mv = 4 * m2;
203    let mp = 4 * m2 + 2;
204    // Implicit bool -> int conversion. True is 1, false is 0.
205    let mm_shift = (ieee_mantissa != 0 || ieee_exponent <= 1) as u32;
206    let mm = 4 * m2 - 1 - mm_shift;
207
208    // Step 3: Convert to a decimal power base using 64-bit arithmetic.
209    let mut vr: u32;
210    let mut vp: u32;
211    let mut vm: u32;
212    let e10: i32;
213    let mut vm_is_trailing_zeros = false;
214    let mut vr_is_trailing_zeros = false;
215    let mut last_removed_digit = 0u8;
216    if e2 >= 0 {
217        let q = log10_pow2(e2);
218        e10 = q as i32;
219        let k = FLOAT_POW5_INV_BITCOUNT + pow5bits(q as i32) - 1;
220        let i = -e2 + q as i32 + k;
221        vr = mul_pow5_inv_div_pow2(mv, q, i);
222        vp = mul_pow5_inv_div_pow2(mp, q, i);
223        vm = mul_pow5_inv_div_pow2(mm, q, i);
224        if q != 0 && (vp - 1) / 10 <= vm / 10 {
225            // We need to know one removed digit even if we are not going to loop below. We could use
226            // q = X - 1 above, except that would require 33 bits for the result, and we've found that
227            // 32-bit arithmetic is faster even on 64-bit machines.
228            let l = FLOAT_POW5_INV_BITCOUNT + pow5bits(q as i32 - 1) - 1;
229            last_removed_digit =
230                (mul_pow5_inv_div_pow2(mv, q - 1, -e2 + q as i32 - 1 + l) % 10) as u8;
231        }
232        if q <= 9 {
233            // The largest power of 5 that fits in 24 bits is 5^10, but q <= 9 seems to be safe as well.
234            // Only one of mp, mv, and mm can be a multiple of 5, if any.
235            if mv % 5 == 0 {
236                vr_is_trailing_zeros = multiple_of_power_of_5(mv, q);
237            } else if accept_bounds {
238                vm_is_trailing_zeros = multiple_of_power_of_5(mm, q);
239            } else {
240                vp -= multiple_of_power_of_5(mp, q) as u32;
241            }
242        }
243    } else {
244        let q = log10_pow5(-e2);
245        e10 = q as i32 + e2;
246        let i = -e2 - q as i32;
247        let k = pow5bits(i) - FLOAT_POW5_BITCOUNT;
248        let mut j = q as i32 - k;
249        vr = mul_pow5_div_pow2(mv, i as u32, j);
250        vp = mul_pow5_div_pow2(mp, i as u32, j);
251        vm = mul_pow5_div_pow2(mm, i as u32, j);
252        if q != 0 && (vp - 1) / 10 <= vm / 10 {
253            j = q as i32 - 1 - (pow5bits(i + 1) - FLOAT_POW5_BITCOUNT);
254            last_removed_digit = (mul_pow5_div_pow2(mv, (i + 1) as u32, j) % 10) as u8;
255        }
256        if q <= 1 {
257            // {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits.
258            // mv = 4 * m2, so it always has at least two trailing 0 bits.
259            vr_is_trailing_zeros = true;
260            if accept_bounds {
261                // mm = mv - 1 - mm_shift, so it has 1 trailing 0 bit iff mm_shift == 1.
262                vm_is_trailing_zeros = mm_shift == 1;
263            } else {
264                // mp = mv + 2, so it always has at least one trailing 0 bit.
265                vp -= 1;
266            }
267        } else if q < 31 {
268            // TODO(ulfjack): Use a tighter bound here.
269            vr_is_trailing_zeros = multiple_of_power_of_2(mv, q - 1);
270        }
271    }
272
273    // Step 4: Find the shortest decimal representation in the interval of valid representations.
274    let mut removed = 0i32;
275    let output = if vm_is_trailing_zeros || vr_is_trailing_zeros {
276        // General case, which happens rarely (~4.0%).
277        while vp / 10 > vm / 10 {
278            vm_is_trailing_zeros &= vm - (vm / 10) * 10 == 0;
279            vr_is_trailing_zeros &= last_removed_digit == 0;
280            last_removed_digit = (vr % 10) as u8;
281            vr /= 10;
282            vp /= 10;
283            vm /= 10;
284            removed += 1;
285        }
286        if vm_is_trailing_zeros {
287            while vm % 10 == 0 {
288                vr_is_trailing_zeros &= last_removed_digit == 0;
289                last_removed_digit = (vr % 10) as u8;
290                vr /= 10;
291                vp /= 10;
292                vm /= 10;
293                removed += 1;
294            }
295        }
296        if vr_is_trailing_zeros && last_removed_digit == 5 && vr % 2 == 0 {
297            // Round even if the exact number is .....50..0.
298            last_removed_digit = 4;
299        }
300        // We need to take vr + 1 if vr is outside bounds or we need to round up.
301        vr + ((vr == vm && (!accept_bounds || !vm_is_trailing_zeros)) || last_removed_digit >= 5)
302            as u32
303    } else {
304        // Specialized for the common case (~96.0%). Percentages below are relative to this.
305        // Loop iterations below (approximately):
306        // 0: 13.6%, 1: 70.7%, 2: 14.1%, 3: 1.39%, 4: 0.14%, 5+: 0.01%
307        while vp / 10 > vm / 10 {
308            last_removed_digit = (vr % 10) as u8;
309            vr /= 10;
310            vp /= 10;
311            vm /= 10;
312            removed += 1;
313        }
314        // We need to take vr + 1 if vr is outside bounds or we need to round up.
315        vr + (vr == vm || last_removed_digit >= 5) as u32
316    };
317    let exp = e10 + removed;
318
319    FloatingDecimal32 {
320        exponent: exp,
321        mantissa: output,
322    }
323}
ryu/f2s.rs

ryu/
f2s.rs
