ryu/
f2s.rs

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 crate::common::{log10_pow2, log10_pow5, pow5bits};
22use crate::f2s_intrinsics::{
23    mul_pow5_div_pow2, mul_pow5_inv_div_pow2, multiple_of_power_of_2_32, multiple_of_power_of_5_32,
24};
25
26pub const FLOAT_MANTISSA_BITS: u32 = 23;
27pub const FLOAT_EXPONENT_BITS: u32 = 8;
28const FLOAT_BIAS: i32 = 127;
29pub use crate::f2s_intrinsics::{FLOAT_POW5_BITCOUNT, FLOAT_POW5_INV_BITCOUNT};
30
31// A floating decimal representing m * 10^e.
32pub struct FloatingDecimal32 {
33    pub mantissa: u32,
34    // Decimal exponent's range is -45 to 38
35    // inclusive, and can fit in i16 if needed.
36    pub exponent: i32,
37}
38
39#[cfg_attr(feature = "no-panic", inline)]
40pub fn f2d(ieee_mantissa: u32, ieee_exponent: u32) -> FloatingDecimal32 {
41    let (e2, m2) = if ieee_exponent == 0 {
42        (
43            // We subtract 2 so that the bounds computation has 2 additional bits.
44            1 - FLOAT_BIAS - FLOAT_MANTISSA_BITS as i32 - 2,
45            ieee_mantissa,
46        )
47    } else {
48        (
49            ieee_exponent as i32 - FLOAT_BIAS - FLOAT_MANTISSA_BITS as i32 - 2,
50            (1u32 << FLOAT_MANTISSA_BITS) | ieee_mantissa,
51        )
52    };
53    let even = (m2 & 1) == 0;
54    let accept_bounds = even;
55
56    // Step 2: Determine the interval of valid decimal representations.
57    let mv = 4 * m2;
58    let mp = 4 * m2 + 2;
59    // Implicit bool -> int conversion. True is 1, false is 0.
60    let mm_shift = (ieee_mantissa != 0 || ieee_exponent <= 1) as u32;
61    let mm = 4 * m2 - 1 - mm_shift;
62
63    // Step 3: Convert to a decimal power base using 64-bit arithmetic.
64    let mut vr: u32;
65    let mut vp: u32;
66    let mut vm: u32;
67    let e10: i32;
68    let mut vm_is_trailing_zeros = false;
69    let mut vr_is_trailing_zeros = false;
70    let mut last_removed_digit = 0u8;
71    if e2 >= 0 {
72        let q = log10_pow2(e2);
73        e10 = q as i32;
74        let k = FLOAT_POW5_INV_BITCOUNT + pow5bits(q as i32) - 1;
75        let i = -e2 + q as i32 + k;
76        vr = mul_pow5_inv_div_pow2(mv, q, i);
77        vp = mul_pow5_inv_div_pow2(mp, q, i);
78        vm = mul_pow5_inv_div_pow2(mm, q, i);
79        if q != 0 && (vp - 1) / 10 <= vm / 10 {
80            // We need to know one removed digit even if we are not going to loop below. We could use
81            // q = X - 1 above, except that would require 33 bits for the result, and we've found that
82            // 32-bit arithmetic is faster even on 64-bit machines.
83            let l = FLOAT_POW5_INV_BITCOUNT + pow5bits(q as i32 - 1) - 1;
84            last_removed_digit =
85                (mul_pow5_inv_div_pow2(mv, q - 1, -e2 + q as i32 - 1 + l) % 10) as u8;
86        }
87        if q <= 9 {
88            // The largest power of 5 that fits in 24 bits is 5^10, but q <= 9 seems to be safe as well.
89            // Only one of mp, mv, and mm can be a multiple of 5, if any.
90            if mv % 5 == 0 {
91                vr_is_trailing_zeros = multiple_of_power_of_5_32(mv, q);
92            } else if accept_bounds {
93                vm_is_trailing_zeros = multiple_of_power_of_5_32(mm, q);
94            } else {
95                vp -= multiple_of_power_of_5_32(mp, q) as u32;
96            }
97        }
98    } else {
99        let q = log10_pow5(-e2);
100        e10 = q as i32 + e2;
101        let i = -e2 - q as i32;
102        let k = pow5bits(i) - FLOAT_POW5_BITCOUNT;
103        let mut j = q as i32 - k;
104        vr = mul_pow5_div_pow2(mv, i as u32, j);
105        vp = mul_pow5_div_pow2(mp, i as u32, j);
106        vm = mul_pow5_div_pow2(mm, i as u32, j);
107        if q != 0 && (vp - 1) / 10 <= vm / 10 {
108            j = q as i32 - 1 - (pow5bits(i + 1) - FLOAT_POW5_BITCOUNT);
109            last_removed_digit = (mul_pow5_div_pow2(mv, (i + 1) as u32, j) % 10) as u8;
110        }
111        if q <= 1 {
112            // {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits.
113            // mv = 4 * m2, so it always has at least two trailing 0 bits.
114            vr_is_trailing_zeros = true;
115            if accept_bounds {
116                // mm = mv - 1 - mm_shift, so it has 1 trailing 0 bit iff mm_shift == 1.
117                vm_is_trailing_zeros = mm_shift == 1;
118            } else {
119                // mp = mv + 2, so it always has at least one trailing 0 bit.
120                vp -= 1;
121            }
122        } else if q < 31 {
123            // TODO(ulfjack): Use a tighter bound here.
124            vr_is_trailing_zeros = multiple_of_power_of_2_32(mv, q - 1);
125        }
126    }
127
128    // Step 4: Find the shortest decimal representation in the interval of valid representations.
129    let mut removed = 0i32;
130    let output = if vm_is_trailing_zeros || vr_is_trailing_zeros {
131        // General case, which happens rarely (~4.0%).
132        while vp / 10 > vm / 10 {
133            vm_is_trailing_zeros &= vm - (vm / 10) * 10 == 0;
134            vr_is_trailing_zeros &= last_removed_digit == 0;
135            last_removed_digit = (vr % 10) as u8;
136            vr /= 10;
137            vp /= 10;
138            vm /= 10;
139            removed += 1;
140        }
141        if vm_is_trailing_zeros {
142            while vm % 10 == 0 {
143                vr_is_trailing_zeros &= last_removed_digit == 0;
144                last_removed_digit = (vr % 10) as u8;
145                vr /= 10;
146                vp /= 10;
147                vm /= 10;
148                removed += 1;
149            }
150        }
151        if vr_is_trailing_zeros && last_removed_digit == 5 && vr % 2 == 0 {
152            // Round even if the exact number is .....50..0.
153            last_removed_digit = 4;
154        }
155        // We need to take vr + 1 if vr is outside bounds or we need to round up.
156        vr + ((vr == vm && (!accept_bounds || !vm_is_trailing_zeros)) || last_removed_digit >= 5)
157            as u32
158    } else {
159        // Specialized for the common case (~96.0%). Percentages below are relative to this.
160        // Loop iterations below (approximately):
161        // 0: 13.6%, 1: 70.7%, 2: 14.1%, 3: 1.39%, 4: 0.14%, 5+: 0.01%
162        while vp / 10 > vm / 10 {
163            last_removed_digit = (vr % 10) as u8;
164            vr /= 10;
165            vp /= 10;
166            vm /= 10;
167            removed += 1;
168        }
169        // We need to take vr + 1 if vr is outside bounds or we need to round up.
170        vr + (vr == vm || last_removed_digit >= 5) as u32
171    };
172    let exp = e10 + removed;
173
174    FloatingDecimal32 {
175        exponent: exp,
176        mantissa: output,
177    }
178}