436 lines
11 KiB
V
436 lines
11 KiB
V
/**********************************************************************
|
||
*
|
||
* f32 to string
|
||
*
|
||
* Copyright (c) 2019-2020 Dario Deledda. All rights reserved.
|
||
* Use of this source code is governed by an MIT license
|
||
* that can be found in the LICENSE file.
|
||
*
|
||
* This file contains the f64 to string functions
|
||
*
|
||
* These functions are based on the work of:
|
||
* Publication:PLDI 2018: Proceedings of the 39th ACM SIGPLAN
|
||
* Conference on Programming Language Design and ImplementationJune 2018
|
||
* Pages 270–282 https://doi.org/10.1145/3192366.3192369
|
||
*
|
||
* inspired by the Go version here:
|
||
* https://github.com/cespare/ryu/tree/ba56a33f39e3bbbfa409095d0f9ae168a595feea
|
||
*
|
||
**********************************************************************/
|
||
module ftoa
|
||
|
||
struct Uint128 {
|
||
mut:
|
||
lo u64 = u64(0)
|
||
hi u64 = u64(0)
|
||
}
|
||
|
||
// dec64 is a floating decimal type representing m * 10^e.
|
||
struct Dec64 {
|
||
mut:
|
||
m u64 = 0
|
||
e int = 0
|
||
}
|
||
|
||
// support union for convert f64 to u64
|
||
union Uf64 {
|
||
mut:
|
||
f f64 = 0
|
||
u u64
|
||
}
|
||
|
||
// pow of ten table used by n_digit reduction
|
||
const(
|
||
ten_pow_table_64 = [
|
||
u64(1),
|
||
u64(10),
|
||
u64(100),
|
||
u64(1000),
|
||
u64(10000),
|
||
u64(100000),
|
||
u64(1000000),
|
||
u64(10000000),
|
||
u64(100000000),
|
||
u64(1000000000),
|
||
u64(10000000000),
|
||
u64(100000000000),
|
||
u64(1000000000000),
|
||
u64(10000000000000),
|
||
u64(100000000000000),
|
||
u64(1000000000000000),
|
||
u64(10000000000000000),
|
||
u64(100000000000000000),
|
||
u64(1000000000000000000),
|
||
u64(10000000000000000000),
|
||
]
|
||
)
|
||
|
||
/******************************************************************************
|
||
*
|
||
* Conversion Functions
|
||
*
|
||
******************************************************************************/
|
||
const(
|
||
mantbits64 = u32(52)
|
||
expbits64 = u32(11)
|
||
bias64 = u32(1023) // f64 exponent bias
|
||
maxexp64 = 2047
|
||
)
|
||
|
||
fn (d Dec64) get_string_64(neg bool, i_n_digit int, i_pad_digit int) string {
|
||
mut n_digit := i_n_digit + 1
|
||
pad_digit := i_pad_digit + 1
|
||
mut out := d.m
|
||
mut d_exp := d.e
|
||
mut out_len := decimal_len_64(out)
|
||
out_len_original := out_len
|
||
|
||
mut fw_zeros := 0
|
||
if pad_digit > out_len {
|
||
fw_zeros = pad_digit - out_len
|
||
}
|
||
|
||
mut buf := [byte(0)].repeat(out_len + 6 + 1 +1 + fw_zeros) // sign + mant_len + . + e + e_sign + exp_len(2) + \0
|
||
mut i := 0
|
||
|
||
if neg {
|
||
buf[i]=`-`
|
||
i++
|
||
}
|
||
|
||
mut disp := 0
|
||
if out_len <= 1 {
|
||
disp = 1
|
||
}
|
||
|
||
// rounding last used digit
|
||
if n_digit < out_len {
|
||
//println("out:[$out]")
|
||
out += ten_pow_table_64[out_len - n_digit - 1] * 5 // round to up
|
||
out /= ten_pow_table_64[out_len - n_digit ]
|
||
//println("out1:[$out] ${d.m / ten_pow_table_64[out_len - n_digit ]}")
|
||
if d.m / ten_pow_table_64[out_len - n_digit ] < out {
|
||
d_exp += 1
|
||
n_digit += 1
|
||
}
|
||
|
||
//println("cmp: ${d.m/ten_pow_table_64[out_len - n_digit ]} ${out/ten_pow_table_64[out_len - n_digit ]}")
|
||
|
||
out_len = n_digit
|
||
//println("orig: ${out_len_original} new len: ${out_len} out:[$out]")
|
||
}
|
||
|
||
y := i + out_len
|
||
mut x := 0
|
||
for x < (out_len-disp-1) {
|
||
buf[y - x] = `0` + byte(out%10)
|
||
out /= 10
|
||
i++
|
||
x++
|
||
}
|
||
|
||
if out_len >= 1 {
|
||
buf[y - x] = `.`
|
||
x++
|
||
i++
|
||
}
|
||
|
||
if y-x >= 0 {
|
||
buf[y - x] = `0` + byte(out%10)
|
||
i++
|
||
}
|
||
|
||
for fw_zeros > 0 {
|
||
buf[i++] = `0`
|
||
fw_zeros--
|
||
}
|
||
|
||
/*
|
||
x=0
|
||
for x<buf.len {
|
||
C.printf("d:%c\n",buf[x])
|
||
x++
|
||
}
|
||
C.printf("\n")
|
||
*/
|
||
|
||
buf[i]=`e`
|
||
i++
|
||
|
||
mut exp := d_exp + out_len_original - 1
|
||
if exp < 0 {
|
||
buf[i]=`-`
|
||
i++
|
||
exp = -exp
|
||
} else {
|
||
buf[i]=`+`
|
||
i++
|
||
}
|
||
|
||
// Always print at least two digits to match strconv's formatting.
|
||
d2 := exp % 10
|
||
exp /= 10
|
||
d1 := exp % 10
|
||
d0 := exp / 10
|
||
if d0 > 0 {
|
||
buf[i]=`0` + byte(d0)
|
||
i++
|
||
}
|
||
buf[i]=`0` + byte(d1)
|
||
i++
|
||
buf[i]=`0` + byte(d2)
|
||
i++
|
||
buf[i]=0
|
||
|
||
|
||
/*
|
||
x=0
|
||
for x<buf.len {
|
||
C.printf("d:%c\n",buf[x])
|
||
x++
|
||
}
|
||
*/
|
||
return tos(byteptr(&buf[0]), i)
|
||
}
|
||
|
||
fn f64_to_decimal_exact_int(i_mant u64, exp u64) (Dec64, bool) {
|
||
mut d := Dec64{}
|
||
e := exp - bias64
|
||
if e > mantbits64 {
|
||
return d, false
|
||
}
|
||
shift := mantbits64 - e
|
||
mant := i_mant | u64(0x0010_0000_0000_0000) // implicit 1
|
||
//mant := i_mant | (1 << mantbits64) // implicit 1
|
||
d.m = mant >> shift
|
||
if (d.m << shift) != mant {
|
||
return d, false
|
||
}
|
||
|
||
for (d.m % 10) == 0 {
|
||
d.m /= 10
|
||
d.e++
|
||
}
|
||
return d, true
|
||
}
|
||
|
||
fn f64_to_decimal(mant u64, exp u64) Dec64 {
|
||
mut e2 := 0
|
||
mut m2 := u64(0)
|
||
if exp == 0 {
|
||
// We subtract 2 so that the bounds computation has
|
||
// 2 additional bits.
|
||
e2 = 1 - bias64 - mantbits64 - 2
|
||
m2 = mant
|
||
} else {
|
||
e2 = int(exp) - bias64 - mantbits64 - 2
|
||
m2 = (u64(1)<<mantbits64) | mant
|
||
}
|
||
even := (m2 & 1) == 0
|
||
accept_bounds := even
|
||
|
||
// Step 2: Determine the interval of valid decimal representations.
|
||
mv := u64(4 * m2)
|
||
mm_shift := bool_to_u64(mant != 0 || exp <= 1)
|
||
|
||
// Step 3: Convert to a decimal power base uing 128-bit arithmetic.
|
||
mut vr := u64(0)
|
||
mut vp := u64(0)
|
||
mut vm := u64(0)
|
||
mut e10 := 0
|
||
mut vm_is_trailing_zeros := false
|
||
mut vr_is_trailing_zeros := false
|
||
|
||
if e2 >= 0 {
|
||
// This expression is slightly faster than max(0, log10Pow2(e2) - 1).
|
||
q := log10_pow2(e2) - bool_to_u32(e2 > 3)
|
||
e10 = int(q)
|
||
k := pow5_inv_num_bits_64 + pow5_bits(int(q)) - 1
|
||
i := -e2 + int(q) + k
|
||
|
||
mul := pow5_inv_split_64[q]
|
||
vr = mul_shift_64(u64(4) * m2 , mul, i)
|
||
vp = mul_shift_64(u64(4) * m2 + u64(2) , mul, i)
|
||
vm = mul_shift_64(u64(4) * m2 - u64(1) - mm_shift, mul, i)
|
||
if q <= 21 {
|
||
// This should use q <= 22, but I think 21 is also safe.
|
||
// Smaller values may still be safe, but it's more
|
||
// difficult to reason about them. Only one of mp, mv,
|
||
// and mm can be a multiple of 5, if any.
|
||
if mv%5 == 0 {
|
||
vr_is_trailing_zeros = multiple_of_power_of_five_64(mv, q)
|
||
} else if accept_bounds {
|
||
// Same as min(e2 + (^mm & 1), pow5Factor64(mm)) >= q
|
||
// <=> e2 + (^mm & 1) >= q && pow5Factor64(mm) >= q
|
||
// <=> true && pow5Factor64(mm) >= q, since e2 >= q.
|
||
vm_is_trailing_zeros = multiple_of_power_of_five_64(mv-1-mm_shift, q)
|
||
} else if multiple_of_power_of_five_64(mv+2, q) {
|
||
vp--
|
||
}
|
||
}
|
||
} else {
|
||
// This expression is slightly faster than max(0, log10Pow5(-e2) - 1).
|
||
q := log10_pow5(-e2) - bool_to_u32(-e2 > 1)
|
||
e10 = int(q) + e2
|
||
i := -e2 - int(q)
|
||
k := pow5_bits(i) - pow5_num_bits_64
|
||
mut j := int(q) - k
|
||
mul := pow5_split_64[i]
|
||
vr = mul_shift_64(u64(4) * m2 , mul, j)
|
||
vp = mul_shift_64(u64(4) * m2 + u64(2) , mul, j)
|
||
vm = mul_shift_64(u64(4) * m2 - u64(1) - mm_shift, mul, j)
|
||
if q <= 1 {
|
||
// {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits.
|
||
// mv = 4 * m2, so it always has at least two trailing 0 bits.
|
||
vr_is_trailing_zeros = true
|
||
if accept_bounds {
|
||
// mm = mv - 1 - mmShift, so it has 1 trailing 0 bit iff mmShift == 1.
|
||
vm_is_trailing_zeros = (mm_shift == 1)
|
||
} else {
|
||
// mp = mv + 2, so it always has at least one trailing 0 bit.
|
||
vp--
|
||
}
|
||
} else if q < 63 { // TODO(ulfjack/cespare): Use a tighter bound here.
|
||
// We need to compute min(ntz(mv), pow5Factor64(mv) - e2) >= q - 1
|
||
// <=> ntz(mv) >= q - 1 && pow5Factor64(mv) - e2 >= q - 1
|
||
// <=> ntz(mv) >= q - 1 (e2 is negative and -e2 >= q)
|
||
// <=> (mv & ((1 << (q - 1)) - 1)) == 0
|
||
// We also need to make sure that the left shift does not overflow.
|
||
vr_is_trailing_zeros = multiple_of_power_of_two_64(mv, q - 1)
|
||
}
|
||
}
|
||
|
||
// Step 4: Find the shortest decimal representation
|
||
// in the interval of valid representations.
|
||
mut removed := 0
|
||
mut last_removed_digit := byte(0)
|
||
mut out := u64(0)
|
||
// On average, we remove ~2 digits.
|
||
if vm_is_trailing_zeros || vr_is_trailing_zeros {
|
||
// General case, which happens rarely (~0.7%).
|
||
for {
|
||
vp_div_10 := vp / 10
|
||
vm_div_10 := vm / 10
|
||
if vp_div_10 <= vm_div_10 {
|
||
break
|
||
}
|
||
vm_mod_10 := vm % 10
|
||
vr_div_10 := vr / 10
|
||
vr_mod_10 := vr % 10
|
||
vm_is_trailing_zeros = vm_is_trailing_zeros && vm_mod_10 == 0
|
||
vr_is_trailing_zeros = vr_is_trailing_zeros && (last_removed_digit == 0)
|
||
last_removed_digit = byte(vr_mod_10)
|
||
vr = vr_div_10
|
||
vp = vp_div_10
|
||
vm = vm_div_10
|
||
removed++
|
||
}
|
||
if vm_is_trailing_zeros {
|
||
for {
|
||
vm_div_10 := vm / 10
|
||
vm_mod_10 := vm % 10
|
||
if vm_mod_10 != 0 {
|
||
break
|
||
}
|
||
vp_div_10 := vp / 10
|
||
vr_div_10 := vr / 10
|
||
vr_mod_10 := vr % 10
|
||
vr_is_trailing_zeros = vr_is_trailing_zeros && (last_removed_digit == 0)
|
||
last_removed_digit = byte(vr_mod_10)
|
||
vr = vr_div_10
|
||
vp = vp_div_10
|
||
vm = vm_div_10
|
||
removed++
|
||
}
|
||
}
|
||
if vr_is_trailing_zeros && (last_removed_digit == 5) && (vr % 2) == 0 {
|
||
// Round even if the exact number is .....50..0.
|
||
last_removed_digit = 4
|
||
}
|
||
out = vr
|
||
// We need to take vr + 1 if vr is outside bounds
|
||
// or we need to round up.
|
||
if (vr == vm && (!accept_bounds || !vm_is_trailing_zeros)) || last_removed_digit >= 5 {
|
||
out++
|
||
}
|
||
} else {
|
||
// Specialized for the common case (~99.3%).
|
||
// Percentages below are relative to this.
|
||
mut round_up := false
|
||
for vp / 100 > vm / 100 {
|
||
// Optimization: remove two digits at a time (~86.2%).
|
||
round_up = (vr % 100) >= 50
|
||
vr /= 100
|
||
vp /= 100
|
||
vm /= 100
|
||
removed += 2
|
||
}
|
||
// Loop iterations below (approximately), without optimization above:
|
||
// 0: 0.03%, 1: 13.8%, 2: 70.6%, 3: 14.0%, 4: 1.40%, 5: 0.14%, 6+: 0.02%
|
||
// Loop iterations below (approximately), with optimization above:
|
||
// 0: 70.6%, 1: 27.8%, 2: 1.40%, 3: 0.14%, 4+: 0.02%
|
||
for vp / 10 > vm / 10 {
|
||
round_up = (vr % 10) >= 5
|
||
vr /= 10
|
||
vp /= 10
|
||
vm /= 10
|
||
removed++
|
||
}
|
||
// We need to take vr + 1 if vr is outside bounds
|
||
// or we need to round up.
|
||
out = vr + bool_to_u64(vr == vm || round_up)
|
||
}
|
||
|
||
return Dec64{m: out, e: e10 + removed}
|
||
}
|
||
|
||
// f64_to_str return a string in scientific notation with max n_digit after the dot
|
||
pub fn f64_to_str(f f64, n_digit int) string {
|
||
mut u1 := Uf64{}
|
||
u1.f = f
|
||
u := u1.u
|
||
|
||
neg := (u>>(mantbits64+expbits64)) != 0
|
||
mant := u & ((u64(1)<<mantbits64) - u64(1))
|
||
exp := (u >> mantbits64) & ((u64(1)<<expbits64) - u64(1))
|
||
//println("s:${neg} mant:${mant} exp:${exp} float:${f} byte:${u1.u:016lx}")
|
||
|
||
// Exit early for easy cases.
|
||
if (exp == maxexp64) || (exp == 0 && mant == 0) {
|
||
return get_string_special(neg, exp == 0, mant == 0)
|
||
}
|
||
|
||
mut d, ok := f64_to_decimal_exact_int(mant, exp)
|
||
if !ok {
|
||
//println("to_decimal")
|
||
d = f64_to_decimal(mant, exp)
|
||
}
|
||
//println("${d.m} ${d.e}")
|
||
return d.get_string_64(neg, n_digit, 0)
|
||
}
|
||
|
||
// f64_to_str return a string in scientific notation with max n_digit after the dot
|
||
pub fn f64_to_str_pad(f f64, n_digit int) string {
|
||
mut u1 := Uf64{}
|
||
u1.f = f
|
||
u := u1.u
|
||
|
||
neg := (u>>(mantbits64+expbits64)) != 0
|
||
mant := u & ((u64(1)<<mantbits64) - u64(1))
|
||
exp := (u >> mantbits64) & ((u64(1)<<expbits64) - u64(1))
|
||
//println("s:${neg} mant:${mant} exp:${exp} float:${f} byte:${u1.u:016lx}")
|
||
|
||
// Exit early for easy cases.
|
||
if (exp == maxexp64) || (exp == 0 && mant == 0) {
|
||
return get_string_special(neg, exp == 0, mant == 0)
|
||
}
|
||
|
||
mut d, ok := f64_to_decimal_exact_int(mant, exp)
|
||
if !ok {
|
||
//println("to_decimal")
|
||
d = f64_to_decimal(mant, exp)
|
||
}
|
||
//println("DEBUG: ${d.m} ${d.e}")
|
||
return d.get_string_64(neg, n_digit, n_digit)
|
||
}
|