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::{log10_pow2, log10_pow5, pow5bits}; |
22 | use 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 | |
26 | pub const FLOAT_MANTISSA_BITS: u32 = 23; |
27 | pub const FLOAT_EXPONENT_BITS: u32 = 8; |
28 | const FLOAT_BIAS: i32 = 127; |
29 | pub use crate::f2s_intrinsics::{FLOAT_POW5_BITCOUNT, FLOAT_POW5_INV_BITCOUNT}; |
30 | |
31 | // A floating decimal representing m * 10^e. |
32 | pub 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)] |
40 | pub 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 | } |
179 | |