Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1 | // After editing this file, run "go generate" in the parent directory. |
| 2 | |
| 3 | // Copyright 2020 The Wuffs Authors. |
| 4 | // |
| 5 | // Licensed under the Apache License, Version 2.0 (the "License"); |
| 6 | // you may not use this file except in compliance with the License. |
| 7 | // You may obtain a copy of the License at |
| 8 | // |
| 9 | // https://www.apache.org/licenses/LICENSE-2.0 |
| 10 | // |
| 11 | // Unless required by applicable law or agreed to in writing, software |
| 12 | // distributed under the License is distributed on an "AS IS" BASIS, |
| 13 | // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 14 | // See the License for the specific language governing permissions and |
| 15 | // limitations under the License. |
| 16 | |
| 17 | // ---------------- IEEE 754 Floating Point |
| 18 | |
Nigel Tao | 7bf7cf2 | 2020-07-12 16:23:15 +1000 | [diff] [blame] | 19 | WUFFS_BASE__MAYBE_STATIC wuffs_base__lossy_value_u16 // |
Nigel Tao | a3931d5 | 2020-07-12 21:06:44 +1000 | [diff] [blame] | 20 | wuffs_base__ieee_754_bit_representation__from_f64_to_u16_truncate(double f) { |
Nigel Tao | 7bf7cf2 | 2020-07-12 16:23:15 +1000 | [diff] [blame] | 21 | uint64_t u = 0; |
| 22 | if (sizeof(uint64_t) == sizeof(double)) { |
| 23 | memcpy(&u, &f, sizeof(uint64_t)); |
| 24 | } |
Nigel Tao | 56d9096 | 2020-07-12 21:11:49 +1000 | [diff] [blame] | 25 | uint16_t neg = ((uint16_t)((u >> 63) << 15)); |
Nigel Tao | 7bf7cf2 | 2020-07-12 16:23:15 +1000 | [diff] [blame] | 26 | u &= 0x7FFFFFFFFFFFFFFF; |
| 27 | uint64_t exp = u >> 52; |
| 28 | uint64_t man = u & 0x000FFFFFFFFFFFFF; |
| 29 | |
| 30 | if (exp == 0x7FF) { |
| 31 | if (man == 0) { // Infinity. |
| 32 | wuffs_base__lossy_value_u16 ret; |
| 33 | ret.value = neg | 0x7C00; |
| 34 | ret.lossy = false; |
| 35 | return ret; |
| 36 | } |
| 37 | // NaN. Shift the 52 mantissa bits to 10 mantissa bits, keeping the most |
| 38 | // significant mantissa bit (quiet vs signaling NaNs). Also set the low 9 |
| 39 | // bits of ret.value so that the 10-bit mantissa is non-zero. |
| 40 | wuffs_base__lossy_value_u16 ret; |
| 41 | ret.value = neg | 0x7DFF | ((uint16_t)(man >> 42)); |
| 42 | ret.lossy = false; |
| 43 | return ret; |
| 44 | |
| 45 | } else if (exp > 0x40E) { // Truncate to the largest finite f16. |
| 46 | wuffs_base__lossy_value_u16 ret; |
| 47 | ret.value = neg | 0x7BFF; |
| 48 | ret.lossy = true; |
| 49 | return ret; |
| 50 | |
| 51 | } else if (exp <= 0x3E6) { // Truncate to zero. |
| 52 | wuffs_base__lossy_value_u16 ret; |
| 53 | ret.value = neg; |
| 54 | ret.lossy = (u != 0); |
| 55 | return ret; |
| 56 | |
| 57 | } else if (exp <= 0x3F0) { // Normal f64, subnormal f16. |
| 58 | // Convert from a 53-bit mantissa (after realizing the implicit bit) to a |
| 59 | // 10-bit mantissa and then adjust for the exponent. |
| 60 | man |= 0x0010000000000000; |
Nigel Tao | 56d9096 | 2020-07-12 21:11:49 +1000 | [diff] [blame] | 61 | uint32_t shift = ((uint32_t)(1051 - exp)); // 1051 = 0x3F0 + 53 - 10. |
Nigel Tao | 7bf7cf2 | 2020-07-12 16:23:15 +1000 | [diff] [blame] | 62 | uint64_t shifted_man = man >> shift; |
| 63 | wuffs_base__lossy_value_u16 ret; |
| 64 | ret.value = neg | ((uint16_t)shifted_man); |
| 65 | ret.lossy = (shifted_man << shift) != man; |
| 66 | return ret; |
| 67 | } |
| 68 | |
| 69 | // Normal f64, normal f16. |
| 70 | |
| 71 | // Re-bias from 1023 to 15 and shift above f16's 10 mantissa bits. |
| 72 | exp = (exp - 1008) << 10; // 1008 = 1023 - 15 = 0x3FF - 0xF. |
| 73 | |
| 74 | // Convert from a 52-bit mantissa (excluding the implicit bit) to a 10-bit |
| 75 | // mantissa (again excluding the implicit bit). We lose some information if |
| 76 | // any of the bottom 42 bits are non-zero. |
| 77 | wuffs_base__lossy_value_u16 ret; |
| 78 | ret.value = neg | ((uint16_t)exp) | ((uint16_t)(man >> 42)); |
| 79 | ret.lossy = (man << 22) != 0; |
| 80 | return ret; |
| 81 | } |
| 82 | |
| 83 | WUFFS_BASE__MAYBE_STATIC wuffs_base__lossy_value_u32 // |
Nigel Tao | a3931d5 | 2020-07-12 21:06:44 +1000 | [diff] [blame] | 84 | wuffs_base__ieee_754_bit_representation__from_f64_to_u32_truncate(double f) { |
Nigel Tao | 7bf7cf2 | 2020-07-12 16:23:15 +1000 | [diff] [blame] | 85 | uint64_t u = 0; |
| 86 | if (sizeof(uint64_t) == sizeof(double)) { |
| 87 | memcpy(&u, &f, sizeof(uint64_t)); |
| 88 | } |
| 89 | uint32_t neg = ((uint32_t)(u >> 63)) << 31; |
| 90 | u &= 0x7FFFFFFFFFFFFFFF; |
| 91 | uint64_t exp = u >> 52; |
| 92 | uint64_t man = u & 0x000FFFFFFFFFFFFF; |
| 93 | |
| 94 | if (exp == 0x7FF) { |
| 95 | if (man == 0) { // Infinity. |
| 96 | wuffs_base__lossy_value_u32 ret; |
| 97 | ret.value = neg | 0x7F800000; |
| 98 | ret.lossy = false; |
| 99 | return ret; |
| 100 | } |
| 101 | // NaN. Shift the 52 mantissa bits to 23 mantissa bits, keeping the most |
| 102 | // significant mantissa bit (quiet vs signaling NaNs). Also set the low 22 |
| 103 | // bits of ret.value so that the 23-bit mantissa is non-zero. |
| 104 | wuffs_base__lossy_value_u32 ret; |
| 105 | ret.value = neg | 0x7FBFFFFF | ((uint32_t)(man >> 29)); |
| 106 | ret.lossy = false; |
| 107 | return ret; |
| 108 | |
| 109 | } else if (exp > 0x47E) { // Truncate to the largest finite f32. |
| 110 | wuffs_base__lossy_value_u32 ret; |
| 111 | ret.value = neg | 0x7F7FFFFF; |
| 112 | ret.lossy = true; |
| 113 | return ret; |
| 114 | |
| 115 | } else if (exp <= 0x369) { // Truncate to zero. |
| 116 | wuffs_base__lossy_value_u32 ret; |
| 117 | ret.value = neg; |
| 118 | ret.lossy = (u != 0); |
| 119 | return ret; |
| 120 | |
| 121 | } else if (exp <= 0x380) { // Normal f64, subnormal f32. |
| 122 | // Convert from a 53-bit mantissa (after realizing the implicit bit) to a |
| 123 | // 23-bit mantissa and then adjust for the exponent. |
| 124 | man |= 0x0010000000000000; |
Nigel Tao | 56d9096 | 2020-07-12 21:11:49 +1000 | [diff] [blame] | 125 | uint32_t shift = ((uint32_t)(926 - exp)); // 926 = 0x380 + 53 - 23. |
Nigel Tao | 7bf7cf2 | 2020-07-12 16:23:15 +1000 | [diff] [blame] | 126 | uint64_t shifted_man = man >> shift; |
| 127 | wuffs_base__lossy_value_u32 ret; |
| 128 | ret.value = neg | ((uint32_t)shifted_man); |
| 129 | ret.lossy = (shifted_man << shift) != man; |
| 130 | return ret; |
| 131 | } |
| 132 | |
| 133 | // Normal f64, normal f32. |
| 134 | |
| 135 | // Re-bias from 1023 to 127 and shift above f32's 23 mantissa bits. |
| 136 | exp = (exp - 896) << 23; // 896 = 1023 - 127 = 0x3FF - 0x7F. |
| 137 | |
| 138 | // Convert from a 52-bit mantissa (excluding the implicit bit) to a 23-bit |
| 139 | // mantissa (again excluding the implicit bit). We lose some information if |
| 140 | // any of the bottom 29 bits are non-zero. |
| 141 | wuffs_base__lossy_value_u32 ret; |
| 142 | ret.value = neg | ((uint32_t)exp) | ((uint32_t)(man >> 29)); |
| 143 | ret.lossy = (man << 35) != 0; |
| 144 | return ret; |
| 145 | } |
| 146 | |
| 147 | // -------- |
| 148 | |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 149 | #define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE 2047 |
| 150 | #define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION 800 |
| 151 | |
| 152 | // WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL is the largest N |
| 153 | // such that ((10 << N) < (1 << 64)). |
| 154 | #define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL 60 |
| 155 | |
| 156 | // wuffs_base__private_implementation__high_prec_dec (abbreviated as HPD) is a |
| 157 | // fixed precision floating point decimal number, augmented with ±infinity |
| 158 | // values, but it cannot represent NaN (Not a Number). |
| 159 | // |
| 160 | // "High precision" means that the mantissa holds 800 decimal digits. 800 is |
| 161 | // WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION. |
| 162 | // |
| 163 | // An HPD isn't for general purpose arithmetic, only for conversions to and |
| 164 | // from IEEE 754 double-precision floating point, where the largest and |
| 165 | // smallest positive, finite values are approximately 1.8e+308 and 4.9e-324. |
| 166 | // HPD exponents above +2047 mean infinity, below -2047 mean zero. The ±2047 |
| 167 | // bounds are further away from zero than ±(324 + 800), where 800 and 2047 is |
| 168 | // WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION and |
| 169 | // WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE. |
| 170 | // |
| 171 | // digits[.. num_digits] are the number's digits in big-endian order. The |
| 172 | // uint8_t values are in the range [0 ..= 9], not ['0' ..= '9'], where e.g. '7' |
| 173 | // is the ASCII value 0x37. |
| 174 | // |
| 175 | // decimal_point is the index (within digits) of the decimal point. It may be |
| 176 | // negative or be larger than num_digits, in which case the explicit digits are |
| 177 | // padded with implicit zeroes. |
| 178 | // |
| 179 | // For example, if num_digits is 3 and digits is "\x07\x08\x09": |
| 180 | // - A decimal_point of -2 means ".00789" |
| 181 | // - A decimal_point of -1 means ".0789" |
| 182 | // - A decimal_point of +0 means ".789" |
| 183 | // - A decimal_point of +1 means "7.89" |
| 184 | // - A decimal_point of +2 means "78.9" |
| 185 | // - A decimal_point of +3 means "789." |
| 186 | // - A decimal_point of +4 means "7890." |
| 187 | // - A decimal_point of +5 means "78900." |
| 188 | // |
| 189 | // As above, a decimal_point higher than +2047 means that the overall value is |
| 190 | // infinity, lower than -2047 means zero. |
| 191 | // |
| 192 | // negative is a sign bit. An HPD can distinguish positive and negative zero. |
| 193 | // |
| 194 | // truncated is whether there are more than |
| 195 | // WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION digits, and at |
| 196 | // least one of those extra digits are non-zero. The existence of long-tail |
| 197 | // digits can affect rounding. |
| 198 | // |
| 199 | // The "all fields are zero" value is valid, and represents the number +0. |
| 200 | typedef struct { |
| 201 | uint32_t num_digits; |
| 202 | int32_t decimal_point; |
| 203 | bool negative; |
| 204 | bool truncated; |
| 205 | uint8_t digits[WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION]; |
| 206 | } wuffs_base__private_implementation__high_prec_dec; |
| 207 | |
| 208 | // wuffs_base__private_implementation__high_prec_dec__trim trims trailing |
| 209 | // zeroes from the h->digits[.. h->num_digits] slice. They have no benefit, |
| 210 | // since we explicitly track h->decimal_point. |
| 211 | // |
| 212 | // Preconditions: |
| 213 | // - h is non-NULL. |
| 214 | static inline void // |
| 215 | wuffs_base__private_implementation__high_prec_dec__trim( |
| 216 | wuffs_base__private_implementation__high_prec_dec* h) { |
| 217 | while ((h->num_digits > 0) && (h->digits[h->num_digits - 1] == 0)) { |
| 218 | h->num_digits--; |
| 219 | } |
| 220 | } |
| 221 | |
| 222 | // wuffs_base__private_implementation__high_prec_dec__assign sets h to |
| 223 | // represent the number x. |
| 224 | // |
| 225 | // Preconditions: |
| 226 | // - h is non-NULL. |
| 227 | static void // |
| 228 | wuffs_base__private_implementation__high_prec_dec__assign( |
| 229 | wuffs_base__private_implementation__high_prec_dec* h, |
| 230 | uint64_t x, |
| 231 | bool negative) { |
| 232 | uint32_t n = 0; |
| 233 | |
| 234 | // Set h->digits. |
| 235 | if (x > 0) { |
| 236 | // Calculate the digits, working right-to-left. After we determine n (how |
| 237 | // many digits there are), copy from buf to h->digits. |
| 238 | // |
| 239 | // UINT64_MAX, 18446744073709551615, is 20 digits long. It can be faster to |
| 240 | // copy a constant number of bytes than a variable number (20 instead of |
| 241 | // n). Make buf large enough (and start writing to it from the middle) so |
| 242 | // that can we always copy 20 bytes: the slice buf[(20-n) .. (40-n)]. |
| 243 | uint8_t buf[40] = {0}; |
| 244 | uint8_t* ptr = &buf[20]; |
| 245 | do { |
| 246 | uint64_t remaining = x / 10; |
| 247 | x -= remaining * 10; |
| 248 | ptr--; |
| 249 | *ptr = (uint8_t)x; |
| 250 | n++; |
| 251 | x = remaining; |
| 252 | } while (x > 0); |
| 253 | memcpy(h->digits, ptr, 20); |
| 254 | } |
| 255 | |
| 256 | // Set h's other fields. |
| 257 | h->num_digits = n; |
| 258 | h->decimal_point = (int32_t)n; |
| 259 | h->negative = negative; |
| 260 | h->truncated = false; |
| 261 | wuffs_base__private_implementation__high_prec_dec__trim(h); |
| 262 | } |
| 263 | |
| 264 | static wuffs_base__status // |
| 265 | wuffs_base__private_implementation__high_prec_dec__parse( |
| 266 | wuffs_base__private_implementation__high_prec_dec* h, |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 267 | wuffs_base__slice_u8 s, |
| 268 | uint32_t options) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 269 | if (!h) { |
| 270 | return wuffs_base__make_status(wuffs_base__error__bad_receiver); |
| 271 | } |
| 272 | h->num_digits = 0; |
| 273 | h->decimal_point = 0; |
| 274 | h->negative = false; |
| 275 | h->truncated = false; |
| 276 | |
| 277 | uint8_t* p = s.ptr; |
| 278 | uint8_t* q = s.ptr + s.len; |
| 279 | |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 280 | if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) { |
| 281 | for (;; p++) { |
| 282 | if (p >= q) { |
| 283 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 284 | } else if (*p != '_') { |
| 285 | break; |
| 286 | } |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 287 | } |
| 288 | } |
| 289 | |
| 290 | // Parse sign. |
| 291 | do { |
| 292 | if (*p == '+') { |
| 293 | p++; |
| 294 | } else if (*p == '-') { |
| 295 | h->negative = true; |
| 296 | p++; |
| 297 | } else { |
| 298 | break; |
| 299 | } |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 300 | if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) { |
| 301 | for (;; p++) { |
| 302 | if (p >= q) { |
| 303 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 304 | } else if (*p != '_') { |
| 305 | break; |
| 306 | } |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 307 | } |
| 308 | } |
| 309 | } while (0); |
| 310 | |
| 311 | // Parse digits, up to (and including) a '.', 'E' or 'e'. Examples for each |
| 312 | // limb in this if-else chain: |
| 313 | // - "0.789" |
| 314 | // - "1002.789" |
| 315 | // - ".789" |
| 316 | // - Other (invalid input). |
| 317 | uint32_t nd = 0; |
| 318 | int32_t dp = 0; |
| 319 | bool no_digits_before_separator = false; |
Nigel Tao | e82bc8e | 2020-07-11 12:49:15 +1000 | [diff] [blame] | 320 | if (('0' == *p) && |
| 321 | !(options & |
| 322 | WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_MULTIPLE_LEADING_ZEROES)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 323 | p++; |
| 324 | for (;; p++) { |
| 325 | if (p >= q) { |
| 326 | goto after_all; |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 327 | } else if (*p == |
| 328 | ((options & |
| 329 | WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA) |
| 330 | ? ',' |
| 331 | : '.')) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 332 | p++; |
| 333 | goto after_sep; |
| 334 | } else if ((*p == 'E') || (*p == 'e')) { |
| 335 | p++; |
| 336 | goto after_exp; |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 337 | } else if ((*p != '_') || |
| 338 | !(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 339 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 340 | } |
| 341 | } |
| 342 | |
Nigel Tao | e82bc8e | 2020-07-11 12:49:15 +1000 | [diff] [blame] | 343 | } else if (('0' <= *p) && (*p <= '9')) { |
| 344 | if (*p == '0') { |
| 345 | for (; (p < q) && (*p == '0'); p++) { |
| 346 | } |
| 347 | } else { |
| 348 | h->digits[nd++] = (uint8_t)(*p - '0'); |
| 349 | dp = (int32_t)nd; |
| 350 | p++; |
| 351 | } |
| 352 | |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 353 | for (;; p++) { |
| 354 | if (p >= q) { |
| 355 | goto after_all; |
| 356 | } else if (('0' <= *p) && (*p <= '9')) { |
| 357 | if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 358 | h->digits[nd++] = (uint8_t)(*p - '0'); |
| 359 | dp = (int32_t)nd; |
| 360 | } else if ('0' != *p) { |
| 361 | // Long-tail non-zeroes set the truncated bit. |
| 362 | h->truncated = true; |
| 363 | } |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 364 | } else if (*p == |
| 365 | ((options & |
| 366 | WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA) |
| 367 | ? ',' |
| 368 | : '.')) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 369 | p++; |
| 370 | goto after_sep; |
| 371 | } else if ((*p == 'E') || (*p == 'e')) { |
| 372 | p++; |
| 373 | goto after_exp; |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 374 | } else if ((*p != '_') || |
| 375 | !(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 376 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 377 | } |
| 378 | } |
| 379 | |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 380 | } else if (*p == ((options & |
| 381 | WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA) |
| 382 | ? ',' |
| 383 | : '.')) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 384 | p++; |
| 385 | no_digits_before_separator = true; |
| 386 | |
| 387 | } else { |
| 388 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 389 | } |
| 390 | |
| 391 | after_sep: |
| 392 | for (;; p++) { |
| 393 | if (p >= q) { |
| 394 | goto after_all; |
| 395 | } else if ('0' == *p) { |
| 396 | if (nd == 0) { |
| 397 | // Track leading zeroes implicitly. |
| 398 | dp--; |
| 399 | } else if (nd < |
| 400 | WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 401 | h->digits[nd++] = (uint8_t)(*p - '0'); |
| 402 | } |
| 403 | } else if (('0' < *p) && (*p <= '9')) { |
| 404 | if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 405 | h->digits[nd++] = (uint8_t)(*p - '0'); |
| 406 | } else { |
| 407 | // Long-tail non-zeroes set the truncated bit. |
| 408 | h->truncated = true; |
| 409 | } |
| 410 | } else if ((*p == 'E') || (*p == 'e')) { |
| 411 | p++; |
| 412 | goto after_exp; |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 413 | } else if ((*p != '_') || |
| 414 | !(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 415 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 416 | } |
| 417 | } |
| 418 | |
| 419 | after_exp: |
| 420 | do { |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 421 | if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) { |
| 422 | for (;; p++) { |
| 423 | if (p >= q) { |
| 424 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 425 | } else if (*p != '_') { |
| 426 | break; |
| 427 | } |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 428 | } |
| 429 | } |
| 430 | |
| 431 | int32_t exp_sign = +1; |
| 432 | if (*p == '+') { |
| 433 | p++; |
| 434 | } else if (*p == '-') { |
| 435 | exp_sign = -1; |
| 436 | p++; |
| 437 | } |
| 438 | |
| 439 | int32_t exp = 0; |
| 440 | const int32_t exp_large = |
| 441 | WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE + |
| 442 | WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION; |
| 443 | bool saw_exp_digits = false; |
| 444 | for (; p < q; p++) { |
Nigel Tao | c5c9885 | 2020-07-11 13:10:14 +1000 | [diff] [blame] | 445 | if ((*p == '_') && |
| 446 | (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 447 | // No-op. |
| 448 | } else if (('0' <= *p) && (*p <= '9')) { |
| 449 | saw_exp_digits = true; |
| 450 | if (exp < exp_large) { |
| 451 | exp = (10 * exp) + ((int32_t)(*p - '0')); |
| 452 | } |
| 453 | } else { |
| 454 | break; |
| 455 | } |
| 456 | } |
| 457 | if (!saw_exp_digits) { |
| 458 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 459 | } |
| 460 | dp += exp_sign * exp; |
| 461 | } while (0); |
| 462 | |
| 463 | after_all: |
| 464 | if (p != q) { |
| 465 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 466 | } |
| 467 | h->num_digits = nd; |
| 468 | if (nd == 0) { |
| 469 | if (no_digits_before_separator) { |
| 470 | return wuffs_base__make_status(wuffs_base__error__bad_argument); |
| 471 | } |
| 472 | h->decimal_point = 0; |
| 473 | } else if (dp < |
| 474 | -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) { |
| 475 | h->decimal_point = |
| 476 | -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE - 1; |
| 477 | } else if (dp > |
| 478 | +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) { |
| 479 | h->decimal_point = |
| 480 | +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE + 1; |
| 481 | } else { |
| 482 | h->decimal_point = dp; |
| 483 | } |
| 484 | wuffs_base__private_implementation__high_prec_dec__trim(h); |
| 485 | return wuffs_base__make_status(NULL); |
| 486 | } |
| 487 | |
| 488 | // -------- |
| 489 | |
| 490 | // wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits |
| 491 | // returns the number of additional decimal digits when left-shifting by shift. |
| 492 | // |
| 493 | // See below for preconditions. |
| 494 | static uint32_t // |
| 495 | wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits( |
| 496 | wuffs_base__private_implementation__high_prec_dec* h, |
| 497 | uint32_t shift) { |
| 498 | // Masking with 0x3F should be unnecessary (assuming the preconditions) but |
| 499 | // it's cheap and ensures that we don't overflow the |
| 500 | // wuffs_base__private_implementation__hpd_left_shift array. |
| 501 | shift &= 63; |
| 502 | |
| 503 | uint32_t x_a = wuffs_base__private_implementation__hpd_left_shift[shift]; |
| 504 | uint32_t x_b = wuffs_base__private_implementation__hpd_left_shift[shift + 1]; |
| 505 | uint32_t num_new_digits = x_a >> 11; |
| 506 | uint32_t pow5_a = 0x7FF & x_a; |
| 507 | uint32_t pow5_b = 0x7FF & x_b; |
| 508 | |
| 509 | const uint8_t* pow5 = |
| 510 | &wuffs_base__private_implementation__powers_of_5[pow5_a]; |
| 511 | uint32_t i = 0; |
| 512 | uint32_t n = pow5_b - pow5_a; |
| 513 | for (; i < n; i++) { |
| 514 | if (i >= h->num_digits) { |
| 515 | return num_new_digits - 1; |
| 516 | } else if (h->digits[i] == pow5[i]) { |
| 517 | continue; |
| 518 | } else if (h->digits[i] < pow5[i]) { |
| 519 | return num_new_digits - 1; |
| 520 | } else { |
| 521 | return num_new_digits; |
| 522 | } |
| 523 | } |
| 524 | return num_new_digits; |
| 525 | } |
| 526 | |
| 527 | // -------- |
| 528 | |
| 529 | // wuffs_base__private_implementation__high_prec_dec__rounded_integer returns |
| 530 | // the integral (non-fractional) part of h, provided that it is 18 or fewer |
| 531 | // decimal digits. For 19 or more digits, it returns UINT64_MAX. Note that: |
| 532 | // - (1 << 53) is 9007199254740992, which has 16 decimal digits. |
| 533 | // - (1 << 56) is 72057594037927936, which has 17 decimal digits. |
| 534 | // - (1 << 59) is 576460752303423488, which has 18 decimal digits. |
| 535 | // - (1 << 63) is 9223372036854775808, which has 19 decimal digits. |
| 536 | // and that IEEE 754 double precision has 52 mantissa bits. |
| 537 | // |
| 538 | // That integral part is rounded-to-even: rounding 7.5 or 8.5 both give 8. |
| 539 | // |
| 540 | // h's negative bit is ignored: rounding -8.6 returns 9. |
| 541 | // |
| 542 | // See below for preconditions. |
| 543 | static uint64_t // |
| 544 | wuffs_base__private_implementation__high_prec_dec__rounded_integer( |
| 545 | wuffs_base__private_implementation__high_prec_dec* h) { |
| 546 | if ((h->num_digits == 0) || (h->decimal_point < 0)) { |
| 547 | return 0; |
| 548 | } else if (h->decimal_point > 18) { |
| 549 | return UINT64_MAX; |
| 550 | } |
| 551 | |
| 552 | uint32_t dp = (uint32_t)(h->decimal_point); |
| 553 | uint64_t n = 0; |
| 554 | uint32_t i = 0; |
| 555 | for (; i < dp; i++) { |
| 556 | n = (10 * n) + ((i < h->num_digits) ? h->digits[i] : 0); |
| 557 | } |
| 558 | |
| 559 | bool round_up = false; |
| 560 | if (dp < h->num_digits) { |
| 561 | round_up = h->digits[dp] >= 5; |
| 562 | if ((h->digits[dp] == 5) && (dp + 1 == h->num_digits)) { |
| 563 | // We are exactly halfway. If we're truncated, round up, otherwise round |
| 564 | // to even. |
| 565 | round_up = h->truncated || // |
| 566 | ((dp > 0) && (1 & h->digits[dp - 1])); |
| 567 | } |
| 568 | } |
| 569 | if (round_up) { |
| 570 | n++; |
| 571 | } |
| 572 | |
| 573 | return n; |
| 574 | } |
| 575 | |
| 576 | // wuffs_base__private_implementation__high_prec_dec__small_xshift shifts h's |
| 577 | // number (where 'x' is 'l' or 'r' for left or right) by a small shift value. |
| 578 | // |
| 579 | // Preconditions: |
| 580 | // - h is non-NULL. |
| 581 | // - h->decimal_point is "not extreme". |
| 582 | // - shift is non-zero. |
| 583 | // - shift is "a small shift". |
| 584 | // |
| 585 | // "Not extreme" means within |
| 586 | // ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE. |
| 587 | // |
| 588 | // "A small shift" means not more than |
| 589 | // WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL. |
| 590 | // |
| 591 | // wuffs_base__private_implementation__high_prec_dec__rounded_integer and |
| 592 | // wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits |
| 593 | // have the same preconditions. |
| 594 | // |
| 595 | // wuffs_base__private_implementation__high_prec_dec__lshift keeps the first |
| 596 | // two preconditions but not the last two. Its shift argument is signed and |
| 597 | // does not need to be "small": zero is a no-op, positive means left shift and |
| 598 | // negative means right shift. |
| 599 | |
| 600 | static void // |
| 601 | wuffs_base__private_implementation__high_prec_dec__small_lshift( |
| 602 | wuffs_base__private_implementation__high_prec_dec* h, |
| 603 | uint32_t shift) { |
| 604 | if (h->num_digits == 0) { |
| 605 | return; |
| 606 | } |
| 607 | uint32_t num_new_digits = |
| 608 | wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits( |
| 609 | h, shift); |
| 610 | uint32_t rx = h->num_digits - 1; // Read index. |
| 611 | uint32_t wx = h->num_digits - 1 + num_new_digits; // Write index. |
| 612 | uint64_t n = 0; |
| 613 | |
| 614 | // Repeat: pick up a digit, put down a digit, right to left. |
| 615 | while (((int32_t)rx) >= 0) { |
| 616 | n += ((uint64_t)(h->digits[rx])) << shift; |
| 617 | uint64_t quo = n / 10; |
| 618 | uint64_t rem = n - (10 * quo); |
| 619 | if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 620 | h->digits[wx] = (uint8_t)rem; |
| 621 | } else if (rem > 0) { |
| 622 | h->truncated = true; |
| 623 | } |
| 624 | n = quo; |
| 625 | wx--; |
| 626 | rx--; |
| 627 | } |
| 628 | |
| 629 | // Put down leading digits, right to left. |
| 630 | while (n > 0) { |
| 631 | uint64_t quo = n / 10; |
| 632 | uint64_t rem = n - (10 * quo); |
| 633 | if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 634 | h->digits[wx] = (uint8_t)rem; |
| 635 | } else if (rem > 0) { |
| 636 | h->truncated = true; |
| 637 | } |
| 638 | n = quo; |
| 639 | wx--; |
| 640 | } |
| 641 | |
| 642 | // Finish. |
| 643 | h->num_digits += num_new_digits; |
| 644 | if (h->num_digits > |
| 645 | WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 646 | h->num_digits = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION; |
| 647 | } |
| 648 | h->decimal_point += (int32_t)num_new_digits; |
| 649 | wuffs_base__private_implementation__high_prec_dec__trim(h); |
| 650 | } |
| 651 | |
| 652 | static void // |
| 653 | wuffs_base__private_implementation__high_prec_dec__small_rshift( |
| 654 | wuffs_base__private_implementation__high_prec_dec* h, |
| 655 | uint32_t shift) { |
| 656 | uint32_t rx = 0; // Read index. |
| 657 | uint32_t wx = 0; // Write index. |
| 658 | uint64_t n = 0; |
| 659 | |
| 660 | // Pick up enough leading digits to cover the first shift. |
| 661 | while ((n >> shift) == 0) { |
| 662 | if (rx < h->num_digits) { |
| 663 | // Read a digit. |
| 664 | n = (10 * n) + h->digits[rx++]; |
| 665 | } else if (n == 0) { |
| 666 | // h's number used to be zero and remains zero. |
| 667 | return; |
| 668 | } else { |
| 669 | // Read sufficient implicit trailing zeroes. |
| 670 | while ((n >> shift) == 0) { |
| 671 | n = 10 * n; |
| 672 | rx++; |
| 673 | } |
| 674 | break; |
| 675 | } |
| 676 | } |
| 677 | h->decimal_point -= ((int32_t)(rx - 1)); |
| 678 | if (h->decimal_point < |
| 679 | -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) { |
| 680 | // After the shift, h's number is effectively zero. |
| 681 | h->num_digits = 0; |
| 682 | h->decimal_point = 0; |
| 683 | h->negative = false; |
| 684 | h->truncated = false; |
| 685 | return; |
| 686 | } |
| 687 | |
| 688 | // Repeat: pick up a digit, put down a digit, left to right. |
| 689 | uint64_t mask = (((uint64_t)(1)) << shift) - 1; |
| 690 | while (rx < h->num_digits) { |
| 691 | uint8_t new_digit = ((uint8_t)(n >> shift)); |
| 692 | n = (10 * (n & mask)) + h->digits[rx++]; |
| 693 | h->digits[wx++] = new_digit; |
| 694 | } |
| 695 | |
| 696 | // Put down trailing digits, left to right. |
| 697 | while (n > 0) { |
| 698 | uint8_t new_digit = ((uint8_t)(n >> shift)); |
| 699 | n = 10 * (n & mask); |
| 700 | if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) { |
| 701 | h->digits[wx++] = new_digit; |
| 702 | } else if (new_digit > 0) { |
| 703 | h->truncated = true; |
| 704 | } |
| 705 | } |
| 706 | |
| 707 | // Finish. |
| 708 | h->num_digits = wx; |
| 709 | wuffs_base__private_implementation__high_prec_dec__trim(h); |
| 710 | } |
| 711 | |
| 712 | static void // |
| 713 | wuffs_base__private_implementation__high_prec_dec__lshift( |
| 714 | wuffs_base__private_implementation__high_prec_dec* h, |
| 715 | int32_t shift) { |
| 716 | if (shift > 0) { |
| 717 | while (shift > +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) { |
| 718 | wuffs_base__private_implementation__high_prec_dec__small_lshift( |
| 719 | h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL); |
| 720 | shift -= WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL; |
| 721 | } |
| 722 | wuffs_base__private_implementation__high_prec_dec__small_lshift( |
| 723 | h, ((uint32_t)(+shift))); |
| 724 | } else if (shift < 0) { |
| 725 | while (shift < -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) { |
| 726 | wuffs_base__private_implementation__high_prec_dec__small_rshift( |
| 727 | h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL); |
| 728 | shift += WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL; |
| 729 | } |
| 730 | wuffs_base__private_implementation__high_prec_dec__small_rshift( |
| 731 | h, ((uint32_t)(-shift))); |
| 732 | } |
| 733 | } |
| 734 | |
| 735 | // -------- |
| 736 | |
| 737 | // wuffs_base__private_implementation__high_prec_dec__round_etc rounds h's |
| 738 | // number. For those functions that take an n argument, rounding produces at |
| 739 | // most n digits (which is not necessarily at most n decimal places). Negative |
| 740 | // n values are ignored, as well as any n greater than or equal to h's number |
| 741 | // of digits. The etc__round_just_enough function implicitly chooses an n to |
| 742 | // implement WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION. |
| 743 | // |
| 744 | // Preconditions: |
| 745 | // - h is non-NULL. |
| 746 | // - h->decimal_point is "not extreme". |
| 747 | // |
| 748 | // "Not extreme" means within |
| 749 | // ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE. |
| 750 | |
| 751 | static void // |
| 752 | wuffs_base__private_implementation__high_prec_dec__round_down( |
| 753 | wuffs_base__private_implementation__high_prec_dec* h, |
| 754 | int32_t n) { |
| 755 | if ((n < 0) || (h->num_digits <= (uint32_t)n)) { |
| 756 | return; |
| 757 | } |
| 758 | h->num_digits = (uint32_t)(n); |
| 759 | wuffs_base__private_implementation__high_prec_dec__trim(h); |
| 760 | } |
| 761 | |
| 762 | static void // |
| 763 | wuffs_base__private_implementation__high_prec_dec__round_up( |
| 764 | wuffs_base__private_implementation__high_prec_dec* h, |
| 765 | int32_t n) { |
| 766 | if ((n < 0) || (h->num_digits <= (uint32_t)n)) { |
| 767 | return; |
| 768 | } |
| 769 | |
| 770 | for (n--; n >= 0; n--) { |
| 771 | if (h->digits[n] < 9) { |
| 772 | h->digits[n]++; |
| 773 | h->num_digits = (uint32_t)(n + 1); |
| 774 | return; |
| 775 | } |
| 776 | } |
| 777 | |
| 778 | // The number is all 9s. Change to a single 1 and adjust the decimal point. |
| 779 | h->digits[0] = 1; |
| 780 | h->num_digits = 1; |
| 781 | h->decimal_point++; |
| 782 | } |
| 783 | |
| 784 | static void // |
| 785 | wuffs_base__private_implementation__high_prec_dec__round_nearest( |
| 786 | wuffs_base__private_implementation__high_prec_dec* h, |
| 787 | int32_t n) { |
| 788 | if ((n < 0) || (h->num_digits <= (uint32_t)n)) { |
| 789 | return; |
| 790 | } |
| 791 | bool up = h->digits[n] >= 5; |
| 792 | if ((h->digits[n] == 5) && ((n + 1) == ((int32_t)(h->num_digits)))) { |
| 793 | up = h->truncated || // |
| 794 | ((n > 0) && ((h->digits[n - 1] & 1) != 0)); |
| 795 | } |
| 796 | |
| 797 | if (up) { |
| 798 | wuffs_base__private_implementation__high_prec_dec__round_up(h, n); |
| 799 | } else { |
| 800 | wuffs_base__private_implementation__high_prec_dec__round_down(h, n); |
| 801 | } |
| 802 | } |
| 803 | |
| 804 | static void // |
| 805 | wuffs_base__private_implementation__high_prec_dec__round_just_enough( |
| 806 | wuffs_base__private_implementation__high_prec_dec* h, |
| 807 | int32_t exp2, |
| 808 | uint64_t mantissa) { |
| 809 | // The magic numbers 52 and 53 in this function are because IEEE 754 double |
| 810 | // precision has 52 mantissa bits. |
| 811 | // |
| 812 | // Let f be the floating point number represented by exp2 and mantissa (and |
| 813 | // also the number in h): the number (mantissa * (2 ** (exp2 - 52))). |
| 814 | // |
| 815 | // If f is zero or a small integer, we can return early. |
| 816 | if ((mantissa == 0) || |
| 817 | ((exp2 < 53) && (h->decimal_point >= ((int32_t)(h->num_digits))))) { |
| 818 | return; |
| 819 | } |
| 820 | |
| 821 | // The smallest normal f has an exp2 of -1022 and a mantissa of (1 << 52). |
| 822 | // Subnormal numbers have the same exp2 but a smaller mantissa. |
| 823 | static const int32_t min_incl_normal_exp2 = -1022; |
| 824 | static const uint64_t min_incl_normal_mantissa = 0x0010000000000000ul; |
| 825 | |
| 826 | // Compute lower and upper bounds such that any number between them (possibly |
| 827 | // inclusive) will round to f. First, the lower bound. Our number f is: |
| 828 | // ((mantissa + 0) * (2 ** ( exp2 - 52))) |
| 829 | // |
| 830 | // The next lowest floating point number is: |
| 831 | // ((mantissa - 1) * (2 ** ( exp2 - 52))) |
| 832 | // unless (mantissa - 1) drops the (1 << 52) bit and exp2 is not the |
| 833 | // min_incl_normal_exp2. Either way, call it: |
| 834 | // ((l_mantissa) * (2 ** (l_exp2 - 52))) |
| 835 | // |
| 836 | // The lower bound is halfway between them (noting that 52 became 53): |
| 837 | // (((2 * l_mantissa) + 1) * (2 ** (l_exp2 - 53))) |
| 838 | int32_t l_exp2 = exp2; |
| 839 | uint64_t l_mantissa = mantissa - 1; |
| 840 | if ((exp2 > min_incl_normal_exp2) && (mantissa <= min_incl_normal_mantissa)) { |
| 841 | l_exp2 = exp2 - 1; |
| 842 | l_mantissa = (2 * mantissa) - 1; |
| 843 | } |
| 844 | wuffs_base__private_implementation__high_prec_dec lower; |
| 845 | wuffs_base__private_implementation__high_prec_dec__assign( |
| 846 | &lower, (2 * l_mantissa) + 1, false); |
| 847 | wuffs_base__private_implementation__high_prec_dec__lshift(&lower, |
| 848 | l_exp2 - 53); |
| 849 | |
| 850 | // Next, the upper bound. Our number f is: |
| 851 | // ((mantissa + 0) * (2 ** (exp2 - 52))) |
| 852 | // |
| 853 | // The next highest floating point number is: |
| 854 | // ((mantissa + 1) * (2 ** (exp2 - 52))) |
| 855 | // |
| 856 | // The upper bound is halfway between them (noting that 52 became 53): |
| 857 | // (((2 * mantissa) + 1) * (2 ** (exp2 - 53))) |
| 858 | wuffs_base__private_implementation__high_prec_dec upper; |
| 859 | wuffs_base__private_implementation__high_prec_dec__assign( |
| 860 | &upper, (2 * mantissa) + 1, false); |
| 861 | wuffs_base__private_implementation__high_prec_dec__lshift(&upper, exp2 - 53); |
| 862 | |
| 863 | // The lower and upper bounds are possible outputs only if the original |
| 864 | // mantissa is even, so that IEEE round-to-even would round to the original |
| 865 | // mantissa and not its neighbors. |
| 866 | bool inclusive = (mantissa & 1) == 0; |
| 867 | |
| 868 | // As we walk the digits, we want to know whether rounding up would fall |
| 869 | // within the upper bound. This is tracked by upper_delta: |
| 870 | // - When -1, the digits of h and upper are the same so far. |
| 871 | // - When +0, we saw a difference of 1 between h and upper on a previous |
| 872 | // digit and subsequently only 9s for h and 0s for upper. Thus, rounding |
| 873 | // up may fall outside of the bound if !inclusive. |
| 874 | // - When +1, the difference is greater than 1 and we know that rounding up |
| 875 | // falls within the bound. |
| 876 | // |
| 877 | // This is a state machine with three states. The numerical value for each |
| 878 | // state (-1, +0 or +1) isn't important, other than their order. |
| 879 | int upper_delta = -1; |
| 880 | |
| 881 | // We can now figure out the shortest number of digits required. Walk the |
| 882 | // digits until h has distinguished itself from lower or upper. |
| 883 | // |
| 884 | // The zi and zd variables are indexes and digits, for z in l (lower), h (the |
| 885 | // number) and u (upper). |
| 886 | // |
| 887 | // The lower, h and upper numbers may have their decimal points at different |
| 888 | // places. In this case, upper is the longest, so we iterate ui starting from |
| 889 | // 0 and iterate li and hi starting from either 0 or -1. |
| 890 | int32_t ui = 0; |
| 891 | for (;; ui++) { |
| 892 | // Calculate hd, the middle number's digit. |
| 893 | int32_t hi = ui - upper.decimal_point + h->decimal_point; |
| 894 | if (hi >= ((int32_t)(h->num_digits))) { |
| 895 | break; |
| 896 | } |
| 897 | uint8_t hd = (((uint32_t)hi) < h->num_digits) ? h->digits[hi] : 0; |
| 898 | |
| 899 | // Calculate ld, the lower bound's digit. |
| 900 | int32_t li = ui - upper.decimal_point + lower.decimal_point; |
| 901 | uint8_t ld = (((uint32_t)li) < lower.num_digits) ? lower.digits[li] : 0; |
| 902 | |
| 903 | // We can round down (truncate) if lower has a different digit than h or if |
| 904 | // lower is inclusive and is exactly the result of rounding down (i.e. we |
| 905 | // have reached the final digit of lower). |
| 906 | bool can_round_down = |
| 907 | (ld != hd) || // |
| 908 | (inclusive && ((li + 1) == ((int32_t)(lower.num_digits)))); |
| 909 | |
| 910 | // Calculate ud, the upper bound's digit, and update upper_delta. |
| 911 | uint8_t ud = (((uint32_t)ui) < upper.num_digits) ? upper.digits[ui] : 0; |
| 912 | if (upper_delta < 0) { |
| 913 | if ((hd + 1) < ud) { |
| 914 | // For example: |
| 915 | // h = 12345??? |
| 916 | // upper = 12347??? |
| 917 | upper_delta = +1; |
| 918 | } else if (hd != ud) { |
| 919 | // For example: |
| 920 | // h = 12345??? |
| 921 | // upper = 12346??? |
| 922 | upper_delta = +0; |
| 923 | } |
| 924 | } else if (upper_delta == 0) { |
| 925 | if ((hd != 9) || (ud != 0)) { |
| 926 | // For example: |
| 927 | // h = 1234598? |
| 928 | // upper = 1234600? |
| 929 | upper_delta = +1; |
| 930 | } |
| 931 | } |
| 932 | |
| 933 | // We can round up if upper has a different digit than h and either upper |
| 934 | // is inclusive or upper is bigger than the result of rounding up. |
| 935 | bool can_round_up = |
| 936 | (upper_delta > 0) || // |
| 937 | ((upper_delta == 0) && // |
| 938 | (inclusive || ((ui + 1) < ((int32_t)(upper.num_digits))))); |
| 939 | |
| 940 | // If we can round either way, round to nearest. If we can round only one |
| 941 | // way, do it. If we can't round, continue the loop. |
| 942 | if (can_round_down) { |
| 943 | if (can_round_up) { |
| 944 | wuffs_base__private_implementation__high_prec_dec__round_nearest( |
| 945 | h, hi + 1); |
| 946 | return; |
| 947 | } else { |
| 948 | wuffs_base__private_implementation__high_prec_dec__round_down(h, |
| 949 | hi + 1); |
| 950 | return; |
| 951 | } |
| 952 | } else { |
| 953 | if (can_round_up) { |
| 954 | wuffs_base__private_implementation__high_prec_dec__round_up(h, hi + 1); |
| 955 | return; |
| 956 | } |
| 957 | } |
| 958 | } |
| 959 | } |
| 960 | |
| 961 | // -------- |
| 962 | |
| 963 | // wuffs_base__private_implementation__parse_number_f64_eisel produces the IEEE |
Nigel Tao | b15a0fc | 2020-07-08 10:50:14 +1000 | [diff] [blame] | 964 | // 754 double-precision value for an exact mantissa and base-10 exponent. For |
| 965 | // example: |
| 966 | // - when parsing "12345.678e+02", man is 12345678 and exp10 is -1. |
| 967 | // - when parsing "-12", man is 12 and exp10 is 0. Processing the leading |
| 968 | // minus sign is the responsibility of the caller, not this function. |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 969 | // |
| 970 | // On success, it returns a non-negative int64_t such that the low 63 bits hold |
| 971 | // the 11-bit exponent and 52-bit mantissa. |
| 972 | // |
| 973 | // On failure, it returns a negative value. |
| 974 | // |
| 975 | // The algorithm is based on an original idea by Michael Eisel. See |
| 976 | // https://lemire.me/blog/2020/03/10/fast-float-parsing-in-practice/ |
| 977 | // |
| 978 | // Preconditions: |
| 979 | // - man is non-zero. |
| 980 | // - exp10 is in the range -326 ..= 310, the same range of the |
| 981 | // wuffs_base__private_implementation__powers_of_10 array. |
| 982 | static int64_t // |
| 983 | wuffs_base__private_implementation__parse_number_f64_eisel(uint64_t man, |
| 984 | int32_t exp10) { |
| 985 | // Look up the (possibly truncated) base-2 representation of (10 ** exp10). |
| 986 | // The look-up table was constructed so that it is already normalized: the |
| 987 | // table entry's mantissa's MSB (most significant bit) is on. |
| 988 | const uint32_t* po10 = |
| 989 | &wuffs_base__private_implementation__powers_of_10[5 * (exp10 + 326)]; |
| 990 | |
| 991 | // Normalize the man argument. The (man != 0) precondition means that a |
| 992 | // non-zero bit exists. |
| 993 | uint32_t clz = wuffs_base__count_leading_zeroes_u64(man); |
| 994 | man <<= clz; |
| 995 | |
| 996 | // Calculate the return value's base-2 exponent. We might tweak it by ±1 |
| 997 | // later, but its initial value comes from the look-up table and clz. |
| 998 | uint64_t ret_exp2 = ((uint64_t)po10[4]) - ((uint64_t)clz); |
| 999 | |
| 1000 | // Multiply the two mantissas. Normalization means that both mantissas are at |
| 1001 | // least (1<<63), so the 128-bit product must be at least (1<<126). The high |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1002 | // 64 bits of the product, x_hi, must therefore be at least (1<<62). |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1003 | // |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1004 | // As a consequence, x_hi has either 0 or 1 leading zeroes. Shifting x_hi |
| 1005 | // right by either 9 or 10 bits (depending on x_hi's MSB) will therefore |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1006 | // leave the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. |
Nigel Tao | e722bb0 | 2020-07-10 20:33:09 +1000 | [diff] [blame] | 1007 | #if defined(__SIZEOF_INT128__) |
Nigel Tao | bff01d7 | 2020-07-10 12:15:09 +1000 | [diff] [blame] | 1008 | // See commit 18449ad75d582dd015c236abc85a16f333b796f3 "Optimize 128-bit muls |
| 1009 | // in parse_number_f64_eisel" for benchmark numbers. |
Nigel Tao | 18449ad | 2020-07-10 11:48:52 +1000 | [diff] [blame] | 1010 | __uint128_t x = |
| 1011 | ((__uint128_t)man) * (((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32)); |
| 1012 | uint64_t x_hi = ((uint64_t)(x >> 64)); |
| 1013 | uint64_t x_lo = ((uint64_t)(x)); |
| 1014 | #else |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1015 | wuffs_base__multiply_u64__output x = wuffs_base__multiply_u64( |
| 1016 | man, ((uint64_t)po10[2]) | (((uint64_t)po10[3]) << 32)); |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1017 | uint64_t x_hi = x.hi; |
| 1018 | uint64_t x_lo = x.lo; |
Nigel Tao | 18449ad | 2020-07-10 11:48:52 +1000 | [diff] [blame] | 1019 | #endif |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1020 | |
| 1021 | // Before we shift right by at least 9 bits, recall that the look-up table |
| 1022 | // entry was possibly truncated. We have so far only calculated a lower bound |
| 1023 | // for the product (man * e), where e is (10 ** exp10). The upper bound would |
| 1024 | // add a further (man * 1) to the 128-bit product, which overflows the lower |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1025 | // 64-bit limb if ((x_lo + man) < man). |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1026 | // |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1027 | // If overflow occurs, that adds 1 to x_hi. Since we're about to shift right |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1028 | // by at least 9 bits, that carried 1 can be ignored unless the higher 64-bit |
| 1029 | // limb's low 9 bits are all on. |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1030 | if (((x_hi & 0x1FF) == 0x1FF) && ((x_lo + man) < man)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1031 | // Refine our calculation of (man * e). Before, our approximation of e used |
| 1032 | // a "low resolution" 64-bit mantissa. Now use a "high resolution" 128-bit |
| 1033 | // mantissa. We've already calculated x = (man * bits_0_to_63_incl_of_e). |
| 1034 | // Now calculate y = (man * bits_64_to_127_incl_of_e). |
Nigel Tao | e722bb0 | 2020-07-10 20:33:09 +1000 | [diff] [blame] | 1035 | #if defined(__SIZEOF_INT128__) |
Nigel Tao | bff01d7 | 2020-07-10 12:15:09 +1000 | [diff] [blame] | 1036 | // See commit 18449ad75d582dd015c236abc85a16f333b796f3 "Optimize 128-bit |
| 1037 | // muls in parse_number_f64_eisel" for benchmark numbers. |
Nigel Tao | 18449ad | 2020-07-10 11:48:52 +1000 | [diff] [blame] | 1038 | __uint128_t y = ((__uint128_t)man) * |
| 1039 | (((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32)); |
| 1040 | uint64_t y_hi = ((uint64_t)(y >> 64)); |
| 1041 | uint64_t y_lo = ((uint64_t)(y)); |
| 1042 | #else |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1043 | wuffs_base__multiply_u64__output y = wuffs_base__multiply_u64( |
| 1044 | man, ((uint64_t)po10[0]) | (((uint64_t)po10[1]) << 32)); |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1045 | uint64_t y_hi = y.hi; |
| 1046 | uint64_t y_lo = y.lo; |
Nigel Tao | 18449ad | 2020-07-10 11:48:52 +1000 | [diff] [blame] | 1047 | #endif |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1048 | |
| 1049 | // Merge the 128-bit x and 128-bit y, which overlap by 64 bits, to |
| 1050 | // calculate the 192-bit product of the 64-bit man by the 128-bit e. |
| 1051 | // As we exit this if-block, we only care about the high 128 bits |
| 1052 | // (merged_hi and merged_lo) of that 192-bit product. |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1053 | uint64_t merged_hi = x_hi; |
| 1054 | uint64_t merged_lo = x_lo + y_hi; |
| 1055 | if (merged_lo < x_lo) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1056 | merged_hi++; // Carry the overflow bit. |
| 1057 | } |
| 1058 | |
| 1059 | // The "high resolution" approximation of e is still a lower bound. Once |
| 1060 | // again, see if the upper bound is large enough to produce a different |
| 1061 | // result. This time, if it does, give up instead of reaching for an even |
| 1062 | // more precise approximation to e. |
| 1063 | // |
| 1064 | // This three-part check is similar to the two-part check that guarded the |
| 1065 | // if block that we're now in, but it has an extra term for the middle 64 |
| 1066 | // bits (checking that adding 1 to merged_lo would overflow). |
| 1067 | if (((merged_hi & 0x1FF) == 0x1FF) && ((merged_lo + 1) == 0) && |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1068 | (y_lo + man < man)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1069 | return -1; |
| 1070 | } |
| 1071 | |
| 1072 | // Replace the 128-bit x with merged. |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1073 | x_hi = merged_hi; |
| 1074 | x_lo = merged_lo; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1075 | } |
| 1076 | |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1077 | // As mentioned above, shifting x_hi right by either 9 or 10 bits will leave |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1078 | // the top 10 MSBs (bits 54 ..= 63) off and the 11th MSB (bit 53) on. If the |
| 1079 | // MSB (before shifting) was on, adjust ret_exp2 for the larger shift. |
| 1080 | // |
| 1081 | // Having bit 53 on (and higher bits off) means that ret_mantissa is a 54-bit |
| 1082 | // number. |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1083 | uint64_t msb = x_hi >> 63; |
| 1084 | uint64_t ret_mantissa = x_hi >> (msb + 9); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1085 | ret_exp2 -= 1 ^ msb; |
| 1086 | |
| 1087 | // IEEE 754 rounds to-nearest with ties rounded to-even. Rounding to-even can |
| 1088 | // be tricky. If we're half-way between two exactly representable numbers |
| 1089 | // (x's low 73 bits are zero and the next 2 bits that matter are "01"), give |
| 1090 | // up instead of trying to pick the winner. |
| 1091 | // |
| 1092 | // Technically, we could tighten the condition by changing "73" to "73 or 74, |
| 1093 | // depending on msb", but a flat "73" is simpler. |
Nigel Tao | 74d4af6 | 2020-07-10 11:27:17 +1000 | [diff] [blame] | 1094 | if ((x_lo == 0) && ((x_hi & 0x1FF) == 0) && ((ret_mantissa & 3) == 1)) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1095 | return -1; |
| 1096 | } |
| 1097 | |
| 1098 | // If we're not halfway then it's rounding to-nearest. Starting with a 54-bit |
| 1099 | // number, carry the lowest bit (bit 0) up if it's on. Regardless of whether |
| 1100 | // it was on or off, shifting right by one then produces a 53-bit number. If |
| 1101 | // carrying up overflowed, shift again. |
| 1102 | ret_mantissa += ret_mantissa & 1; |
| 1103 | ret_mantissa >>= 1; |
| 1104 | if ((ret_mantissa >> 53) > 0) { |
| 1105 | ret_mantissa >>= 1; |
| 1106 | ret_exp2++; |
| 1107 | } |
| 1108 | |
| 1109 | // Starting with a 53-bit number, IEEE 754 double-precision normal numbers |
| 1110 | // have an implicit mantissa bit. Mask that away and keep the low 52 bits. |
| 1111 | ret_mantissa &= 0x000FFFFFFFFFFFFF; |
| 1112 | |
| 1113 | // IEEE 754 double-precision floating point has 11 exponent bits. All off (0) |
| 1114 | // means subnormal numbers. All on (2047) means infinity or NaN. |
| 1115 | if ((ret_exp2 <= 0) || (2047 <= ret_exp2)) { |
| 1116 | return -1; |
| 1117 | } |
| 1118 | |
| 1119 | // Pack the bits and return. |
| 1120 | return ((int64_t)(ret_mantissa | (ret_exp2 << 52))); |
| 1121 | } |
| 1122 | |
| 1123 | // -------- |
| 1124 | |
| 1125 | static wuffs_base__result_f64 // |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 1126 | wuffs_base__private_implementation__parse_number_f64_special( |
| 1127 | wuffs_base__slice_u8 s, |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1128 | uint32_t options) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1129 | do { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1130 | if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) { |
| 1131 | goto fail; |
| 1132 | } |
| 1133 | |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1134 | uint8_t* p = s.ptr; |
| 1135 | uint8_t* q = s.ptr + s.len; |
| 1136 | |
| 1137 | for (; (p < q) && (*p == '_'); p++) { |
| 1138 | } |
| 1139 | if (p >= q) { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1140 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1141 | } |
| 1142 | |
| 1143 | // Parse sign. |
| 1144 | bool negative = false; |
| 1145 | do { |
| 1146 | if (*p == '+') { |
| 1147 | p++; |
| 1148 | } else if (*p == '-') { |
| 1149 | negative = true; |
| 1150 | p++; |
| 1151 | } else { |
| 1152 | break; |
| 1153 | } |
| 1154 | for (; (p < q) && (*p == '_'); p++) { |
| 1155 | } |
| 1156 | } while (0); |
| 1157 | if (p >= q) { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1158 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1159 | } |
| 1160 | |
| 1161 | bool nan = false; |
| 1162 | switch (p[0]) { |
| 1163 | case 'I': |
| 1164 | case 'i': |
| 1165 | if (((q - p) < 3) || // |
| 1166 | ((p[1] != 'N') && (p[1] != 'n')) || // |
| 1167 | ((p[2] != 'F') && (p[2] != 'f'))) { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1168 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1169 | } |
| 1170 | p += 3; |
| 1171 | |
| 1172 | if ((p >= q) || (*p == '_')) { |
| 1173 | break; |
| 1174 | } else if (((q - p) < 5) || // |
| 1175 | ((p[0] != 'I') && (p[0] != 'i')) || // |
| 1176 | ((p[1] != 'N') && (p[1] != 'n')) || // |
| 1177 | ((p[2] != 'I') && (p[2] != 'i')) || // |
| 1178 | ((p[3] != 'T') && (p[3] != 't')) || // |
| 1179 | ((p[4] != 'Y') && (p[4] != 'y'))) { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1180 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1181 | } |
| 1182 | p += 5; |
| 1183 | |
| 1184 | if ((p >= q) || (*p == '_')) { |
| 1185 | break; |
| 1186 | } |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1187 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1188 | |
| 1189 | case 'N': |
| 1190 | case 'n': |
| 1191 | if (((q - p) < 3) || // |
| 1192 | ((p[1] != 'A') && (p[1] != 'a')) || // |
| 1193 | ((p[2] != 'N') && (p[2] != 'n'))) { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1194 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1195 | } |
| 1196 | p += 3; |
| 1197 | |
| 1198 | if ((p >= q) || (*p == '_')) { |
| 1199 | nan = true; |
| 1200 | break; |
| 1201 | } |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1202 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1203 | |
| 1204 | default: |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1205 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1206 | } |
| 1207 | |
| 1208 | // Finish. |
| 1209 | for (; (p < q) && (*p == '_'); p++) { |
| 1210 | } |
| 1211 | if (p != q) { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1212 | goto fail; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1213 | } |
| 1214 | wuffs_base__result_f64 ret; |
| 1215 | ret.status.repr = NULL; |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1216 | ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64( |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1217 | (nan ? 0x7FFFFFFFFFFFFFFF : 0x7FF0000000000000) | |
| 1218 | (negative ? 0x8000000000000000 : 0)); |
| 1219 | return ret; |
| 1220 | } while (0); |
| 1221 | |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1222 | fail: |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1223 | do { |
| 1224 | wuffs_base__result_f64 ret; |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1225 | ret.status.repr = wuffs_base__error__bad_argument; |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1226 | ret.value = 0; |
| 1227 | return ret; |
| 1228 | } while (0); |
| 1229 | } |
| 1230 | |
| 1231 | WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 // |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 1232 | wuffs_base__private_implementation__high_prec_dec__to_f64( |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1233 | wuffs_base__private_implementation__high_prec_dec* h, |
| 1234 | uint32_t options) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1235 | do { |
| 1236 | // powers converts decimal powers of 10 to binary powers of 2. For example, |
| 1237 | // (10000 >> 13) is 1. It stops before the elements exceed 60, also known |
| 1238 | // as WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL. |
| 1239 | static const uint32_t num_powers = 19; |
| 1240 | static const uint8_t powers[19] = { |
| 1241 | 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, // |
| 1242 | 33, 36, 39, 43, 46, 49, 53, 56, 59, // |
| 1243 | }; |
| 1244 | |
| 1245 | // Handle zero and obvious extremes. The largest and smallest positive |
| 1246 | // finite f64 values are approximately 1.8e+308 and 4.9e-324. |
| 1247 | if ((h->num_digits == 0) || (h->decimal_point < -326)) { |
| 1248 | goto zero; |
| 1249 | } else if (h->decimal_point > 310) { |
| 1250 | goto infinity; |
| 1251 | } |
| 1252 | |
| 1253 | // Try the fast Eisel algorithm again. Calculating the (man, exp10) pair |
| 1254 | // from the high_prec_dec h is more correct but slower than the approach |
| 1255 | // taken in wuffs_base__parse_number_f64. The latter is optimized for the |
| 1256 | // common cases (e.g. assuming no underscores or a leading '+' sign) rather |
| 1257 | // than the full set of cases allowed by the Wuffs API. |
| 1258 | if (h->num_digits <= 19) { |
| 1259 | uint64_t man = 0; |
| 1260 | uint32_t i; |
| 1261 | for (i = 0; i < h->num_digits; i++) { |
| 1262 | man = (10 * man) + h->digits[i]; |
| 1263 | } |
| 1264 | int32_t exp10 = h->decimal_point - ((int32_t)(h->num_digits)); |
| 1265 | if ((man != 0) && (-326 <= exp10) && (exp10 <= 310)) { |
| 1266 | int64_t r = wuffs_base__private_implementation__parse_number_f64_eisel( |
| 1267 | man, exp10); |
| 1268 | if (r >= 0) { |
| 1269 | wuffs_base__result_f64 ret; |
| 1270 | ret.status.repr = NULL; |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1271 | ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64( |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1272 | ((uint64_t)r) | (((uint64_t)(h->negative)) << 63)); |
| 1273 | return ret; |
| 1274 | } |
| 1275 | } |
| 1276 | } |
| 1277 | |
| 1278 | // Scale by powers of 2 until we're in the range [½ .. 1], which gives us |
| 1279 | // our exponent (in base-2). First we shift right, possibly a little too |
| 1280 | // far, ending with a value certainly below 1 and possibly below ½... |
| 1281 | const int32_t f64_bias = -1023; |
| 1282 | int32_t exp2 = 0; |
| 1283 | while (h->decimal_point > 0) { |
| 1284 | uint32_t n = (uint32_t)(+h->decimal_point); |
| 1285 | uint32_t shift = |
| 1286 | (n < num_powers) |
| 1287 | ? powers[n] |
| 1288 | : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL; |
| 1289 | |
| 1290 | wuffs_base__private_implementation__high_prec_dec__small_rshift(h, shift); |
| 1291 | if (h->decimal_point < |
| 1292 | -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) { |
| 1293 | goto zero; |
| 1294 | } |
| 1295 | exp2 += (int32_t)shift; |
| 1296 | } |
| 1297 | // ...then we shift left, putting us in [½ .. 1]. |
| 1298 | while (h->decimal_point <= 0) { |
| 1299 | uint32_t shift; |
| 1300 | if (h->decimal_point == 0) { |
| 1301 | if (h->digits[0] >= 5) { |
| 1302 | break; |
| 1303 | } |
| 1304 | shift = (h->digits[0] <= 2) ? 2 : 1; |
| 1305 | } else { |
| 1306 | uint32_t n = (uint32_t)(-h->decimal_point); |
| 1307 | shift = (n < num_powers) |
| 1308 | ? powers[n] |
| 1309 | : WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL; |
| 1310 | } |
| 1311 | |
| 1312 | wuffs_base__private_implementation__high_prec_dec__small_lshift(h, shift); |
| 1313 | if (h->decimal_point > |
| 1314 | +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) { |
| 1315 | goto infinity; |
| 1316 | } |
| 1317 | exp2 -= (int32_t)shift; |
| 1318 | } |
| 1319 | |
| 1320 | // We're in the range [½ .. 1] but f64 uses [1 .. 2]. |
| 1321 | exp2--; |
| 1322 | |
| 1323 | // The minimum normal exponent is (f64_bias + 1). |
| 1324 | while ((f64_bias + 1) > exp2) { |
| 1325 | uint32_t n = (uint32_t)((f64_bias + 1) - exp2); |
| 1326 | if (n > WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) { |
| 1327 | n = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL; |
| 1328 | } |
| 1329 | wuffs_base__private_implementation__high_prec_dec__small_rshift(h, n); |
| 1330 | exp2 += (int32_t)n; |
| 1331 | } |
| 1332 | |
| 1333 | // Check for overflow. |
| 1334 | if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1. |
| 1335 | goto infinity; |
| 1336 | } |
| 1337 | |
| 1338 | // Extract 53 bits for the mantissa (in base-2). |
| 1339 | wuffs_base__private_implementation__high_prec_dec__small_lshift(h, 53); |
| 1340 | uint64_t man2 = |
| 1341 | wuffs_base__private_implementation__high_prec_dec__rounded_integer(h); |
| 1342 | |
| 1343 | // Rounding might have added one bit. If so, shift and re-check overflow. |
| 1344 | if ((man2 >> 53) != 0) { |
| 1345 | man2 >>= 1; |
| 1346 | exp2++; |
| 1347 | if ((exp2 - f64_bias) >= 0x07FF) { // (1 << 11) - 1. |
| 1348 | goto infinity; |
| 1349 | } |
| 1350 | } |
| 1351 | |
| 1352 | // Handle subnormal numbers. |
| 1353 | if ((man2 >> 52) == 0) { |
| 1354 | exp2 = f64_bias; |
| 1355 | } |
| 1356 | |
| 1357 | // Pack the bits and return. |
| 1358 | uint64_t exp2_bits = |
| 1359 | (uint64_t)((exp2 - f64_bias) & 0x07FF); // (1 << 11) - 1. |
| 1360 | uint64_t bits = (man2 & 0x000FFFFFFFFFFFFF) | // (1 << 52) - 1. |
| 1361 | (exp2_bits << 52) | // |
| 1362 | (h->negative ? 0x8000000000000000 : 0); // (1 << 63). |
| 1363 | |
| 1364 | wuffs_base__result_f64 ret; |
| 1365 | ret.status.repr = NULL; |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1366 | ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1367 | return ret; |
| 1368 | } while (0); |
| 1369 | |
| 1370 | zero: |
| 1371 | do { |
| 1372 | uint64_t bits = h->negative ? 0x8000000000000000 : 0; |
| 1373 | |
| 1374 | wuffs_base__result_f64 ret; |
| 1375 | ret.status.repr = NULL; |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1376 | ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1377 | return ret; |
| 1378 | } while (0); |
| 1379 | |
| 1380 | infinity: |
| 1381 | do { |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1382 | if (options & WUFFS_BASE__PARSE_NUMBER_FXX__REJECT_INF_AND_NAN) { |
| 1383 | wuffs_base__result_f64 ret; |
| 1384 | ret.status.repr = wuffs_base__error__bad_argument; |
| 1385 | ret.value = 0; |
| 1386 | return ret; |
| 1387 | } |
| 1388 | |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1389 | uint64_t bits = h->negative ? 0xFFF0000000000000 : 0x7FF0000000000000; |
| 1390 | |
| 1391 | wuffs_base__result_f64 ret; |
| 1392 | ret.status.repr = NULL; |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1393 | ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64(bits); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1394 | return ret; |
| 1395 | } while (0); |
| 1396 | } |
| 1397 | |
| 1398 | static inline bool // |
| 1399 | wuffs_base__private_implementation__is_decimal_digit(uint8_t c) { |
| 1400 | return ('0' <= c) && (c <= '9'); |
| 1401 | } |
| 1402 | |
| 1403 | WUFFS_BASE__MAYBE_STATIC wuffs_base__result_f64 // |
| 1404 | wuffs_base__parse_number_f64(wuffs_base__slice_u8 s, uint32_t options) { |
| 1405 | // In practice, almost all "dd.ddddE±xxx" numbers can be represented |
| 1406 | // losslessly by a uint64_t mantissa "dddddd" and an int32_t base-10 |
| 1407 | // exponent, adjusting "xxx" for the position (if present) of the decimal |
| 1408 | // separator '.' or ','. |
| 1409 | // |
| 1410 | // This (u64 man, i32 exp10) data structure is superficially similar to the |
| 1411 | // "Do It Yourself Floating Point" type from Loitsch (†), but the exponent |
| 1412 | // here is base-10, not base-2. |
| 1413 | // |
| 1414 | // If s's number fits in a (man, exp10), parse that pair with the Eisel |
| 1415 | // algorithm. If not, or if Eisel fails, parsing s with the fallback |
| 1416 | // algorithm is slower but comprehensive. |
| 1417 | // |
| 1418 | // † "Printing Floating-Point Numbers Quickly and Accurately with Integers" |
| 1419 | // (https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf). |
| 1420 | // Florian Loitsch is also the primary contributor to |
| 1421 | // https://github.com/google/double-conversion |
| 1422 | do { |
| 1423 | // Calculating that (man, exp10) pair needs to stay within s's bounds. |
| 1424 | // Provided that s isn't extremely long, work on a NUL-terminated copy of |
| 1425 | // s's contents. The NUL byte isn't a valid part of "±dd.ddddE±xxx". |
| 1426 | // |
| 1427 | // As the pointer p walks the contents, it's faster to repeatedly check "is |
| 1428 | // *p a valid digit" than "is p within bounds and *p a valid digit". |
| 1429 | if (s.len >= 256) { |
| 1430 | goto fallback; |
| 1431 | } |
| 1432 | uint8_t z[256]; |
| 1433 | memcpy(&z[0], s.ptr, s.len); |
| 1434 | z[s.len] = 0; |
| 1435 | const uint8_t* p = &z[0]; |
| 1436 | |
| 1437 | // Look for a leading minus sign. Technically, we could also look for an |
| 1438 | // optional plus sign, but the "script/process-json-numbers.c with -p" |
| 1439 | // benchmark is noticably slower if we do. It's optional and, in practice, |
| 1440 | // usually absent. Let the fallback catch it. |
| 1441 | bool negative = (*p == '-'); |
| 1442 | if (negative) { |
| 1443 | p++; |
| 1444 | } |
| 1445 | |
| 1446 | // After walking "dd.dddd", comparing p later with p now will produce the |
| 1447 | // number of "d"s and "."s. |
| 1448 | const uint8_t* const start_of_digits_ptr = p; |
| 1449 | |
| 1450 | // Walk the "d"s before a '.', 'E', NUL byte, etc. If it starts with '0', |
| 1451 | // it must be a single '0'. If it starts with a non-zero decimal digit, it |
| 1452 | // can be a sequence of decimal digits. |
| 1453 | // |
| 1454 | // Update the man variable during the walk. It's OK if man overflows now. |
| 1455 | // We'll detect that later. |
| 1456 | uint64_t man; |
| 1457 | if (*p == '0') { |
| 1458 | man = 0; |
| 1459 | p++; |
| 1460 | if (wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1461 | goto fallback; |
| 1462 | } |
| 1463 | } else if (wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1464 | man = ((uint8_t)(*p - '0')); |
| 1465 | p++; |
| 1466 | for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) { |
| 1467 | man = (10 * man) + ((uint8_t)(*p - '0')); |
| 1468 | } |
| 1469 | } else { |
| 1470 | goto fallback; |
| 1471 | } |
| 1472 | |
| 1473 | // Walk the "d"s after the optional decimal separator ('.' or ','), |
| 1474 | // updating the man and exp10 variables. |
| 1475 | int32_t exp10 = 0; |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 1476 | if (*p == |
| 1477 | ((options & WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA) |
| 1478 | ? ',' |
| 1479 | : '.')) { |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1480 | p++; |
| 1481 | const uint8_t* first_after_separator_ptr = p; |
| 1482 | if (!wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1483 | goto fallback; |
| 1484 | } |
| 1485 | man = (10 * man) + ((uint8_t)(*p - '0')); |
| 1486 | p++; |
| 1487 | for (; wuffs_base__private_implementation__is_decimal_digit(*p); p++) { |
| 1488 | man = (10 * man) + ((uint8_t)(*p - '0')); |
| 1489 | } |
| 1490 | exp10 = ((int32_t)(first_after_separator_ptr - p)); |
| 1491 | } |
| 1492 | |
| 1493 | // Count the number of digits: |
| 1494 | // - for an input of "314159", digit_count is 6. |
| 1495 | // - for an input of "3.14159", digit_count is 7. |
| 1496 | // |
| 1497 | // This is off-by-one if there is a decimal separator. That's OK for now. |
| 1498 | // We'll correct for that later. The "script/process-json-numbers.c with |
| 1499 | // -p" benchmark is noticably slower if we try to correct for that now. |
| 1500 | uint32_t digit_count = (uint32_t)(p - start_of_digits_ptr); |
| 1501 | |
| 1502 | // Update exp10 for the optional exponent, starting with 'E' or 'e'. |
| 1503 | if ((*p | 0x20) == 'e') { |
| 1504 | p++; |
| 1505 | int32_t exp_sign = +1; |
| 1506 | if (*p == '-') { |
| 1507 | p++; |
| 1508 | exp_sign = -1; |
| 1509 | } else if (*p == '+') { |
| 1510 | p++; |
| 1511 | } |
| 1512 | if (!wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1513 | goto fallback; |
| 1514 | } |
| 1515 | int32_t exp_num = ((uint8_t)(*p - '0')); |
| 1516 | p++; |
| 1517 | // The rest of the exp_num walking has a peculiar control flow but, once |
| 1518 | // again, the "script/process-json-numbers.c with -p" benchmark is |
| 1519 | // sensitive to alternative formulations. |
| 1520 | if (wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1521 | exp_num = (10 * exp_num) + ((uint8_t)(*p - '0')); |
| 1522 | p++; |
| 1523 | } |
| 1524 | if (wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1525 | exp_num = (10 * exp_num) + ((uint8_t)(*p - '0')); |
| 1526 | p++; |
| 1527 | } |
| 1528 | while (wuffs_base__private_implementation__is_decimal_digit(*p)) { |
| 1529 | if (exp_num > 0x1000000) { |
| 1530 | goto fallback; |
| 1531 | } |
| 1532 | exp_num = (10 * exp_num) + ((uint8_t)(*p - '0')); |
| 1533 | p++; |
| 1534 | } |
| 1535 | exp10 += exp_sign * exp_num; |
| 1536 | } |
| 1537 | |
| 1538 | // The Wuffs API is that the original slice has no trailing data. It also |
| 1539 | // allows underscores, which we don't catch here but the fallback should. |
| 1540 | if (p != &z[s.len]) { |
| 1541 | goto fallback; |
| 1542 | } |
| 1543 | |
| 1544 | // Check that the uint64_t typed man variable has not overflowed, based on |
| 1545 | // digit_count. |
| 1546 | // |
| 1547 | // For reference: |
| 1548 | // - (1 << 63) is 9223372036854775808, which has 19 decimal digits. |
| 1549 | // - (1 << 64) is 18446744073709551616, which has 20 decimal digits. |
| 1550 | // - 19 nines, 9999999999999999999, is 0x8AC7230489E7FFFF, which has 64 |
| 1551 | // bits and 16 hexadecimal digits. |
| 1552 | // - 20 nines, 99999999999999999999, is 0x56BC75E2D630FFFFF, which has 67 |
| 1553 | // bits and 17 hexadecimal digits. |
| 1554 | if (digit_count > 19) { |
| 1555 | // Even if we have more than 19 pseudo-digits, it's not yet definitely an |
| 1556 | // overflow. Recall that digit_count might be off-by-one (too large) if |
| 1557 | // there's a decimal separator. It will also over-report the number of |
| 1558 | // meaningful digits if the input looks something like "0.000dddExxx". |
| 1559 | // |
| 1560 | // We adjust by the number of leading '0's and '.'s and re-compare to 19. |
| 1561 | // Once again, technically, we could skip ','s too, but that perturbs the |
| 1562 | // "script/process-json-numbers.c with -p" benchmark. |
| 1563 | const uint8_t* q = start_of_digits_ptr; |
| 1564 | for (; (*q == '0') || (*q == '.'); q++) { |
| 1565 | } |
| 1566 | digit_count -= (uint32_t)(q - start_of_digits_ptr); |
| 1567 | if (digit_count > 19) { |
| 1568 | goto fallback; |
| 1569 | } |
| 1570 | } |
| 1571 | |
| 1572 | // The wuffs_base__private_implementation__parse_number_f64_eisel |
| 1573 | // preconditions include that exp10 is in the range -326 ..= 310. |
| 1574 | if ((exp10 < -326) || (310 < exp10)) { |
| 1575 | goto fallback; |
| 1576 | } |
| 1577 | |
| 1578 | // If man and exp10 are small enough, all three of (man), (10 ** exp10) and |
| 1579 | // (man ** (10 ** exp10)) are exactly representable by a double. We don't |
| 1580 | // need to run the Eisel algorithm. |
| 1581 | if ((-22 <= exp10) && (exp10 <= 22) && ((man >> 53) == 0)) { |
| 1582 | double d = (double)man; |
| 1583 | if (exp10 >= 0) { |
| 1584 | d *= wuffs_base__private_implementation__f64_powers_of_10[+exp10]; |
| 1585 | } else { |
| 1586 | d /= wuffs_base__private_implementation__f64_powers_of_10[-exp10]; |
| 1587 | } |
| 1588 | wuffs_base__result_f64 ret; |
| 1589 | ret.status.repr = NULL; |
| 1590 | ret.value = negative ? -d : +d; |
| 1591 | return ret; |
| 1592 | } |
| 1593 | |
| 1594 | // The wuffs_base__private_implementation__parse_number_f64_eisel |
| 1595 | // preconditions include that man is non-zero. Parsing "0" should be caught |
| 1596 | // by the "If man and exp10 are small enough" above, but "0e99" might not. |
| 1597 | if (man == 0) { |
| 1598 | goto fallback; |
| 1599 | } |
| 1600 | |
| 1601 | // Our man and exp10 are in range. Run the Eisel algorithm. |
| 1602 | int64_t r = |
| 1603 | wuffs_base__private_implementation__parse_number_f64_eisel(man, exp10); |
| 1604 | if (r < 0) { |
| 1605 | goto fallback; |
| 1606 | } |
| 1607 | wuffs_base__result_f64 ret; |
| 1608 | ret.status.repr = NULL; |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1609 | ret.value = wuffs_base__ieee_754_bit_representation__from_u64_to_f64( |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1610 | ((uint64_t)r) | (((uint64_t)negative) << 63)); |
| 1611 | return ret; |
| 1612 | } while (0); |
| 1613 | |
| 1614 | fallback: |
| 1615 | do { |
| 1616 | wuffs_base__private_implementation__high_prec_dec h; |
| 1617 | wuffs_base__status status = |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 1618 | wuffs_base__private_implementation__high_prec_dec__parse(&h, s, |
| 1619 | options); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1620 | if (status.repr) { |
Nigel Tao | e0c5de9 | 2020-07-11 11:48:17 +1000 | [diff] [blame] | 1621 | return wuffs_base__private_implementation__parse_number_f64_special( |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1622 | s, options); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1623 | } |
Nigel Tao | 4d61a05 | 2020-07-11 12:34:40 +1000 | [diff] [blame] | 1624 | return wuffs_base__private_implementation__high_prec_dec__to_f64(&h, |
| 1625 | options); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1626 | } while (0); |
| 1627 | } |
| 1628 | |
| 1629 | // -------- |
| 1630 | |
| 1631 | static inline size_t // |
| 1632 | wuffs_base__private_implementation__render_inf(wuffs_base__slice_u8 dst, |
| 1633 | bool neg, |
| 1634 | uint32_t options) { |
| 1635 | if (neg) { |
| 1636 | if (dst.len < 4) { |
| 1637 | return 0; |
| 1638 | } |
| 1639 | wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492D); // '-Inf'le. |
| 1640 | return 4; |
| 1641 | } |
| 1642 | |
| 1643 | if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) { |
| 1644 | if (dst.len < 4) { |
| 1645 | return 0; |
| 1646 | } |
| 1647 | wuffs_base__store_u32le__no_bounds_check(dst.ptr, 0x666E492B); // '+Inf'le. |
| 1648 | return 4; |
| 1649 | } |
| 1650 | |
| 1651 | if (dst.len < 3) { |
| 1652 | return 0; |
| 1653 | } |
| 1654 | wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x666E49); // 'Inf'le. |
| 1655 | return 3; |
| 1656 | } |
| 1657 | |
| 1658 | static inline size_t // |
| 1659 | wuffs_base__private_implementation__render_nan(wuffs_base__slice_u8 dst) { |
| 1660 | if (dst.len < 3) { |
| 1661 | return 0; |
| 1662 | } |
| 1663 | wuffs_base__store_u24le__no_bounds_check(dst.ptr, 0x4E614E); // 'NaN'le. |
| 1664 | return 3; |
| 1665 | } |
| 1666 | |
| 1667 | static size_t // |
| 1668 | wuffs_base__private_implementation__high_prec_dec__render_exponent_absent( |
| 1669 | wuffs_base__slice_u8 dst, |
| 1670 | wuffs_base__private_implementation__high_prec_dec* h, |
| 1671 | uint32_t precision, |
| 1672 | uint32_t options) { |
| 1673 | size_t n = (h->negative || |
| 1674 | (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN)) |
| 1675 | ? 1 |
| 1676 | : 0; |
| 1677 | if (h->decimal_point <= 0) { |
| 1678 | n += 1; |
| 1679 | } else { |
| 1680 | n += (size_t)(h->decimal_point); |
| 1681 | } |
| 1682 | if (precision > 0) { |
| 1683 | n += precision + 1; // +1 for the '.'. |
| 1684 | } |
| 1685 | |
| 1686 | // Don't modify dst if the formatted number won't fit. |
| 1687 | if (n > dst.len) { |
| 1688 | return 0; |
| 1689 | } |
| 1690 | |
| 1691 | // Align-left or align-right. |
| 1692 | uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT) |
| 1693 | ? &dst.ptr[dst.len - n] |
| 1694 | : &dst.ptr[0]; |
| 1695 | |
| 1696 | // Leading "±". |
| 1697 | if (h->negative) { |
| 1698 | *ptr++ = '-'; |
| 1699 | } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) { |
| 1700 | *ptr++ = '+'; |
| 1701 | } |
| 1702 | |
| 1703 | // Integral digits. |
| 1704 | if (h->decimal_point <= 0) { |
| 1705 | *ptr++ = '0'; |
| 1706 | } else { |
| 1707 | uint32_t m = |
| 1708 | wuffs_base__u32__min(h->num_digits, (uint32_t)(h->decimal_point)); |
| 1709 | uint32_t i = 0; |
| 1710 | for (; i < m; i++) { |
| 1711 | *ptr++ = (uint8_t)('0' | h->digits[i]); |
| 1712 | } |
| 1713 | for (; i < (uint32_t)(h->decimal_point); i++) { |
| 1714 | *ptr++ = '0'; |
| 1715 | } |
| 1716 | } |
| 1717 | |
| 1718 | // Separator and then fractional digits. |
| 1719 | if (precision > 0) { |
| 1720 | *ptr++ = |
| 1721 | (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA) |
| 1722 | ? ',' |
| 1723 | : '.'; |
| 1724 | uint32_t i = 0; |
| 1725 | for (; i < precision; i++) { |
| 1726 | uint32_t j = ((uint32_t)(h->decimal_point)) + i; |
| 1727 | *ptr++ = (uint8_t)('0' | ((j < h->num_digits) ? h->digits[j] : 0)); |
| 1728 | } |
| 1729 | } |
| 1730 | |
| 1731 | return n; |
| 1732 | } |
| 1733 | |
| 1734 | static size_t // |
| 1735 | wuffs_base__private_implementation__high_prec_dec__render_exponent_present( |
| 1736 | wuffs_base__slice_u8 dst, |
| 1737 | wuffs_base__private_implementation__high_prec_dec* h, |
| 1738 | uint32_t precision, |
| 1739 | uint32_t options) { |
| 1740 | int32_t exp = 0; |
| 1741 | if (h->num_digits > 0) { |
| 1742 | exp = h->decimal_point - 1; |
| 1743 | } |
| 1744 | bool negative_exp = exp < 0; |
| 1745 | if (negative_exp) { |
| 1746 | exp = -exp; |
| 1747 | } |
| 1748 | |
| 1749 | size_t n = (h->negative || |
| 1750 | (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN)) |
| 1751 | ? 4 |
| 1752 | : 3; // Mininum 3 bytes: first digit and then "e±". |
| 1753 | if (precision > 0) { |
| 1754 | n += precision + 1; // +1 for the '.'. |
| 1755 | } |
| 1756 | n += (exp < 100) ? 2 : 3; |
| 1757 | |
| 1758 | // Don't modify dst if the formatted number won't fit. |
| 1759 | if (n > dst.len) { |
| 1760 | return 0; |
| 1761 | } |
| 1762 | |
| 1763 | // Align-left or align-right. |
| 1764 | uint8_t* ptr = (options & WUFFS_BASE__RENDER_NUMBER_XXX__ALIGN_RIGHT) |
| 1765 | ? &dst.ptr[dst.len - n] |
| 1766 | : &dst.ptr[0]; |
| 1767 | |
| 1768 | // Leading "±". |
| 1769 | if (h->negative) { |
| 1770 | *ptr++ = '-'; |
| 1771 | } else if (options & WUFFS_BASE__RENDER_NUMBER_XXX__LEADING_PLUS_SIGN) { |
| 1772 | *ptr++ = '+'; |
| 1773 | } |
| 1774 | |
| 1775 | // Integral digit. |
| 1776 | if (h->num_digits > 0) { |
| 1777 | *ptr++ = (uint8_t)('0' | h->digits[0]); |
| 1778 | } else { |
| 1779 | *ptr++ = '0'; |
| 1780 | } |
| 1781 | |
| 1782 | // Separator and then fractional digits. |
| 1783 | if (precision > 0) { |
| 1784 | *ptr++ = |
| 1785 | (options & WUFFS_BASE__RENDER_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA) |
| 1786 | ? ',' |
| 1787 | : '.'; |
| 1788 | uint32_t i = 1; |
| 1789 | uint32_t j = wuffs_base__u32__min(h->num_digits, precision + 1); |
| 1790 | for (; i < j; i++) { |
| 1791 | *ptr++ = (uint8_t)('0' | h->digits[i]); |
| 1792 | } |
| 1793 | for (; i <= precision; i++) { |
| 1794 | *ptr++ = '0'; |
| 1795 | } |
| 1796 | } |
| 1797 | |
| 1798 | // Exponent: "e±" and then 2 or 3 digits. |
| 1799 | *ptr++ = 'e'; |
| 1800 | *ptr++ = negative_exp ? '-' : '+'; |
| 1801 | if (exp < 10) { |
| 1802 | *ptr++ = '0'; |
| 1803 | *ptr++ = (uint8_t)('0' | exp); |
| 1804 | } else if (exp < 100) { |
| 1805 | *ptr++ = (uint8_t)('0' | (exp / 10)); |
| 1806 | *ptr++ = (uint8_t)('0' | (exp % 10)); |
| 1807 | } else { |
| 1808 | int32_t e = exp / 100; |
| 1809 | exp -= e * 100; |
| 1810 | *ptr++ = (uint8_t)('0' | e); |
| 1811 | *ptr++ = (uint8_t)('0' | (exp / 10)); |
| 1812 | *ptr++ = (uint8_t)('0' | (exp % 10)); |
| 1813 | } |
| 1814 | |
| 1815 | return n; |
| 1816 | } |
| 1817 | |
| 1818 | WUFFS_BASE__MAYBE_STATIC size_t // |
| 1819 | wuffs_base__render_number_f64(wuffs_base__slice_u8 dst, |
| 1820 | double x, |
| 1821 | uint32_t precision, |
| 1822 | uint32_t options) { |
| 1823 | // Decompose x (64 bits) into negativity (1 bit), base-2 exponent (11 bits |
| 1824 | // with a -1023 bias) and mantissa (52 bits). |
Nigel Tao | 4d449dc | 2020-07-12 11:00:47 +1000 | [diff] [blame] | 1825 | uint64_t bits = wuffs_base__ieee_754_bit_representation__from_f64_to_u64(x); |
Nigel Tao | 2a7e1ed | 2020-07-07 21:50:06 +1000 | [diff] [blame] | 1826 | bool neg = (bits >> 63) != 0; |
| 1827 | int32_t exp2 = ((int32_t)(bits >> 52)) & 0x7FF; |
| 1828 | uint64_t man = bits & 0x000FFFFFFFFFFFFFul; |
| 1829 | |
| 1830 | // Apply the exponent bias and set the implicit top bit of the mantissa, |
| 1831 | // unless x is subnormal. Also take care of Inf and NaN. |
| 1832 | if (exp2 == 0x7FF) { |
| 1833 | if (man != 0) { |
| 1834 | return wuffs_base__private_implementation__render_nan(dst); |
| 1835 | } |
| 1836 | return wuffs_base__private_implementation__render_inf(dst, neg, options); |
| 1837 | } else if (exp2 == 0) { |
| 1838 | exp2 = -1022; |
| 1839 | } else { |
| 1840 | exp2 -= 1023; |
| 1841 | man |= 0x0010000000000000ul; |
| 1842 | } |
| 1843 | |
| 1844 | // Ensure that precision isn't too large. |
| 1845 | if (precision > 4095) { |
| 1846 | precision = 4095; |
| 1847 | } |
| 1848 | |
| 1849 | // Convert from the (neg, exp2, man) tuple to an HPD. |
| 1850 | wuffs_base__private_implementation__high_prec_dec h; |
| 1851 | wuffs_base__private_implementation__high_prec_dec__assign(&h, man, neg); |
| 1852 | if (h.num_digits > 0) { |
| 1853 | wuffs_base__private_implementation__high_prec_dec__lshift( |
| 1854 | &h, exp2 - 52); // 52 mantissa bits. |
| 1855 | } |
| 1856 | |
| 1857 | // Handle the "%e" and "%f" formats. |
| 1858 | switch (options & (WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT | |
| 1859 | WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT)) { |
| 1860 | case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_ABSENT: // The "%"f" format. |
| 1861 | if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) { |
| 1862 | wuffs_base__private_implementation__high_prec_dec__round_just_enough( |
| 1863 | &h, exp2, man); |
| 1864 | int32_t p = ((int32_t)(h.num_digits)) - h.decimal_point; |
| 1865 | precision = ((uint32_t)(wuffs_base__i32__max(0, p))); |
| 1866 | } else { |
| 1867 | wuffs_base__private_implementation__high_prec_dec__round_nearest( |
| 1868 | &h, ((int32_t)precision) + h.decimal_point); |
| 1869 | } |
| 1870 | return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent( |
| 1871 | dst, &h, precision, options); |
| 1872 | |
| 1873 | case WUFFS_BASE__RENDER_NUMBER_FXX__EXPONENT_PRESENT: // The "%e" format. |
| 1874 | if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) { |
| 1875 | wuffs_base__private_implementation__high_prec_dec__round_just_enough( |
| 1876 | &h, exp2, man); |
| 1877 | precision = (h.num_digits > 0) ? (h.num_digits - 1) : 0; |
| 1878 | } else { |
| 1879 | wuffs_base__private_implementation__high_prec_dec__round_nearest( |
| 1880 | &h, ((int32_t)precision) + 1); |
| 1881 | } |
| 1882 | return wuffs_base__private_implementation__high_prec_dec__render_exponent_present( |
| 1883 | dst, &h, precision, options); |
| 1884 | } |
| 1885 | |
| 1886 | // We have the "%g" format and so precision means the number of significant |
| 1887 | // digits, not the number of digits after the decimal separator. Perform |
| 1888 | // rounding and determine whether to use "%e" or "%f". |
| 1889 | int32_t e_threshold = 0; |
| 1890 | if (options & WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION) { |
| 1891 | wuffs_base__private_implementation__high_prec_dec__round_just_enough( |
| 1892 | &h, exp2, man); |
| 1893 | precision = h.num_digits; |
| 1894 | e_threshold = 6; |
| 1895 | } else { |
| 1896 | if (precision == 0) { |
| 1897 | precision = 1; |
| 1898 | } |
| 1899 | wuffs_base__private_implementation__high_prec_dec__round_nearest( |
| 1900 | &h, ((int32_t)precision)); |
| 1901 | e_threshold = ((int32_t)precision); |
| 1902 | int32_t nd = ((int32_t)(h.num_digits)); |
| 1903 | if ((e_threshold > nd) && (nd >= h.decimal_point)) { |
| 1904 | e_threshold = nd; |
| 1905 | } |
| 1906 | } |
| 1907 | |
| 1908 | // Use the "%e" format if the exponent is large. |
| 1909 | int32_t e = h.decimal_point - 1; |
| 1910 | if ((e < -4) || (e_threshold <= e)) { |
| 1911 | uint32_t p = wuffs_base__u32__min(precision, h.num_digits); |
| 1912 | return wuffs_base__private_implementation__high_prec_dec__render_exponent_present( |
| 1913 | dst, &h, (p > 0) ? (p - 1) : 0, options); |
| 1914 | } |
| 1915 | |
| 1916 | // Use the "%f" format otherwise. |
| 1917 | int32_t p = ((int32_t)precision); |
| 1918 | if (p > h.decimal_point) { |
| 1919 | p = ((int32_t)(h.num_digits)); |
| 1920 | } |
| 1921 | precision = ((uint32_t)(wuffs_base__i32__max(0, p - h.decimal_point))); |
| 1922 | return wuffs_base__private_implementation__high_prec_dec__render_exponent_absent( |
| 1923 | dst, &h, precision, options); |
| 1924 | } |