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