From 200f25a38f4fe6db2cdbb96511a2d504f67cba89 Mon Sep 17 00:00:00 2001 From: penguindark <57967770+penguindark@users.noreply.github.com> Date: Tue, 25 Feb 2020 11:12:37 +0100 Subject: [PATCH] ftoa in V (#3831) --- vlib/strconv/ftoa/f32_f64_to_string_test.v | 168 +++++ vlib/strconv/ftoa/f32_str.v | 322 ++++++++++ vlib/strconv/ftoa/f64_str.v | 367 +++++++++++ vlib/strconv/ftoa/ftoa.v | 40 ++ vlib/strconv/ftoa/tables.v | 680 +++++++++++++++++++++ vlib/strconv/ftoa/utilities.v | 336 ++++++++++ 6 files changed, 1913 insertions(+) create mode 100644 vlib/strconv/ftoa/f32_f64_to_string_test.v create mode 100644 vlib/strconv/ftoa/f32_str.v create mode 100644 vlib/strconv/ftoa/f64_str.v create mode 100644 vlib/strconv/ftoa/ftoa.v create mode 100644 vlib/strconv/ftoa/tables.v create mode 100644 vlib/strconv/ftoa/utilities.v diff --git a/vlib/strconv/ftoa/f32_f64_to_string_test.v b/vlib/strconv/ftoa/f32_f64_to_string_test.v new file mode 100644 index 0000000000..85a8bbe37c --- /dev/null +++ b/vlib/strconv/ftoa/f32_f64_to_string_test.v @@ -0,0 +1,168 @@ +/********************************************************************** +* +* Float to string Test +* +**********************************************************************/ +module ftoa +import math + +union Ufloat32 { +mut: + f f32 = f32(0) + b u32 +} + +union Ufloat64 { +mut: + f f64 = f64(0) + b u64 +} + +fn f64_from_bits1(b u64) f64 { + mut x := Ufloat64{} + x.b = b + //C.printf("bin: %016llx\n",x.f) + return x.f +} + +fn f32_from_bits1(b u32) f32 { + mut x := Ufloat32{} + x.b = b + //C.printf("bin: %08x\n",x.f) + return x.f +} + +const( +test_cases_f32 = [ + f32_from_bits1(0x0000_0000), // +0 + f32_from_bits1(0x8000_0000), // -0 + f32_from_bits1(0xFFC0_0001), // sNan + f32_from_bits1(0xFF80_0001), // qNan + f32_from_bits1(0x7F80_0000), // +inf + f32_from_bits1(0xFF80_0000), // -inf + 1, + -1, + 10, + -10, + 0.3, + -0.3, + 1000000, + 123456.7, + 123e35, + -123.45, + 1e23, + f32_from_bits1(0x0080_0000), // smallest float32 + math.max_f32, + 383260575764816448, +] + +exp_result_f32 = [ + "0e+00", + "-0e+00", + "NaN", + "NaN", + "+inf", + "-inf", + "1.e+00", + "-1.e+00", + "1.e+01", + "-1.e+01", + "3.e-01", + "-3.e-01", + "1.e+06", + "1.234567e+05", + "1.23e+37", + "-1.2345e+02", + "1.e+23", + "1.1754944e-38", // aprox from 1.1754943508 × 10−38, + "3.4028235e+38", + "3.8326058e+17", +] + +test_cases_f64 = [ + f64_from_bits1(0x0000_0000_0000_0000), // +0 + f64_from_bits1(0x8000_0000_0000_0000), // -0 + f64_from_bits1(0x7FF0_0000_0000_0001), // sNan + f64_from_bits1(0x7FF8_0000_0000_0001), // qNan + f64_from_bits1(0x7FF0_0000_0000_0000), // +inf + f64_from_bits1(0xFFF0_0000_0000_0000), // -inf + 1, + -1, + 10, + -10, + 0.3, + -0.3, + 1000000, + 123456.7, + 123e45, + -123.45, + 1e23, + f64_from_bits1(0x0010_0000_0000_0000), // smallest float64 + math.max_f32, + 383260575764816448, + 383260575764816448, + // C failing cases + 123e300, + 123e-300, + 5.e-324, + -5.e-324, +] + +exp_result_f64 = [ + "0e+00", + "-0e+00", + "NaN", + "NaN", + "+inf", + "-inf", + "1.e+00", + "-1.e+00", + "1.e+01", + "-1.e+01", + "3.e-01", + "-3.e-01", + "1.e+06", + "1.234567e+05", + "1.23e+47", + "-1.2345e+02", + "1.e+23", + "2.2250738585072014e-308", + "3.4028234663852886e+38", + "3.8326057576481645e+17", + "3.8326057576481645e+17", + + "1.23e+302", // this test is failed from C sprintf!! + "1.23e-298", + "5.e-324", + "-5.e-324", +] +) + +fn test_float_to_str(){ + // test f32 + for c,x in test_cases_f32 { + s := ftoa.f32_to_str(x,8) + s1 := exp_result_f32[c] + //println("$s1 $s") + assert s == s1 + } + + // test f64 + for c,x in test_cases_f64 { + s := ftoa.f64_to_str(x,17) + s1 := exp_result_f64[c] + //println("$s1 $s") + assert s == s1 + } + + // test long format + for exp := 1 ; exp < 120 ; exp++ { + a :=ftoa.f64_to_str_l(("1e"+exp.str()).f64()) + //println(a) + assert a.len == exp + 1 + + b :=ftoa.f64_to_str_l(("1e-"+exp.str()).f64()) + //println(b) + assert b.len == exp + 2 + } +} diff --git a/vlib/strconv/ftoa/f32_str.v b/vlib/strconv/ftoa/f32_str.v new file mode 100644 index 0000000000..37d09f37d4 --- /dev/null +++ b/vlib/strconv/ftoa/f32_str.v @@ -0,0 +1,322 @@ +/********************************************************************** +* +* 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 f32 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 + +// dec32 is a floating decimal type representing m * 10^e. +struct Dec32 { +mut: + m u32 = u32(0) + e int = 0 +} + +// support union for convert f32 to u32 +union Uf32 { +mut: + f f32 = f32(0) + u u32 +} + +/****************************************************************************** +* +* Conversion Functions +* +******************************************************************************/ +const( + mantbits32 = u32(23) + expbits32 = u32(8) + bias32 = u32(127) // f32 exponent bias + maxexp32 = 255 +) + +// max 46 char +// -3.40282346638528859811704183484516925440e+38 +fn (d Dec32) get_string_32(neg bool, n_digit int) string { + mut out := d.m + mut out_len := decimal_len_32(out) + + mut buf := [byte(0)].repeat(out_len + 5 + 1 +1) // sign + mant_len + . + e + e_sign + exp_len(2) + \0 + mut i := 0 + + if n_digit > 0 && out_len > n_digit { + out_len = n_digit+1 + } + + if neg { + buf[i]=`-` + i++ + } + + mut disp := 0 + if out_len <= 1 { + disp = 1 + } + + 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++ + } + + /* + x=0 + for x mantbits32 { + return d, false + } + shift := mantbits32 - e + mant := i_mant | 0x0080_0000 // implicit 1 + //mant := i_mant | (1 << mantbits32) // 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 +} + +pub fn f32_to_decimal(mant u32, exp u32) Dec32 { + mut e2 := 0 + mut m2 := u32(0) + if exp == 0 { + // We subtract 2 so that the bounds computation has + // 2 additional bits. + e2 = 1 - bias32 - mantbits32 - 2 + m2 = mant + } else { + e2 = int(exp) - bias32 - mantbits32 - 2 + m2 = (u32(1) << mantbits32) | mant + } + even := (m2 & 1) == 0 + accept_bounds := even + + // Step 2: Determine the interval of valid decimal representations. + mv := u32(4 * m2) + mp := u32(4 * m2 + 2) + mm_shift := bool_to_u32(mant != 0 || exp <= 1) + mm := u32(4 * m2 - 1 - mm_shift) + + mut vr := u32(0) + mut vp := u32(0) + mut vm := u32(0) + mut e10 := 0 + mut vm_is_trailing_zeros := false + mut vr_is_trailing_zeros := false + mut last_removed_digit := byte(0) + + if e2 >= 0 { + q := log10_pow2(e2) + e10 = int(q) + k := pow5_inv_num_bits_32 + pow5_bits(int(q)) - 1 + i := -e2 + int(q) + k + + vr = mul_pow5_invdiv_pow2(mv, q, i) + vp = mul_pow5_invdiv_pow2(mp, q, i) + vm = mul_pow5_invdiv_pow2(mm, q, i) + if q != 0 && (vp-1)/10 <= vm/10 { + // We need to know one removed digit even if we are not + // going to loop below. We could use q = X - 1 above, + // except that would require 33 bits for the result, and + // we've found that 32-bit arithmetic is faster even on + // 64-bit machines. + l := pow5_inv_num_bits_32 + pow5_bits(int(q - 1)) - 1 + last_removed_digit = byte(mul_pow5_invdiv_pow2(mv, q - 1, -e2 + int(q - 1) + l) % 10) + } + if q <= 9 { + // The largest power of 5 that fits in 24 bits is 5^10, + // but q <= 9 seems to be safe as well. 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_32(mv, q) + } else if accept_bounds { + vm_is_trailing_zeros = multiple_of_power_of_five_32(mm, q) + } else if multiple_of_power_of_five_32(mp, q) { + vp-- + } + } + } else { + q := log10_pow5(-e2) + e10 = int(q) + e2 + i := -e2 - int(q) + k := pow5_bits(i) - pow5_num_bits_32 + mut j := int(q) - k + vr = mul_pow5_div_pow2(mv, u32(i), j) + vp = mul_pow5_div_pow2(mp, u32(i), j) + vm = mul_pow5_div_pow2(mm, u32(i), j) + if q != 0 && ((vp-1)/10) <= vm/10 { + j = int(q) - 1 - (pow5_bits(i + 1) - pow5_num_bits_32) + last_removed_digit = byte(mul_pow5_div_pow2(mv, u32(i + 1), j) % 10) + } + 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 - mm_shift, so it has 1 trailing 0 bit + // if mm_shift == 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 < 31 { + vr_is_trailing_zeros = multiple_of_power_of_two_32(mv, q - 1) + } + } + + // Step 4: Find the shortest decimal representation + // in the interval of valid representations. + mut removed := 0 + mut out := u32(0) + if vm_is_trailing_zeros || vr_is_trailing_zeros { + // General case, which happens rarely (~4.0%). + for vp/10 > vm/10 { + vm_is_trailing_zeros = vm_is_trailing_zeros && (vm % 10) == 0 + vr_is_trailing_zeros = vr_is_trailing_zeros && (last_removed_digit == 0) + last_removed_digit = byte(vr % 10) + vr /= 10 + vp /= 10 + vm /= 10 + removed++ + } + if vm_is_trailing_zeros { + for vm%10 == 0 { + vr_is_trailing_zeros = vr_is_trailing_zeros && (last_removed_digit == 0) + last_removed_digit = byte(vr % 10) + vr /= 10 + vp /= 10 + vm /= 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 (~96.0%). Percentages below + // are relative to this. Loop iterations below (approximately): + // 0: 13.6%, 1: 70.7%, 2: 14.1%, 3: 1.39%, 4: 0.14%, 5+: 0.01% + for vp/10 > vm/10 { + last_removed_digit = byte(vr % 10) + 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_u32(vr == vm || last_removed_digit >= 5) + } + + return Dec32{m: out, e: e10 + removed} +} + +// f32_to_str return a string in scientific notation with max n_digit after the dot +pub fn f32_to_str(f f32, n_digit int) string { + mut u1 := Uf32{} + u1.f = f + u := u1.u + + neg := (u>>(mantbits32+expbits32)) != 0 + mant := u & ((u32(1)<> mantbits32) & ((u32(1)< 0 && out_len > n_digit { + out_len = n_digit+1 + } + + if neg { + buf[i]=`-` + i++ + } + + mut disp := 0 + if out_len <= 1 { + disp = 1 + } + + 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++ + } + + /* + x=0 + for x 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 mantbits64 { + return d, false + } + shift := mantbits64 - e + mant := i_mant | 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)<= 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)<= 100000000 { return 9 } + else if u >= 10000000 { return 8 } + else if u >= 1000000 { return 7 } + else if u >= 100000 { return 6 } + else if u >= 10000 { return 5 } + else if u >= 1000 { return 4 } + else if u >= 100 { return 3 } + else if u >= 10 { return 2 } + return 1 +} + +fn mul_shift_32(m u32, mul u64, ishift int) u32 { + assert ishift > 32 + + hi, lo := bits.mul_64(u64(m), mul) + shifted_sum := (lo >> u64(ishift)) + (hi << u64(64-ishift)) + assert1(shifted_sum <= math.max_u32, "shiftedSum <= math.max_u32") + return u32(shifted_sum) +} + +fn mul_pow5_invdiv_pow2(m u32, q u32, j int) u32 { + return mul_shift_32(m, pow5_inv_split_32[q], j) +} + +fn mul_pow5_div_pow2(m u32, i u32, j int) u32 { + return mul_shift_32(m, pow5_split_32[i], j) +} + +fn pow5_factor_32(i_v u32) u32 { + mut v := i_v + for n := u32(0); ; n++ { + q := v/5 + r := v%5 + if r != 0 { + return n + } + v = q + } + return v +} + +// multiple_of_power_of_five_32 reports whether v is divisible by 5^p. +fn multiple_of_power_of_five_32(v u32, p u32) bool { + return pow5_factor_32(v) >= p +} + +// multiple_of_power_of_two_32 reports whether v is divisible by 2^p. +fn multiple_of_power_of_two_32(v u32, p u32) bool { + return bits.trailing_zeros_32(v) >= p +} + +// log10_pow2 returns floor(log_10(2^e)). +fn log10_pow2(e int) u32 { + // The first value this approximation fails for is 2^1651 + // which is just greater than 10^297. + assert1(e >= 0, "e >= 0") + assert1(e <= 1650, "e <= 1650") + return (u32(e) * 78913) >> 18 +} + +// log10_pow5 returns floor(log_10(5^e)). +fn log10_pow5(e int) u32 { + // The first value this approximation fails for is 5^2621 + // which is just greater than 10^1832. + assert1(e >= 0, "e >= 0") + assert1(e <= 2620, "e <= 2620") + return (u32(e) * 732923) >> 20 +} + +// pow5_bits returns ceil(log_2(5^e)), or else 1 if e==0. +fn pow5_bits(e int) int { + // This approximation works up to the point that the multiplication + // overflows at e = 3529. If the multiplication were done in 64 bits, + // it would fail at 5^4004 which is just greater than 2^9297. + assert1(e >= 0, "e >= 0") + assert1(e <= 3528, "e <= 3528") + return int( ((u32(e)*1217359)>>19) + 1) +} + +/****************************************************************************** +* +* 64 bit functions +* +******************************************************************************/ +fn decimal_len_64(u u64) int { + // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLog10 + log2 := 64 - bits.leading_zeros_64(u) - 1 + t := (log2 + 1) * 1233 >> 12 + return t - bool_to_int(u < powers_of_10[t]) + 1 +} + +fn shift_right_128(v Uint128, shift int) u64 { + // The shift value is always modulo 64. + // In the current implementation of the 64-bit version + // of Ryu, the shift value is always < 64. + // (It is in the range [2, 59].) + // Check this here in case a future change requires larger shift + // values. In this case this function needs to be adjusted. + assert1(shift < 64, "shift < 64") + return (v.hi << u64(64 - shift)) | (v.lo >> u32(shift)) +} + +fn mul_shift_64(m u64, mul Uint128, shift int) u64 { + hihi, hilo := bits.mul_64(m, mul.hi) + lohi, _ := bits.mul_64(m, mul.lo) + mut sum := Uint128{hi: hihi, lo: lohi + hilo} + if sum.lo < lohi { + sum.hi++ // overflow + } + return shift_right_128(sum, shift-64) +} + +fn pow5_factor_64(v_i u64) u32 { + mut v := v_i + for n := u32(0); ; n++ { + q := v/5 + r := v%5 + if r != 0 { + return n + } + v = q + } + return u32(0) +} + +fn multiple_of_power_of_five_64(v u64, p u32) bool { + return pow5_factor_64(v) >= p +} + +fn multiple_of_power_of_two_64(v u64, p u32) bool { + return u32(bits.trailing_zeros_64(v)) >= p +} + +/****************************************************************************** +* +* f64 to string with string format +* +******************************************************************************/ + +// f32_to_str_l return a string with the f32 converted in a strign in decimal notation +pub fn f32_to_str_l(f f64) string { + return f64_to_str_l(f32(f)) +} + +// f64_to_str_l return a string with the f64 converted in a strign in decimal notation +pub fn f64_to_str_l(f f64) string { + s := ftoa.f64_to_str(f,18) + + // check for +inf -inf Nan + if s.len > 2 && (s[0] == `N` || s[1] == `i`) { + return s + } + + m_sgn_flag := false + mut sgn := 1 + mut b := [32]byte + mut d_pos := 1 + mut i := 0 + mut i1 := 0 + mut exp := 0 + mut exp_sgn := 1 + + // get sign and deciaml parts + for c in s { + if c == `-` { + sgn = -1 + i++ + } else if c == `+` { + sgn = 1 + i++ + } + else if c >= `0` && c <= `9` { + b[i1++] = c + i++ + } else if c == `.` { + if sgn > 0 { + d_pos = i + } else { + d_pos = i-1 + } + i++ + } else if c == `e` { + i++ + break + } else { + return "Float conversion error!!" + } + } + b[i1] = 0 + + // get exponent + if s[i] == `-` { + exp_sgn = -1 + i++ + } else if s[i] == `+` { + exp_sgn = 1 + i++ + } + for c in s[i..] { + exp = exp * 10 + int(c-`0`) + } + + // allocate exp+32 chars for the return string + mut res := [`0`].repeat(exp+32) // TODO: Slow!! is there other possibilities to allocate this? + mut r_i := 0 // result string buffer index + + //println("s:${sgn} b:${b[0]} es:${exp_sgn} exp:${exp}") + + if sgn == 1 { + if m_sgn_flag { + res[r_i++] = `+` + } + } else { + res[r_i++] = `-` + } + + i = 0 + if exp_sgn >= 0 { + for b[i] != 0 { + res[r_i++] = b[i] + i++ + if i >= d_pos && exp >= 0 { + if exp == 0 { + res[r_i++] = `.` + } + exp-- + } + } + for exp >= 0 { + res[r_i++] = `0` + exp-- + } + } else { + mut dot_p := true + for exp > 0 { + res[r_i++] = `0` + exp-- + if dot_p { + res[r_i++] = `.` + dot_p = false + } + } + for b[i] != 0 { + res[r_i++] = b[i] + i++ + } + } + res[r_i] = 0 + return tos(&res[0],r_i) +}