// This module was based on Rust's dec2flt // https://github.com/rust-lang/rust/blob/1cbc45942d5c0f6eb5d94e3b10762ba541958035/library/core/src/num/dec2flt/lemire.rs // Rust's MIT license is provided below: /* * Permission is hereby granted, free of charge, to any * person obtaining a copy of this software and associated * documentation files (the "Software"), to deal in the * Software without restriction, including without * limitation the rights to use, copy, modify, merge, * publish, distribute, sublicense, and/or sell copies of * the Software, and to permit persons to whom the Software * is furnished to do so, subject to the following * conditions: * * The above copyright notice and this permission notice * shall be included in all copies or substantial portions * of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF * ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED * TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT * SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR * IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. */ @noPervasives module Lemire from "runtime/unsafe/wasmi32" include WasmI32 from "runtime/unsafe/wasmi64" include WasmI64 from "runtime/dataStructures" include DataStructures use DataStructures.{ newInt32, newInt64, newFloat64 } from "runtime/atof/table" include Table use Table.{ get_POWERS5 } from "runtime/atof/common" include Common use Common.{ _SMALLEST_POWER_OF_FIVE, _LARGEST_POWER_OF_FIVE, _INFINITE_POWER, _SMALLEST_POWER_OF_TEN, _LARGEST_POWER_OF_TEN, _MANTISSA_EXPLICIT_BITS_32, _MANTISSA_EXPLICIT_BITS_64, _MINIMUM_EXPONENT, _MIN_EXPONENT_ROUND_TO_EVEN, _MAX_EXPONENT_ROUND_TO_EVEN, fullMultiplication, power, type BiasedFp, fpZero, fpInf, fpErr, } primitive (&&) = "@and" primitive (||) = "@or" primitive (!) = "@not" // From Rust: // This will compute or rather approximate w * 5**q and return a pair of 64-bit words // approximating the result, with the "high" part corresponding to the most significant // bits and the low part corresponding to the least significant bits. @unsafe let computeProductApprox = (q: WasmI64, w: WasmI64, precision: WasmI64) => { use WasmI64.{ (==), (&), gtU as (>), ltU as (<), (>=), (<=), (>>>), (-), (+) } use WasmI32.{ (*) } let mask = if (precision < 64N) { 0xFFFF_FFFF_FFFF_FFFFN >>> precision } else { 0xFFFF_FFFF_FFFF_FFFFN } // From Rust: // 5^q < 2^64, then the multiplication always provides an exact value. // That means whenever we need to round ties to even, we always have // an exact value. let index = q - _SMALLEST_POWER_OF_FIVE let n = 16n * WasmI32.wrapI64(index) use WasmI32.{ (+) } let lo5 = WasmI64.load(get_POWERS5(), n) and hi5 = WasmI64.load(get_POWERS5(), n + 8n) use WasmI64.{ (+) } // From Rust: // Only need one multiplication as long as there is 1 zero but // in the explicit mantissa bits, +1 for the hidden bit, +1 to // determine the rounding direction, +1 for if the computed // product has a leading zero. let (firstLo, firstHi) = fullMultiplication(w, lo5) let mut firstLo = WasmI64.load(WasmI32.fromGrain(firstLo), 8n) let mut firstHi = WasmI64.load(WasmI32.fromGrain(firstHi), 8n) if ((firstHi & mask) == mask) { // From Rust: // Need to do a second multiplication to get better precision // for the lower product. This will always be exact // where q is < 55, since 5^55 < 2^128. If this wraps, // then we need to need to round up the hi product. let (_, secondHi) = fullMultiplication(w, hi5) let secondHi = WasmI64.load(WasmI32.fromGrain(secondHi), 8n) firstLo += secondHi if (secondHi > firstLo) { firstHi += 1N } } ( WasmI32.toGrain(newInt64(firstLo)): Int64, WasmI32.toGrain(newInt64(firstHi)): Int64, ) } // Returns the significant bits and a biased binary exponent that can be directly shifted into exponent bits: // (u64, i32) @unsafe provide let computeFloat = (exponent: WasmI64, mantissa: WasmI64) => { use WasmI64.{ (==), (<), (<=), (>), (>=), (<<), (>>), (&), (+), (>>>) } let w = mantissa let q = exponent if (w == 0N || q < _SMALLEST_POWER_OF_TEN) { fpZero() } else if (q > _LARGEST_POWER_OF_TEN) { fpInf() } else { let lz = WasmI64.clz(w) let w = w << lz let (lo, hi) = computeProductApprox(q, w, _MANTISSA_EXPLICIT_BITS_64 + 3N) let lo = WasmI64.load(WasmI32.fromGrain(lo), 8n) let hi = WasmI64.load(WasmI32.fromGrain(hi), 8n) // From Rust: // If we have failed to approximate w x 5^-q with our 128-bit value. // Since the addition of 1 could lead to an overflow which could then // round up over the half-way point, this can lead to improper rounding // of a float. // // However, this can only occur if q ∈ [-27, 55]. The upper bound of q // is 55 because 5^55 < 2^128, however, this can only happen if 5^q > 2^64, // since otherwise the product can be represented in 64-bits, producing // an exact result. For negative exponents, rounding-to-even can // only occur if 5^-q < 2^64. // // For detailed explanations of rounding for negative exponents, see // . For detailed // explanations of rounding for positive exponents, see // . match (lo == 0xFFFF_FFFF_FFFF_FFFFN && !(q >= -27N && q <= 55N)) { true => fpErr(), false => { let upperbit = WasmI32.wrapI64(hi >>> 63N) let mut mantissa = hi >>> WasmI64.extendI32S({ use WasmI32.{ (+), (-) } upperbit + 64n - _MANTISSA_EXPLICIT_BITS_32 - 3n }) let mut power2 = { use WasmI32.{ (+), (-) } let q = WasmI32.wrapI64(q) let lz = WasmI32.wrapI64(lz) power(q) + upperbit - lz - _MINIMUM_EXPONENT } use WasmI32.{ (<=) } if (power2 <= 0n) { use WasmI32.{ (+), (*), (>=) } // -power2 + 1 >= 64 if (-1n * power2 + 1n >= 64n) { // Have more than 64 bits below the minimum exponent, must be 0. fpZero() } else { use WasmI64.{ (>=) } // Have a subnormal value. mantissa = mantissa >> WasmI64.extendI32S(-1n * power2 + 1n) use WasmI64.{ (+) } mantissa += mantissa & 1N mantissa = mantissa >> 1N power2 = if (mantissa >= 1N << _MANTISSA_EXPLICIT_BITS_64) { 1n } else { 0n } { f: WasmI32.toGrain(newInt64(mantissa)): Int64, e: WasmI32.toGrain(newInt32(power2)): Int32, } } } else { use WasmI64.{ (<=) } // Need to handle rounding ties. Normally, we need to round up, // but if we fall right in between and and we have an even basis, we // need to round down. // // This will only occur if: // 1. The lower 64 bits of the 128-bit representation is 0. // IE, 5^q fits in single 64-bit word. // 2. The least-significant bit prior to truncated mantissa is odd. // 3. All the bits truncated when shifting to mantissa bits + 1 are 0. // // Or, we may fall between two floats: we are exactly halfway. if ( WasmI64.leU(lo, 1N) && q >= _MIN_EXPONENT_ROUND_TO_EVEN && q <= _MAX_EXPONENT_ROUND_TO_EVEN && (mantissa & 3N) == 1N && mantissa << WasmI64.extendI32S({ use WasmI32.{ (+), (-) } upperbit + 64n - _MANTISSA_EXPLICIT_BITS_32 - 3n }) == hi ) { use WasmI64.{ (^) } // Zero the lowest bit, so we don't round up. // mantissa &= !1N; mantissa = mantissa & (1N ^ -1N) } // Round-to-even, then shift the significant digits into place. mantissa += mantissa & 1N use WasmI64.{ (>>>) } mantissa = mantissa >>> 1N if (WasmI64.geU(mantissa, 2N << _MANTISSA_EXPLICIT_BITS_64)) { use WasmI32.{ (+) } // Rounding up overflowed, so the carry bit is set. Set the // mantissa to 1 (only the implicit, hidden bit is set) and // increase the exponent. mantissa = 1N << _MANTISSA_EXPLICIT_BITS_64 power2 += 1n } use WasmI32.{ (>=) } use WasmI64.{ (^) } // Zero out the hidden bit. mantissa = mantissa & (1N << _MANTISSA_EXPLICIT_BITS_64 ^ -1N) if (power2 >= _INFINITE_POWER) { // Exponent is above largest normal value, must be infinite. fpInf() } else { { f: WasmI32.toGrain(newInt64(mantissa)): Int64, e: WasmI32.toGrain(newInt32(power2)): Int32, } } } }, } } }