1 | use crate::common::*; |
---|---|

2 | use crate::d2s; |

3 | use crate::d2s_intrinsics::*; |

4 | use crate::parse::Error; |

5 | #[cfg(feature = "no-panic")] |

6 | use no_panic::no_panic; |

7 | |

8 | const DOUBLE_EXPONENT_BIAS: usize = 1023; |

9 | |

10 | fn floor_log2(value: u64) -> u32 { |

11 | 63_u32.wrapping_sub(value.leading_zeros()) |

12 | } |

13 | |

14 | #[cfg_attr(feature = "no-panic", no_panic)] |

15 | pub fn s2d(buffer: &[u8]) -> Result<f64, Error> { |

16 | let len = buffer.len(); |

17 | if len == 0 { |

18 | return Err(Error::InputTooShort); |

19 | } |

20 | |

21 | let mut m10digits = 0; |

22 | let mut e10digits = 0; |

23 | let mut dot_index = len; |

24 | let mut e_index = len; |

25 | let mut m10 = 0u64; |

26 | let mut e10 = 0i32; |

27 | let mut signed_m = false; |

28 | let mut signed_e = false; |

29 | |

30 | let mut i = 0; |

31 | if unsafe { *buffer.get_unchecked(0) } == b'-'{ |

32 | signed_m = true; |

33 | i += 1; |

34 | } |

35 | |

36 | while let Some(c) = buffer.get(i).copied() { |

37 | if c == b'.'{ |

38 | if dot_index != len { |

39 | return Err(Error::MalformedInput); |

40 | } |

41 | dot_index = i; |

42 | i += 1; |

43 | continue; |

44 | } |

45 | if c < b'0'|| c > b'9'{ |

46 | break; |

47 | } |

48 | if m10digits >= 17 { |

49 | return Err(Error::InputTooLong); |

50 | } |

51 | m10 = 10 * m10 + (c - b'0') as u64; |

52 | if m10 != 0 { |

53 | m10digits += 1; |

54 | } |

55 | i += 1; |

56 | } |

57 | |

58 | if let Some(b'e') | Some( b'E') = buffer.get(i) { |

59 | e_index = i; |

60 | i += 1; |

61 | match buffer.get(i) { |

62 | Some(b'-') => { |

63 | signed_e = true; |

64 | i += 1; |

65 | } |

66 | Some(b'+') => i += 1, |

67 | _ => {} |

68 | } |

69 | while let Some(c) = buffer.get(i).copied() { |

70 | if c < b'0'|| c > b'9'{ |

71 | return Err(Error::MalformedInput); |

72 | } |

73 | if e10digits > 3 { |

74 | // TODO: Be more lenient. Return +/-Infinity or +/-0 instead. |

75 | return Err(Error::InputTooLong); |

76 | } |

77 | e10 = 10 * e10 + (c - b'0') as i32; |

78 | if e10 != 0 { |

79 | e10digits += 1; |

80 | } |

81 | i += 1; |

82 | } |

83 | } |

84 | |

85 | if i < len { |

86 | return Err(Error::MalformedInput); |

87 | } |

88 | if signed_e { |

89 | e10 = -e10; |

90 | } |

91 | e10 -= if dot_index < e_index { |

92 | (e_index - dot_index - 1) as i32 |

93 | } else { |

94 | 0 |

95 | }; |

96 | if m10 == 0 { |

97 | return Ok(if signed_m { -0.0 } else { 0.0 }); |

98 | } |

99 | |

100 | if m10digits + e10 <= -324 || m10 == 0 { |

101 | // Number is less than 1e-324, which should be rounded down to 0; return |

102 | // +/-0.0. |

103 | let ieee = (signed_m as u64) << (d2s::DOUBLE_EXPONENT_BITS + d2s::DOUBLE_MANTISSA_BITS); |

104 | return Ok(f64::from_bits(ieee)); |

105 | } |

106 | if m10digits + e10 >= 310 { |

107 | // Number is larger than 1e+309, which should be rounded to +/-Infinity. |

108 | let ieee = ((signed_m as u64) << (d2s::DOUBLE_EXPONENT_BITS + d2s::DOUBLE_MANTISSA_BITS)) |

109 | | (0x7ff_u64 << d2s::DOUBLE_MANTISSA_BITS); |

110 | return Ok(f64::from_bits(ieee)); |

111 | } |

112 | |

113 | // Convert to binary float m2 * 2^e2, while retaining information about |

114 | // whether the conversion was exact (trailing_zeros). |

115 | let e2: i32; |

116 | let m2: u64; |

117 | let mut trailing_zeros: bool; |

118 | if e10 >= 0 { |

119 | // The length of m * 10^e in bits is: |

120 | // log2(m10 * 10^e10) = log2(m10) + e10 log2(10) = log2(m10) + e10 + e10 * log2(5) |

121 | // |

122 | // We want to compute the DOUBLE_MANTISSA_BITS + 1 top-most bits (+1 for |

123 | // the implicit leading one in IEEE format). We therefore choose a |

124 | // binary output exponent of |

125 | // log2(m10 * 10^e10) - (DOUBLE_MANTISSA_BITS + 1). |

126 | // |

127 | // We use floor(log2(5^e10)) so that we get at least this many bits; |

128 | // better to have an additional bit than to not have enough bits. |

129 | e2 = floor_log2(m10) |

130 | .wrapping_add(e10 as u32) |

131 | .wrapping_add(log2_pow5(e10) as u32) |

132 | .wrapping_sub(d2s::DOUBLE_MANTISSA_BITS + 1) as i32; |

133 | |

134 | // We now compute [m10 * 10^e10 / 2^e2] = [m10 * 5^e10 / 2^(e2-e10)]. |

135 | // To that end, we use the DOUBLE_POW5_SPLIT table. |

136 | let j = e2 |

137 | .wrapping_sub(e10) |

138 | .wrapping_sub(ceil_log2_pow5(e10)) |

139 | .wrapping_add(d2s::DOUBLE_POW5_BITCOUNT); |

140 | debug_assert!(j >= 0); |

141 | debug_assert!(e10 < d2s::DOUBLE_POW5_SPLIT.len() as i32); |

142 | m2 = mul_shift_64( |

143 | m10, |

144 | unsafe { d2s::DOUBLE_POW5_SPLIT.get_unchecked(e10 as usize) }, |

145 | j as u32, |

146 | ); |

147 | |

148 | // We also compute if the result is exact, i.e., |

149 | // [m10 * 10^e10 / 2^e2] == m10 * 10^e10 / 2^e2. |

150 | // This can only be the case if 2^e2 divides m10 * 10^e10, which in turn |

151 | // requires that the largest power of 2 that divides m10 + e10 is |

152 | // greater than e2. If e2 is less than e10, then the result must be |

153 | // exact. Otherwise we use the existing multiple_of_power_of_2 function. |

154 | trailing_zeros = |

155 | e2 < e10 || e2 - e10 < 64 && multiple_of_power_of_2(m10, (e2 - e10) as u32); |

156 | } else { |

157 | e2 = floor_log2(m10) |

158 | .wrapping_add(e10 as u32) |

159 | .wrapping_sub(ceil_log2_pow5(-e10) as u32) |

160 | .wrapping_sub(d2s::DOUBLE_MANTISSA_BITS + 1) as i32; |

161 | let j = e2 |

162 | .wrapping_sub(e10) |

163 | .wrapping_add(ceil_log2_pow5(-e10)) |

164 | .wrapping_sub(1) |

165 | .wrapping_add(d2s::DOUBLE_POW5_INV_BITCOUNT); |

166 | debug_assert!(-e10 < d2s::DOUBLE_POW5_INV_SPLIT.len() as i32); |

167 | m2 = mul_shift_64( |

168 | m10, |

169 | unsafe { d2s::DOUBLE_POW5_INV_SPLIT.get_unchecked(-e10 as usize) }, |

170 | j as u32, |

171 | ); |

172 | trailing_zeros = multiple_of_power_of_5(m10, -e10 as u32); |

173 | } |

174 | |

175 | // Compute the final IEEE exponent. |

176 | let mut ieee_e2 = i32::max(0, e2 + DOUBLE_EXPONENT_BIAS as i32 + floor_log2(m2) as i32) as u32; |

177 | |

178 | if ieee_e2 > 0x7fe { |

179 | // Final IEEE exponent is larger than the maximum representable; return +/-Infinity. |

180 | let ieee = ((signed_m as u64) << (d2s::DOUBLE_EXPONENT_BITS + d2s::DOUBLE_MANTISSA_BITS)) |

181 | | (0x7ff_u64 << d2s::DOUBLE_MANTISSA_BITS); |

182 | return Ok(f64::from_bits(ieee)); |

183 | } |

184 | |

185 | // We need to figure out how much we need to shift m2. The tricky part is |

186 | // that we need to take the final IEEE exponent into account, so we need to |

187 | // reverse the bias and also special-case the value 0. |

188 | let shift = if ieee_e2 == 0 { 1 } else { ieee_e2 as i32 } |

189 | .wrapping_sub(e2) |

190 | .wrapping_sub(DOUBLE_EXPONENT_BIAS as i32) |

191 | .wrapping_sub(d2s::DOUBLE_MANTISSA_BITS as i32); |

192 | debug_assert!(shift >= 0); |

193 | |

194 | // We need to round up if the exact value is more than 0.5 above the value |

195 | // we computed. That's equivalent to checking if the last removed bit was 1 |

196 | // and either the value was not just trailing zeros or the result would |

197 | // otherwise be odd. |

198 | // |

199 | // We need to update trailing_zeros given that we have the exact output |

200 | // exponent ieee_e2 now. |

201 | trailing_zeros &= (m2 & ((1_u64 << (shift - 1)) - 1)) == 0; |

202 | let last_removed_bit = (m2 >> (shift - 1)) & 1; |

203 | let round_up = last_removed_bit != 0 && (!trailing_zeros || ((m2 >> shift) & 1) != 0); |

204 | |

205 | let mut ieee_m2 = (m2 >> shift).wrapping_add(round_up as u64); |

206 | debug_assert!(ieee_m2 <= 1_u64 << (d2s::DOUBLE_MANTISSA_BITS + 1)); |

207 | ieee_m2 &= (1_u64 << d2s::DOUBLE_MANTISSA_BITS) - 1; |

208 | if ieee_m2 == 0 && round_up { |

209 | // Due to how the IEEE represents +/-Infinity, we don't need to check |

210 | // for overflow here. |

211 | ieee_e2 += 1; |

212 | } |

213 | let ieee = ((((signed_m as u64) << d2s::DOUBLE_EXPONENT_BITS) | ieee_e2 as u64) |

214 | << d2s::DOUBLE_MANTISSA_BITS) |

215 | | ieee_m2; |

216 | Ok(f64::from_bits(ieee)) |

217 | } |

218 |