lemire.gr

// 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
    // <https://arxiv.org/pdf/2101.11408.pdf#section.9.1>. For detailed
    // explanations of rounding for positive exponents, see
    // <https://arxiv.org/pdf/2101.11408.pdf#section.8>.
    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,
            }
          }
        }
      },
    }
  }
}