hare

[hare] The Hare programming language
git clone https://git.torresjrjr.com/hare.git
Log | Files | Refs | README | LICENSE

commit a65b9f6db75f600909f35cb1ee7b59d4cf836b6c
parent 239e16420ea1b809a413abc8ca897c31c8c175cb
Author: Joe Finney <me@spxtr.net>
Date:   Wed,  4 Oct 2023 12:52:59 -0700

strconv: implement ftosf.

For void precision, uses Ryu algorithm. For uint precision, uses a
multiprecision implementation heavily based on golang's strconv.

The rounding and encoding functions are highly nontrivial, largely
because I decided to add a SHOW_POINT flag.

Signed-off-by: Joe Finney <me@spxtr.net>

Diffstat:
Mmakefiles/freebsd.aarch64.mk | 4++--
Mmakefiles/freebsd.riscv64.mk | 4++--
Mmakefiles/freebsd.x86_64.mk | 4++--
Mmakefiles/linux.aarch64.mk | 4++--
Mmakefiles/linux.riscv64.mk | 4++--
Mmakefiles/linux.x86_64.mk | 4++--
Astrconv/+test/ftos_test.ha | 232+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Mstrconv/ftos.ha | 934++++++++++++++++++++++++-------------------------------------------------------
Astrconv/ftos_multiprecision.ha | 305++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Astrconv/ftos_ryu.ha | 503+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dstrconv/ftostof+test.ha | 33---------------------------------
11 files changed, 1332 insertions(+), 699 deletions(-)

diff --git a/makefiles/freebsd.aarch64.mk b/makefiles/freebsd.aarch64.mk @@ -117,8 +117,8 @@ $(HARECACHE)/os.ssa: $(os_ha) $(HARECACHE)/bufio.td $(HARECACHE)/encoding_utf8.t @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/os.ssa -t $(HARECACHE)/os.td.tmp -N os $(os_ha) -strconv_ha = strconv/ftos.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha -$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/math.td $(HARECACHE)/strings.td $(HARECACHE)/types.td +strconv_ha = strconv/ftos.ha strconv/ftos_multiprecision.ha strconv/ftos_ryu.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha +$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/io.td $(HARECACHE)/math.td $(HARECACHE)/memio.td $(HARECACHE)/strings.td $(HARECACHE)/types.td @mkdir -p -- "$(HARECACHE)" @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/strconv.ssa -t $(HARECACHE)/strconv.td.tmp -N strconv $(strconv_ha) diff --git a/makefiles/freebsd.riscv64.mk b/makefiles/freebsd.riscv64.mk @@ -117,8 +117,8 @@ $(HARECACHE)/os.ssa: $(os_ha) $(HARECACHE)/bufio.td $(HARECACHE)/encoding_utf8.t @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/os.ssa -t $(HARECACHE)/os.td.tmp -N os $(os_ha) -strconv_ha = strconv/ftos.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha -$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/math.td $(HARECACHE)/strings.td $(HARECACHE)/types.td +strconv_ha = strconv/ftos.ha strconv/ftos_multiprecision.ha strconv/ftos_ryu.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha +$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/io.td $(HARECACHE)/math.td $(HARECACHE)/memio.td $(HARECACHE)/strings.td $(HARECACHE)/types.td @mkdir -p -- "$(HARECACHE)" @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/strconv.ssa -t $(HARECACHE)/strconv.td.tmp -N strconv $(strconv_ha) diff --git a/makefiles/freebsd.x86_64.mk b/makefiles/freebsd.x86_64.mk @@ -117,8 +117,8 @@ $(HARECACHE)/os.ssa: $(os_ha) $(HARECACHE)/bufio.td $(HARECACHE)/encoding_utf8.t @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/os.ssa -t $(HARECACHE)/os.td.tmp -N os $(os_ha) -strconv_ha = strconv/ftos.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha -$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/math.td $(HARECACHE)/strings.td $(HARECACHE)/types.td +strconv_ha = strconv/ftos.ha strconv/ftos_multiprecision.ha strconv/ftos_ryu.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha +$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/io.td $(HARECACHE)/math.td $(HARECACHE)/memio.td $(HARECACHE)/strings.td $(HARECACHE)/types.td @mkdir -p -- "$(HARECACHE)" @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/strconv.ssa -t $(HARECACHE)/strconv.td.tmp -N strconv $(strconv_ha) diff --git a/makefiles/linux.aarch64.mk b/makefiles/linux.aarch64.mk @@ -135,8 +135,8 @@ $(HARECACHE)/os.ssa: $(os_ha) $(HARECACHE)/bufio.td $(HARECACHE)/encoding_utf8.t @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/os.ssa -t $(HARECACHE)/os.td.tmp -N os $(os_ha) -strconv_ha = strconv/ftos.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha -$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/math.td $(HARECACHE)/strings.td $(HARECACHE)/types.td +strconv_ha = strconv/ftos.ha strconv/ftos_multiprecision.ha strconv/ftos_ryu.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha +$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/io.td $(HARECACHE)/math.td $(HARECACHE)/memio.td $(HARECACHE)/strings.td $(HARECACHE)/types.td @mkdir -p -- "$(HARECACHE)" @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/strconv.ssa -t $(HARECACHE)/strconv.td.tmp -N strconv $(strconv_ha) diff --git a/makefiles/linux.riscv64.mk b/makefiles/linux.riscv64.mk @@ -135,8 +135,8 @@ $(HARECACHE)/os.ssa: $(os_ha) $(HARECACHE)/bufio.td $(HARECACHE)/encoding_utf8.t @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/os.ssa -t $(HARECACHE)/os.td.tmp -N os $(os_ha) -strconv_ha = strconv/ftos.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha -$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/math.td $(HARECACHE)/strings.td $(HARECACHE)/types.td +strconv_ha = strconv/ftos.ha strconv/ftos_multiprecision.ha strconv/ftos_ryu.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha +$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/io.td $(HARECACHE)/math.td $(HARECACHE)/memio.td $(HARECACHE)/strings.td $(HARECACHE)/types.td @mkdir -p -- "$(HARECACHE)" @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/strconv.ssa -t $(HARECACHE)/strconv.td.tmp -N strconv $(strconv_ha) diff --git a/makefiles/linux.x86_64.mk b/makefiles/linux.x86_64.mk @@ -135,8 +135,8 @@ $(HARECACHE)/os.ssa: $(os_ha) $(HARECACHE)/bufio.td $(HARECACHE)/encoding_utf8.t @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/os.ssa -t $(HARECACHE)/os.td.tmp -N os $(os_ha) -strconv_ha = strconv/ftos.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha -$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/math.td $(HARECACHE)/strings.td $(HARECACHE)/types.td +strconv_ha = strconv/ftos.ha strconv/ftos_multiprecision.ha strconv/ftos_ryu.ha strconv/itos.ha strconv/numeric.ha strconv/stof.ha strconv/stof_data.ha strconv/stoi.ha strconv/stou.ha strconv/types.ha strconv/utos.ha +$(HARECACHE)/strconv.ssa: $(strconv_ha) $(HARECACHE)/ascii.td $(HARECACHE)/bytes.td $(HARECACHE)/encoding_utf8.td $(HARECACHE)/io.td $(HARECACHE)/math.td $(HARECACHE)/memio.td $(HARECACHE)/strings.td $(HARECACHE)/types.td @mkdir -p -- "$(HARECACHE)" @printf 'HAREC\t%s\n' "$@" @$(TDENV) $(HAREC) $(HARECFLAGS) -o $(HARECACHE)/strconv.ssa -t $(HARECACHE)/strconv.td.tmp -N strconv $(strconv_ha) diff --git a/strconv/+test/ftos_test.ha b/strconv/+test/ftos_test.ha @@ -0,0 +1,232 @@ +// SPDX-License-Identifier: MPL-2.0 +// (c) Hare authors <https://harelang.org> + +use math; + +@test fn ftosf() void = { + // These tests should pass for both f32 and f64. + const tcs: [](f64, ffmt, (void | uint), fflags, str) = [ + // First test special values + (1.0 / 0.0, ffmt::G, void, 0, "infinity"), + (1.0 / 0.0, ffmt::G, void, fflags::SHOW_POS, "+infinity"), + (-1.0 / 0.0, ffmt::F, void, 0, "-infinity"), + (-1.0 / 0.0, ffmt::E, void, fflags::UPPERCASE, "-INFINITY"), + (0.0 / 0.0, ffmt::G, void, 0, "nan"), + (0.0 / 0.0, ffmt::E, void, fflags::UPPERCASE, "NAN"), + + // Then a million tests for zero. + (0.0, ffmt::E, void, 0, "0e0"), + (0.0, ffmt::E, void, fflags::SHOW_TWO_EXP_DIGITS, "0e00"), + (0.0, ffmt::E, void, fflags::UPPER_EXP, "0E0"), + (0.0, ffmt::E, void, fflags::SHOW_POS_EXP, "0e+0"), + (0.0, ffmt::E, void, fflags::SHOW_POINT, "0.0e0"), + (0.0, ffmt::F, void, 0, "0"), + (0.0, ffmt::F, void, fflags::SHOW_POINT, "0.0"), + (0.0, ffmt::G, void, 0, "0"), + (-0.0, ffmt::G, void, fflags::SHOW_POS, "-0"), + (0.0, ffmt::G, void, fflags::SHOW_POS, "+0"), + (0.0, ffmt::G, void, fflags::SHOW_POINT, "0.0"), + + // ... e and f do not cut trailing zeros + (0.0, ffmt::E, 0, 0, "0e0"), + (0.0, ffmt::E, 1, 0, "0.0e0"), + (0.0, ffmt::E, 2, 0, "0.00e0"), + (0.0, ffmt::F, 0, 0, "0"), + (0.0, ffmt::F, 1, 0, "0.0"), + (0.0, ffmt::F, 2, 0, "0.00"), + // ... g cuts trailing zeros + (0.0, ffmt::G, 0, 0, "0"), + (0.0, ffmt::G, 1, 0, "0"), + (0.0, ffmt::G, 2, 0, "0"), + + // ... SHOW_POINT only changes precision 0 + (0.0, ffmt::E, 0, fflags::SHOW_POINT, "0.0e0"), + (0.0, ffmt::E, 1, fflags::SHOW_POINT, "0.0e0"), + (0.0, ffmt::E, 2, fflags::SHOW_POINT, "0.00e0"), + (0.0, ffmt::F, 0, fflags::SHOW_POINT, "0.0"), + (0.0, ffmt::F, 1, fflags::SHOW_POINT, "0.0"), + (0.0, ffmt::F, 2, fflags::SHOW_POINT, "0.00"), + // ... g with SHOW_POINT only has the one extra 0 + (0.0, ffmt::G, 0, fflags::SHOW_POINT, "0.0"), + (0.0, ffmt::G, 1, fflags::SHOW_POINT, "0.0"), + (0.0, ffmt::G, 2, fflags::SHOW_POINT, "0.0"), + + // Now we can test actual numbers. + (10.0, ffmt::F, void, 0, "10"), + (1.0, ffmt::F, void, 0, "1"), + (1.1, ffmt::F, void, 0, "1.1"), + (13.37, ffmt::G, void, 0, "13.37"), + (0.3, ffmt::F, void, 0, "0.3"), + (0.0031415, ffmt::F, void, 0, "0.0031415"), + (-6345.972, ffmt::F, void, 0, "-6345.972"), + (1.414, ffmt::F, void, 0, "1.414"), + (1000000.0e9, ffmt::F, void, 0, "1000000000000000"), + (10.0, ffmt::E, void, 0, "1e1"), + (10.0, ffmt::E, void, fflags::SHOW_TWO_EXP_DIGITS, "1e01"), + (10.0, ffmt::E, void, fflags::UPPER_EXP, "1E1"), + (10.0, ffmt::E, void, fflags::SHOW_POS_EXP, "1e+1"), + (0.1, ffmt::E, void, fflags::SHOW_POS_EXP, "1e-1"), + (1.0, ffmt::E, void, 0, "1e0"), + (0.3, ffmt::E, void, 0, "3e-1"), + (0.0031415, ffmt::E, void, 0, "3.1415e-3"), + (0.12345, ffmt::E, void, 0, "1.2345e-1"), + + // ... g is shortest + (12345.0, ffmt::G, void, 0, "12345"), + (10000.0, ffmt::G, void, 0, "1e4"), + (11000.0, ffmt::G, void, 0, "1.1e4"), + (1000.0, ffmt::G, void, 0, "1e3"), + (1100.0, ffmt::G, void, 0, "1100"), + (100.0, ffmt::G, void, 0, "100"), + (10.0, ffmt::G, void, 0, "10"), + (1.0, ffmt::G, void, 0, "1"), + (0.1, ffmt::G, void, 0, "0.1"), + (0.01, ffmt::G, void, 0, "0.01"), + (0.011, ffmt::G, void, 0, "0.011"), + (0.001, ffmt::G, void, 0, "1e-3"), // one shorter than f + (0.0011, ffmt::G, void, 0, "1.1e-3"), // same length as f + (0.0001, ffmt::G, void, 0, "1e-4"), + + // ... fixed precision stuff + (0.5, ffmt::F, 0, 0, "0"), + (1.0 / 3.0, ffmt::F, 2, 0, "0.33"), + (1.0 / 3.0, ffmt::F, 1, 0, "0.3"), + (1.0 / 3.0, ffmt::F, 0, 0, "0"), + (1.0 / 3.0, ffmt::F, 0, fflags::SHOW_POINT, "0.3"), + (2.0 / 3.0, ffmt::F, 2, 0, "0.67"), + (2.0 / 3.0, ffmt::F, 1, 0, "0.7"), + (2.0 / 3.0, ffmt::F, 0, 0, "1"), + (2.0 / 3.0, ffmt::F, 0, fflags::SHOW_POINT, "0.7"), + (2.0 / 30.0, ffmt::F, 5, 0, "0.06667"), + (2.0 / 30.0, ffmt::F, 2, 0, "0.07"), + (2.0 / 30.0, ffmt::F, 1, 0, "0.1"), + (2.0 / 30.0, ffmt::F, 1, fflags::SHOW_POINT, "0.1"), + (200.0 / 3.0, ffmt::F, 4, 0, "66.6667"), + (100.0 / 3.0, ffmt::F, 4, 0, "33.3333"), + (100.0 / 3.0, ffmt::F, 0, 0, "33"), + (100.0 / 3.0, ffmt::F, 0, fflags::SHOW_POINT, "33.3"), + (0.00001, ffmt::F, 1, 0, "0.0"), + (0.001, ffmt::F, 2, 0, "0.00"), + (0.006, ffmt::F, 2, 0, "0.01"), + (0.001, ffmt::F, 6, 0, "0.001000"), + (1.0, ffmt::F, 6, 0, "1.000000"), + (100.0, ffmt::F, 6, 0, "100.000000"), + + // ... scientific notation stuff + (460.0, ffmt::E, 2, 0, "4.60e2"), + (1.0 / 3.0, ffmt::E, 2, 0, "3.33e-1"), + (1.0 / 3.0, ffmt::E, 1, 0, "3.3e-1"), + (1.0 / 3.0, ffmt::E, 0, 0, "3e-1"), + (1.0 / 3.0, ffmt::E, 0, fflags::SHOW_POINT, "3.3e-1"), + (2.0 / 3.0, ffmt::E, 2, 0, "6.67e-1"), + (2.0 / 3.0, ffmt::E, 1, 0, "6.7e-1"), + (2.0 / 3.0, ffmt::E, 0, 0, "7e-1"), + (2.0 / 3.0, ffmt::E, 0, fflags::SHOW_POINT, "6.7e-1"), + (2.0 / 30.0, ffmt::E, 5, 0, "6.66667e-2"), + (2.0 / 30.0, ffmt::E, 2, 0, "6.67e-2"), + (2.0 / 30.0, ffmt::E, 0, 0, "7e-2"), + (2.0 / 30.0, ffmt::E, 0, fflags::SHOW_POINT, "6.7e-2"), + (200.0 / 3.0, ffmt::E, 5, 0, "6.66667e1"), + (100.0 / 3.0, ffmt::E, 5, 0, "3.33333e1"), + (100.0 / 3.0, ffmt::E, 0, 0, "3e1"), + (100.0 / 3.0, ffmt::E, 0, fflags::SHOW_POINT, "3.3e1"), + (0.001, ffmt::E, 2, 0, "1.00e-3"), + (1.0, ffmt::E, 6, 0, "1.000000e0"), + (100.0, ffmt::E, 6, 0, "1.000000e2"), + + // ... and G. The behavior with SHOW_POINT is gnarly. + (1.0 / 3.0, ffmt::G, 2, 0, "0.33"), + (1.0 / 3.0, ffmt::G, 0, 0, "0.3"), + (0.01, ffmt::G, void, fflags::SHOW_POINT, "0.01"), + (1.0, ffmt::G, void, fflags::SHOW_POINT, "1.0"), + (0.0001, ffmt::G, 0, 0, "1e-4"), + (0.001, ffmt::G, 0, 0, "1e-3"), + (0.00123, ffmt::G, 2, 0, "1.2e-3"), + (0.01, ffmt::G, 5, 0, "0.01"), // trim trailing zeros + (0.1, ffmt::G, 5, 0, "0.1"), + (1.0, ffmt::G, 0, 0, "1"), + (10.0, ffmt::G, 0, 0, "10"), + (120.0, ffmt::G, 2, 0, "120"), + (12000.0, ffmt::G, 2, 0, "1.2e4"), + (0.0001, ffmt::G, 0, fflags::SHOW_POINT, "1.0e-4"), + (0.001, ffmt::G, 0, fflags::SHOW_POINT, "1.0e-3"), + (0.01, ffmt::G, 0, fflags::SHOW_POINT, "0.01"), + (0.1, ffmt::G, 0, fflags::SHOW_POINT, "0.1"), + (1.0, ffmt::G, 0, fflags::SHOW_POINT, "1.0"), + (10.0, ffmt::G, 0, fflags::SHOW_POINT, "10.0"), + (100.0, ffmt::G, 0, fflags::SHOW_POINT, "100.0"), + (1000.0, ffmt::G, 0, fflags::SHOW_POINT, "1.0e3"), + (0.0123, ffmt::G, 2, fflags::SHOW_POINT, "0.012"), + (0.0123, ffmt::G, 5, fflags::SHOW_POINT, "0.0123"), + ]; + for (let i = 0z; i < len(tcs); i += 1) { + const res64 = ftosf(tcs[i].0, tcs[i].1, tcs[i].2, tcs[i].3); + //defer free(res64); + assert(res64 == tcs[i].4, res64); + if (tcs[i].2 is void) { + // void precision should guarantee that it parses back + // to the original number. + const back = stof64(res64)!; + assert(math::isnan(back) == math::isnan(tcs[i].0)); + if (!math::isnan(back)) + assert(back == tcs[i].0); + }; + + const res32 = ftosf(tcs[i].0: f32, tcs[i].1, + tcs[i].2, tcs[i].3); + defer free(res32); + assert(res32 == tcs[i].4); + if (tcs[i].2 is void) { + const back = stof32(res32)!; + assert(math::isnan(back) == math::isnan(tcs[i].0)); + if (!math::isnan(back)) + assert(back == tcs[i].0: f32); + }; + }; + // These tests will only pass for f64 + const tcsf64: [](f64, ffmt, (void | uint), fflags, str) = [ + (9007199254740991.0, ffmt::F, void, 0, "9007199254740991"), + (90071992547409915.0, ffmt::F, void, 0, "90071992547409920"), + (90071992547409925.0, ffmt::F, void, 0, "90071992547409920"), + (math::F64_MIN, ffmt::E, void, 0, "5e-324"), + (math::F64_MIN, ffmt::E, void, fflags::SHOW_TWO_EXP_DIGITS, "5e-324"), + (-math::F64_MIN, ffmt::E, void, 0, "-5e-324"), + (math::F64_MIN_NORMAL, ffmt::E, void, 0, "2.2250738585072014e-308"), + (math::F64_MAX_NORMAL, ffmt::E, void, 0, "1.7976931348623157e308"), + (math::F64_MIN, ffmt::E, 2, 0, "4.94e-324"), + (math::F64_MIN_NORMAL, ffmt::E, 0, 0, "2e-308"), + (math::F64_MAX_NORMAL, ffmt::E, 3, 0, "1.798e308"), + ]; + for (let i = 0z; i < len(tcsf64); i += 1) { + const res64 = ftosf(tcsf64[i].0, tcsf64[i].1, + tcsf64[i].2, tcsf64[i].3); + defer free(res64); + assert(res64 == tcsf64[i].4); + }; + // These tests will only pass for f32 + const tcsf32: [](f32, ffmt, (void | uint), fflags, str) = [ + (math::F32_MIN, ffmt::G, void, 0, "1e-45"), + (math::F32_MIN_NORMAL, ffmt::G, void, 0, "1.1754944e-38"), + (math::F32_MAX_NORMAL, ffmt::G, void, 0, "3.4028235e38"), + ]; + for (let i = 0z; i < len(tcsf32); i += 1) { + const res32 = ftosf(tcsf32[i].0, tcsf32[i].1, + tcsf32[i].2, tcsf32[i].3); + defer free(res32); + assert(res32 == tcsf32[i].4); + }; + // Just make sure we can generate big numbers without breaking anything. + const tcslen: [](f64, ffmt, (void | uint), fflags, size) = [ + (9007199254740991.0, ffmt::F, void, 0, 16), + (-math::F64_MIN, ffmt::E, 100, 0, 108), + (1.0, ffmt::F, 1000, 0, 1002), + (2.22507385850720088902458687609E-308, ffmt::F, 1000, 0, 1002), + ]; + for (let i = 0z; i < len(tcslen); i += 1) { + const res64 = ftosf(tcslen[i].0, tcslen[i].1, + tcslen[i].2, tcslen[i].3); + defer free(res64); + assert(len(res64) == tcslen[i].4); + }; + assert(f64tos(13.37) == "13.37"); +}; diff --git a/strconv/ftos.ha b/strconv/ftos.ha @@ -1,265 +1,50 @@ // SPDX-License-Identifier: MPL-2.0 // (c) Hare authors <https://harelang.org> -// Using Ryū: fast float-to-string conversion algorithm by Ulf Adams. -// https://doi.org/10.1145/3192366.3192369 -// This Hare implementation is translated from the original -// C implementation here: https://github.com/ulfjack/ryu +// Uses Ryū for shortest, falls back to multiprecision for fixed precision. +use ascii; +use io; use math; +use memio; +use strings; use types; -type r128 = struct { - hi: u64, - lo: u64, +export type ffmt = enum { + // Scientific notation. Consists of a number in [1, 10), an 'e' (or 'E', + // if UPPER_EXP flag is present), then an exponent. + E, + // Fixed-point notation. + F, + // General format. Uses whichever of E and F is shortest, not accounting + // for flags. + G, }; -// TODO: use 128-bit integers when implemented -fn u128mul(a: u64, b: u64) r128 = { - const a0 = a: u32: u64, a1 = a >> 32; - const b0 = b: u32: u64, b1 = b >> 32; - const p00 = a0 * b0, p01 = a0 * b1, p10 = a1 * b0, p11 = a1 * b1; - const p00_lo = p00: u32: u64, p00_hi = p00 >> 32; - const mid1 = p10 + p00_hi; - const mid1_lo = mid1: u32: u64, mid1_hi = mid1 >> 32; - const mid2 = p01 + mid1_lo; - const mid2_lo = mid2: u32: u64, mid2_hi = mid2 >> 32; - const r_hi = p11 + mid1_hi + mid2_hi; - const r_lo = (mid2_lo << 32) | p00_lo; - return r128 { hi = r_hi, lo = r_lo }; +export type fflags = enum uint { + // Use a sign for both positive and negative numbers. + SHOW_POS = 1 << 0, + // Include at least one decimal digit. + SHOW_POINT = 1 << 1, + // Uppercase INFINITY and NAN. + UPPERCASE = 1 << 2, + // Uppercase exponent symbols E and P rather than e and p. + UPPER_EXP = 1 << 3, + // Use a sign for both positive and negative exponents. + SHOW_POS_EXP = 1 << 4, + // Show at least two digits of the exponent. + SHOW_TWO_EXP_DIGITS = 1 << 5, }; -// TODO: Same as above -fn u128rshift(lo: u64, hi: u64, s: u32) u64 = { - assert(0 <= s); - assert(s <= 64); - return (hi << (64 - s)) | (lo >> s); -}; - -fn pow5fac(value: u64) u32 = { - const m_inv_5: u64 = 14757395258967641293; // 5 * m_inv_5 = 1 (mod 2^64) - const n_div_5: u64 = 3689348814741910323; - let count: u32 = 0; - for (true) { - assert(value != 0); - value *= m_inv_5; - if (value > n_div_5) break; - count += 1; - }; - return count; -}; - -fn pow5fac32(value: u32) u32 = { - let count: u32 = 0; - for (true) { - assert(value != 0); - const q = value / 5, r = value % 5; - if (r != 0) break; - value = q; - count += 1; - }; - return count; -}; - -fn ibool(b: bool) u8 = if (b) 1 else 0; - -fn pow5multiple(v: u64, p: u32) bool = pow5fac(v) >= p; -fn pow5multiple32(v: u32, p: u32) bool = pow5fac32(v) >= p; - -fn pow2multiple(v: u64, p: u32) bool = { - assert(v > 0); - assert(p < 64); - return (v & ((1u64 << p) - 1)) == 0; -}; - -fn pow2multiple32(v: u32, p: u32) bool = { - assert(v > 0); - assert(p < 32); - return (v & ((1u32 << p) - 1)) == 0; -}; - -fn mulshift64(m: u64, mul: (u64, u64), j: u32) u64 = { - // m is maximum 55 bits - let r0 = u128mul(m, mul.0), r1 = u128mul(m, mul.1); - const sum = r1.lo + r0.hi; - r1.hi += ibool(sum < r0.hi); - return u128rshift(sum, r1.hi, j - 64); -}; - -fn mulshiftall64(m: u64, mul: (u64, u64), j: i32, mm_shift: u32) (u64, u64, u64) = { - m <<= 1; - const r0 = u128mul(m, mul.0), r1 = u128mul(m, mul.1); - const lo = r0.lo, tmp = r0.hi, mid = tmp + r1.lo; - const hi = r1.hi + ibool(mid < tmp); - const lo2 = lo + mul.0; - const mid2 = mid + mul.1 + ibool(lo2 < lo); - const hi2 = hi + ibool(mid2 < mid); - const v_plus = u128rshift(mid2, hi2, (j - 64 - 1): u32); - const v_minus = if (mm_shift == 1) { - const lo3 = lo - mul.0; - const mid3 = mid - mul.1 - ibool(lo3 > lo); - const hi3 = hi - ibool(mid3 > mid); - yield u128rshift(mid3, hi3, (j - 64 - 1): u32); - } else { - const lo3 = lo + lo; - const mid3 = mid + mid + ibool(lo3 < lo); - const hi3 = hi + hi + ibool(mid3 < mid); - const lo4 = lo3 - mul.0; - const mid4 = mid3 - mul.1 - ibool(lo4 > lo3); - const hi4 = hi3 - ibool(mid4 > mid3); - yield u128rshift(mid4, hi4, (j - 64): u32); - }; - const v_rounded = u128rshift(mid, hi, (j - 64 - 1): u32); - return (v_plus, v_rounded, v_minus); -}; - -fn mulshift32(m: u32, a: u64, s: u32) u32 = { - assert(s > 32); - const a_lo = a: u32: u64, a_hi = a >> 32; - const b0 = m * a_lo, b1 = m * a_hi; - const sum = (b0 >> 32) + b1, ss = sum >> (s - 32); - assert(ss <= types::U32_MAX); - return ss: u32; -}; - -fn mulpow5inv_divpow2(m: u32, q: u32, j: i32) u32 = { - const pow5 = f64computeinvpow5(q); - return mulshift32(m, pow5.1 + 1, j: u32); -}; - -fn mulpow5_divpow2(m: u32, i: u32, j: i32) u32 = { - const pow5 = f64computepow5(i); - return mulshift32(m, pow5.1, j: u32); -}; - -fn log2pow5(e: u32) u32 = { - assert(e <= 3528); - return ((e * 1217359) >> 19); -}; - -fn ceil_log2pow5(e: u32) u32 = log2pow5(e) + 1; - -fn pow5bits(e: u32) u32 = ceil_log2pow5(e); - -fn log10pow2(e: u32) u32 = { - assert(e <= 1650); - return ((e * 78913) >> 18); -}; +// Just for convenience... inline functions when? +fn ffpos(f: fflags) bool = f & fflags::SHOW_POS != 0; +fn ffpoint(f: fflags) bool = f & fflags::SHOW_POINT != 0; +fn ffcaps(f: fflags) bool = f & fflags::UPPERCASE != 0; +fn ffcaps_exp(f: fflags) bool = f & fflags::UPPER_EXP != 0; +fn ffpos_exp(f: fflags) bool = f & fflags::SHOW_POS_EXP != 0; +fn fftwodigs(f: fflags) bool = f & fflags::SHOW_TWO_EXP_DIGITS != 0; -fn log10pow5(e: u32) u32 = { - assert(e <= 2620); - return ((e * 732923) >> 20); -}; - -def F64_POW5_INV_BITCOUNT: u8 = 125; -def F64_POW5_BITCOUNT: u8 = 125; - -def F32_POW5_INV_BITCOUNT: u8 = F64_POW5_INV_BITCOUNT - 64; -def F32_POW5_BITCOUNT: u8 = F64_POW5_BITCOUNT - 64; - -const F64_POW5_INV_SPLIT2: [15][2]u64 = [ - [1, 2305843009213693952], - [5955668970331000884, 1784059615882449851], - [8982663654677661702, 1380349269358112757], - [7286864317269821294, 2135987035920910082], - [7005857020398200553, 1652639921975621497], - [17965325103354776697, 1278668206209430417], - [8928596168509315048, 1978643211784836272], - [10075671573058298858, 1530901034580419511], - [597001226353042382, 1184477304306571148], - [1527430471115325346, 1832889850782397517], - [12533209867169019542, 1418129833677084982], - [5577825024675947042, 2194449627517475473], - [11006974540203867551, 1697873161311732311], - [10313493231639821582, 1313665730009899186], - [12701016819766672773, 2032799256770390445], -]; - -const POW5_INV_OFFSETS: [19]u32 = [ - 0x54544554, 0x04055545, 0x10041000, 0x00400414, 0x40010000, 0x41155555, - 0x00000454, 0x00010044, 0x40000000, 0x44000041, 0x50454450, 0x55550054, - 0x51655554, 0x40004000, 0x01000001, 0x00010500, 0x51515411, 0x05555554, - 0x00000000 -]; - -const F64_POW5_SPLIT2: [13][2]u64 = [ - [0, 1152921504606846976], - [0, 1490116119384765625], - [1032610780636961552, 1925929944387235853], - [7910200175544436838, 1244603055572228341], - [16941905809032713930, 1608611746708759036], - [13024893955298202172, 2079081953128979843], - [6607496772837067824, 1343575221513417750], - [17332926989895652603, 1736530273035216783], - [13037379183483547984, 2244412773384604712], - [1605989338741628675, 1450417759929778918], - [9630225068416591280, 1874621017369538693], - [665883850346957067, 1211445438634777304], - [14931890668723713708, 1565756531257009982] -]; - -const POW5_OFFSETS: [21]u32 = [ - 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x40000000, 0x59695995, - 0x55545555, 0x56555515, 0x41150504, 0x40555410, 0x44555145, 0x44504540, - 0x45555550, 0x40004000, 0x96440440, 0x55565565, 0x54454045, 0x40154151, - 0x55559155, 0x51405555, 0x00000105 -]; - -def POW5_TABLE_SZ: u8 = 26; - -const POW5_TABLE: [POW5_TABLE_SZ]u64 = [ - 1u64, 5u64, 25u64, 125u64, 625u64, 3125u64, 15625u64, 78125u64, - 390625u64, 1953125u64, 9765625u64, 48828125u64, 244140625u64, - 1220703125u64, 6103515625u64, 30517578125u64, 152587890625u64, - 762939453125u64, 3814697265625u64, 19073486328125u64, 95367431640625u64, - 476837158203125u64, 2384185791015625u64, 11920928955078125u64, - 59604644775390625u64, 298023223876953125u64 //, 1490116119384765625u64 -]; - -fn f64computeinvpow5(i: u32) (u64, u64) = { - const base = ((i + POW5_TABLE_SZ - 1) / POW5_TABLE_SZ): u32; - const base2 = base * POW5_TABLE_SZ; - const mul = F64_POW5_INV_SPLIT2[base]; - const off = base2 - i; - if (off == 0) { - return (mul[0], mul[1]); - }; - const m = POW5_TABLE[off]; - const r1 = u128mul(m, mul[1]), r0 = u128mul(m, mul[0] - 1); - let high1 = r1.hi, low1 = r1.lo, high0 = r0.hi, low0 = r0.lo; - const sum = high0 + low1; - if (sum < high0) { - high1 += 1; - }; - const delta = pow5bits(base2) - pow5bits(i); - const res0 = u128rshift(low0, sum, delta) + 1 + - ((POW5_INV_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3); - const res1 = u128rshift(sum, high1, delta); - return (res0, res1); -}; - -fn f64computepow5(i: u32) (u64, u64) = { - const base = i / POW5_TABLE_SZ, base2 = base * POW5_TABLE_SZ; - const mul = F64_POW5_SPLIT2[base]; - const off = i - base2; - if (off == 0) { - return (mul[0], mul[1]); - }; - const m = POW5_TABLE[off]; - const r1 = u128mul(m, mul[1]), r0 = u128mul(m, mul[0]); - let high1 = r1.hi, low1 = r1.lo, high0 = r0.hi, low0 = r0.lo; - const sum = high0 + low1; - if (sum < high0) { - high1 += 1; - }; - const delta = pow5bits(i) - pow5bits(base2); - const res0 = u128rshift(low0, sum, delta) + - ((POW5_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3); - const res1 = u128rshift(sum, high1, delta); - return (res0, res1); -}; - -fn declen(n: u64) u8 = { +fn declen(n: u64) uint = { assert(n <= 1e17); return if (n >= 1e17) 18 else if (n >= 1e16) 17 @@ -281,341 +66,270 @@ fn declen(n: u64) u8 = { else 1; }; -type decf64 = struct { - mantissa: u64, - exponent: i32, +fn writestr(h: io::handle, s: str) (void | io::error) = { + io::writeall(h, strings::toutf8(s))?; }; -fn f64todecf64(mantissa: u64, exponent: u32) decf64 = { - let e2 = (math::F64_EXPONENT_BIAS + math::F64_MANTISSA_BITS + 2): i32; - let m2: u64 = 0; - if (exponent == 0) { - e2 = 1 - e2; - m2 = mantissa; - } else { - e2 = (exponent: i32) - e2; - m2 = (1u64 << math::F64_MANTISSA_BITS) | mantissa; - }; - const accept_bounds = (m2 & 1) == 0; - const mv = 4 * m2; - const mm_shift = ibool(mantissa != 0 || exponent <= 1); - let vp: u64 = 0, vr: u64 = 0, vm: u64 = 0; - let e10: i32 = 0; - let vm_trailing_zeros = false, vr_trailing_zeros = false; - if (e2 >= 0) { - const q = log10pow2(e2: u32) - ibool(e2 > 3); - e10 = q: i32; - const k = F64_POW5_INV_BITCOUNT + pow5bits(q) - 1; - const i = -e2 + (q + k): i32; - let pow5 = f64computeinvpow5(q); - const res = mulshiftall64(m2, pow5, i, mm_shift); - vp = res.0; vr = res.1; vm = res.2; - if (q <= 21) { - if ((mv - 5 * (mv / 5)) == 0) { - vr_trailing_zeros = pow5multiple(mv, q); - } else if (accept_bounds) { - vm_trailing_zeros = pow5multiple(mv - 1 - mm_shift, q); - } else { - vp -= ibool(pow5multiple(mv + 2, q)); +// XXX: this can likely be dedup'd with the other encode functions. +fn encode_zero( + h: io::handle, + f: ffmt, + prec: (void | uint), + flag: fflags, +) (void | io::error) = { + memio::appendrune(h, '0')?; + let hasdec = false; + match (prec) { + case void => void; + case let u: uint => + if (u > 0 && f != ffmt::G) { + memio::appendrune(h, '.')?; + for (let i = 0u; i < u; i += 1) { + memio::appendrune(h, '0')?; }; - }; - } else { - const q = log10pow5((-e2): u32) - ibool(-e2 > 1); - e10 = e2 + (q: i32); - const i = -e2 - (q: i32); - const k = pow5bits(i: u32): i32 - F64_POW5_BITCOUNT: i32; - const j = (q: i32) - k; - let pow5 = f64computepow5(i: u32); - const res = mulshiftall64(m2, pow5, j, mm_shift); - vp = res.0; vr = res.1; vm = res.2; - if (q <= 1) { - vr_trailing_zeros = true; - if (accept_bounds) { - vm_trailing_zeros = mm_shift == 1; - } else { - vp -= 1; - }; - } else if (q < 63) { - vr_trailing_zeros = pow2multiple(mv, q); + hasdec = true; }; }; - let removed: i32 = 0, last_removed_digit: u8 = 0; - let output: u64 = 0; - if (vm_trailing_zeros || vr_trailing_zeros) { - for (true) { - const vpby10 = vp / 10, vmby10 = vm / 10; - if (vpby10 <= vmby10) break; - const vmmod10 = (vm: u32) - 10 * (vmby10: u32); - const vrby10 = vr / 10; - const vrmod10 = (vr: u32) - 10 * (vrby10: u32); - vm_trailing_zeros &&= vmmod10 == 0; - vr_trailing_zeros &&= last_removed_digit == 0; - last_removed_digit = vrmod10: u8; - vr = vrby10; vp = vpby10; vm = vmby10; - removed += 1; - }; - if (vm_trailing_zeros) { - for (true) { - const vmby10 = vm / 10; - const vmmod10 = (vm: u32) - 10 * (vmby10: u32); - if (vmmod10 != 0) break; - const vpby10 = vp / 10, vrby10 = vr / 10; - const vrmod10 = (vr: u32) - 10 * (vrby10: u32); - vr_trailing_zeros &&= last_removed_digit == 0; - last_removed_digit = vrmod10: u8; - vr = vrby10; vp = vpby10; vm = vmby10; - removed += 1; + if (!hasdec && ffpoint(flag)) { + memio::appendrune(h, '.')?; + memio::appendrune(h, '0')?; + }; + if (f == ffmt::E) { + memio::appendrune(h, if (ffcaps_exp(flag)) 'E' else 'e')?; + if (ffpos_exp(flag)) memio::appendrune(h, '+')?; + memio::appendrune(h, '0')?; + if (fftwodigs(flag)) memio::appendrune(h, '0')?; + }; +}; + +fn encode_f_mp( + m: *mp, + h: io::handle, + f: ffmt, + prec: (void | uint), + flag: fflags, +) (void | io::error) = { + // we will loop from lo <= i < hi, printing either zeros or a digit. + // lo is simple, but hi depends intricately on f, prec, and the + // SHOW_POINT flag. + const lo = if (m.dp <= 0) m.dp - 1 else 0; + let hi = match (prec) { + case void => + yield if (m.nd: int > m.dp) m.nd: int else m.dp; + case let u: uint => + yield if (m.dp <= 0) lo + u: int + 1 else m.dp + u: int; + }; + // ffmt::G: we need to remove trailing zeros + if (f == ffmt::G) { + // first, make sure we include at least prec digits + if (prec is uint) { + const p = prec as uint; + if (m.dp <= 0 && hi < p: int) { + hi = p: int; }; }; - if (vr_trailing_zeros && last_removed_digit == 5 && (vr & 1 == 0)) { - // round to even - last_removed_digit = 4; + // then, cut back to the decimal point or nd + if (hi > m.nd: int && m.dp <= 0) { + hi = m.nd: int; + } else if (hi > m.dp && m.dp > 0) { + hi = if (m.nd: int > m.dp) m.nd: int else m.dp; }; - output = vr + ibool((vr == vm && - (!accept_bounds || !vm_trailing_zeros)) || - last_removed_digit >= 5); - } else { - let round_up = false; - const vpby100 = vp / 100, vmby100 = vm / 100; - if (vpby100 > vmby100) { - const vrby100 = vr / 100; - const vrmod100 = (vr: u32) - 100 * (vrby100: u32); - round_up = vrmod100 >= 50; - vr = vrby100; vp = vpby100; vm = vmby100; - removed += 2; + }; + // SHOW_POINT: we need to go at least one past the decimal + if (ffpoint(flag) && hi <= m.dp) { + hi = m.dp + 1; + }; + for (let i = lo; i < hi; i += 1) { + if (i == m.dp) { + memio::appendrune(h, '.')?; }; - for (true) { - const vmby10 = vm / 10, vpby10 = vp / 10; - if (vpby10 <= vmby10) break; - const vrby10 = vr / 10; - const vrmod10 = (vr: u32) - 10 * (vrby10: u32); - round_up = vrmod10 >= 5; - vr = vrby10; vp = vpby10; vm = vmby10; - removed += 1; + if (0 <= i && i < m.nd: int) { + memio::appendrune(h, (m.buf[i] + '0'): rune)?; + } else { + memio::appendrune(h, '0')?; }; - output = vr + ibool(vr == vm || round_up); }; - const exp = e10 + removed; - return decf64 { exponent = exp, mantissa = output }; }; -type decf32 = struct { - mantissa: u32, - exponent: i32, +fn encode_e_mp( + m: *mp, + h: io::handle, + f: ffmt, + prec: (void | uint), + flag: fflags, +) (void | io::error) = { + assert(m.nd > 0); + memio::appendrune(h, (m.buf[0] + '0'): rune)?; + const zeros: uint = match (prec) { + case void => + yield 0; + case let u: uint => + yield switch (f) { + case ffmt::G => + yield if (m.nd + 1 < u) u - m.nd + 1 else 0; + case ffmt::E => + yield if (m.nd < u + 1) u - m.nd + 1 else 0; + case => abort(); + }; + }; + if (m.nd <= 1 && ffpoint(flag) && zeros < 1) { + zeros = 1; + }; + if (m.nd > 1 || zeros > 0) { + memio::appendrune(h, '.')?; + }; + for (let i = 1z; i < m.nd; i += 1) { + memio::appendrune(h, (m.buf[i] + '0'): rune)?; + }; + for (let i = 0u; i < zeros; i += 1) { + memio::appendrune(h, '0')?; + }; + memio::appendrune(h, if (ffcaps_exp(flag)) 'E' else 'e')?; + let e = m.dp - 1; + if (e < 0) { + e = -e; + memio::appendrune(h, '-')?; + } else if (ffpos_exp(flag)) { + memio::appendrune(h, '+')?; + }; + let ebuf: [3]u8 = [0...]; // max and min exponents are 3 digits + let l = declen(e: u64); + for (let i = 0z; i < l; i += 1) { + ebuf[2 - i] = (e % 10): u8; + e /= 10; + }; + if (fftwodigs(flag) && l == 1) { + l = 2; + }; + for (let i = 3 - l; i < 3; i += 1) { + memio::appendrune(h, (ebuf[i] + '0'): rune)?; + }; }; -fn f32todecf32(mantissa: u32, exponent: u32) decf32 = { - let e2 = (math::F32_EXPONENT_BIAS + math::F32_MANTISSA_BITS + 2): i32; - let m2: u32 = 0; - if (exponent == 0) { - e2 = 1 - e2; - m2 = mantissa; - } else { - e2 = (exponent: i32) - e2; - m2 = (1u32 << math::F32_MANTISSA_BITS: u32) | mantissa; +// Converts a [[types::floating]] to a string in base 10 and writes the result +// to the provided handle. Format parameters are as in [[ftosf]]. +export fn fftosf( + h: io::handle, + n: types::floating, + f: ffmt, + prec: (void | uint), + flag: fflags, +) (void | io::error) = { + const (mantissa, exponent, sign, special) = match (n) { + case let n: f64 => + const bits = math::f64bits(n); + const mantissa = bits & math::F64_MANTISSA_MASK; + const exponent = ((bits >> math::F64_MANTISSA_BITS) & + math::F64_EXPONENT_MASK): u32; + const sign = bits >> (math::F64_EXPONENT_BITS + + math::F64_MANTISSA_BITS) > 0; + const special = exponent == math::F64_EXPONENT_MASK; + yield (mantissa, exponent, sign, special); + case let n: f32 => + const bits = math::f32bits(n); + const mantissa = bits & math::F32_MANTISSA_MASK; + const exponent = ((bits >> math::F32_MANTISSA_BITS) & + math::F32_EXPONENT_MASK): u32; + const sign = bits >> (math::F32_EXPONENT_BITS + + math::F32_MANTISSA_BITS) > 0; + const special = exponent == math::F32_EXPONENT_MASK; + yield (mantissa, exponent, sign, special); }; - const accept_bounds = (m2 & 1) == 0; - const mv = 4 * m2, mp = mv + 2; - const mm_shift = ibool(mantissa != 0 || exponent <= 1); - const mm = mv - 1 - mm_shift; - let vr: u32 = 0, vp: u32 = 0, vm: u32 = 0; - let e10: i32 = 0; - let vm_trailing_zeroes = false, vr_trailing_zeroes = false; - let last_removed_digit: u8 = 0; - if (e2 >= 0) { - const q = log10pow2(e2: u32); - e10 = q: i32; - const k = F32_POW5_INV_BITCOUNT + pow5bits(q) - 1; - const i = -e2 + (q + k): i32; - vr = mulpow5inv_divpow2(mv, q, i); - vp = mulpow5inv_divpow2(mp, q, i); - vm = mulpow5inv_divpow2(mm, q, i); - if (q != 0 && (vp - 1) / 10 <= vm / 10) { - const l = F32_POW5_INV_BITCOUNT + pow5bits(q - 1) - 1; - last_removed_digit = (mulpow5inv_divpow2(mv, q - 1, - -e2 + ((q + l): i32) - 1) % 10): u8; - }; - if (q <= 9) { - if (mv % 5 == 0) { - vr_trailing_zeroes = pow5multiple32(mv, q); - } else if (accept_bounds) { - vm_trailing_zeroes = pow5multiple32(mm, q); - } else { - vp -= ibool(pow5multiple32(mp, q)); - }; - }; - } else { - const q = log10pow5((-e2): u32); - e10 = (q: i32) + e2; - const i = (-e2 - (q: i32)): u32; - const k = pow5bits(i) - F32_POW5_BITCOUNT; - let j = (q: i32) - k: i32; - vr = mulpow5_divpow2(mv, i, j); - vp = mulpow5_divpow2(mp, i, j); - vm = mulpow5_divpow2(mm, i, j); - if (q != 0 && (vp - 1) / 10 <= vm / 10) { - j = (q: i32) - 1 - (pow5bits(i + 1): i32 - - F32_POW5_BITCOUNT: i32); - last_removed_digit = (mulpow5_divpow2(mv, - (i + 1), j) % 10): u8; - }; - if (q <= 1) { - vr_trailing_zeroes = true; - if (accept_bounds) { - vm_trailing_zeroes = mm_shift == 1; - } else { - vp -= 1; - }; - } else if (q < 31) { - vr_trailing_zeroes = pow2multiple32(mv, q - 1); - }; + + if (special && mantissa != 0) { + return writestr(h, if (ffcaps(flag)) "NAN" else "nan"); }; - let removed: i32 = 0, output: u32 = 0; - if (vm_trailing_zeroes || vr_trailing_zeroes) { - for (vp / 10 > vm / 10) { - vm_trailing_zeroes &&= (vm - (vm / 10) * 10) == 0; - vr_trailing_zeroes &&= last_removed_digit == 0; - last_removed_digit = (vr % 10): u8; - vr /= 10; - vp /= 10; - vm /= 10; - removed += 1; - }; - if (vm_trailing_zeroes) { - for (vm % 10 == 0) { - vr_trailing_zeroes &&= last_removed_digit == 0; - last_removed_digit = (vr % 10): u8; - vr /= 10; - vp /= 10; - vm /= 10; - removed += 1; - }; - }; - if (vr_trailing_zeroes && last_removed_digit == 5 && vr % 2 == 0) { - // round to even - last_removed_digit = 4; + + if (sign) { + memio::appendrune(h, '-')?; + } else if (ffpos(flag)) { + memio::appendrune(h, '+')?; + }; + + if (special) { + return writestr(h, + if (ffcaps(flag)) "INFINITY" else "infinity"); + } else if (exponent == 0 && mantissa == 0) { + return encode_zero(h, f, prec, flag); + }; + + let m = mp { ... }; + let ok = false; + if (prec is void) { + // Shortest via Ryū. It is not correct to use f64todecf64 for + // f32s, they must be handled separately. + const (mdec, edec) = match (n) { + case f64 => + const d = f64todecf64(mantissa, exponent); + yield (d.mantissa, d.exponent); + case f32 => + const d = f32todecf32(mantissa: u32, exponent); + yield (d.mantissa: u64, d.exponent); }; - output = vr + ibool((vr == vm && - (!accept_bounds || !vm_trailing_zeroes)) || - last_removed_digit >= 5); - } else { - for (vp / 10 > vm / 10) { - last_removed_digit = (vr % 10): u8; - vr /= 10; - vp /= 10; - vm /= 10; - removed += 1; + init_mp_dec(&m, mdec, edec); + // If SHOW_POINT and we have too few digits, then we need to + // fall back to multiprecision. + ok = !ffpoint(flag) || m.dp < m.nd: int; + }; + + if (!ok) { + // Fall back to multiprecision. + match (n) { + case f64 => + init_mp(&m, mantissa, exponent, math::F64_EXPONENT_BIAS, + math::F64_MANTISSA_BITS); + case f32 => + init_mp(&m, mantissa, exponent, math::F32_EXPONENT_BIAS, + math::F32_MANTISSA_BITS); }; - output = vr + ibool(vr == vm || last_removed_digit >= 5); + trim_mp(&m); + const nd = compute_round_mp(&m, f, prec, flag); + round_mp(&m, nd); }; - const exp = e10 + removed; - return decf32 { mantissa = output, exponent = exp }; -}; -def F32_DECIMAL_DIGITS: i32 = 9; -def F64_DECIMAL_DIGITS: i32 = 17; + if (f == ffmt::G) { + trim_mp(&m); + }; -fn encode_base10(buf: []u8, mantissa: u64, exponent: i32, digits: i32) size = { - const n = mantissa, e = exponent, olen = declen(n); - const exp = e + olen: i32 - 1; - // use scientific notation for numbers whose exponent is beyond the - // precision available for the underlying type - if (exp > -4 && exp < digits) { - if (e >= 0) { - let k = exp; - buf[k + 1] = '.'; - buf[k + 2] = '0'; - for (let a = e; a > 0; a -= 1) { - buf[k] = '0'; - k -= 1; - }; - let m = n; - for (k >= 0; k -= 1) { - const mby10 = m / 10; - const mmod10 = (m: u32) - 10 * (mby10: u32); - buf[k] = '0' + mmod10: u8; - m = mby10; - }; - return (e + olen: i32 + 2): size; - } else if (exp < 0) { - buf[0] = '0'; - buf[1] = '.'; - let k = -e + 1; - let m = n; - for (let a = olen: i32; a > 0; a -= 1) { - const mby10 = m / 10; - const mmod10 = (m: u32) - 10 * (mby10: u32); - buf[k] = '0' + mmod10: u8; - k -= 1; - m = mby10; - }; - for (k > 1; k -= 1) { - buf[k] = '0'; - }; - return (-e + 2): size; - } else { - let k = olen: i32; - let m = n; - for (let a = -e; a > 0; a -= 1) { - const mby10 = m / 10; - const mmod10 = (m: u32) - 10 * (mby10: u32); - buf[k] = '0' + mmod10: u8; - k -= 1; - m = mby10; - }; - buf[k] = '.'; - k -= 1; - for (k >= 0; k -= 1) { - const mby10 = m / 10; - const mmod10 = (m: u32) - 10 * (mby10: u32); - buf[k] = '0' + mmod10: u8; - m = mby10; - }; - return (olen + 1): size; - }; + if (f == ffmt::G && prec is uint) { + if (prec as uint == 0) + prec = 1u; + }; + + if (m.nd == 0) { + // rounded to zero + encode_zero(h, f, prec, flag)?; + } else if (f == ffmt::E || (f == ffmt::G && + (m.dp < -1 || m.dp - m.nd: int > 2))) { + encode_e_mp(&m, h, f, prec, flag)?; } else { - let h: i32 = 0; - let m = n; - if (olen == 1) { - buf[0] = '0' + m: u8; - h = 1; - } else { - let k = olen: i32; - for (let a = k; a > 1; a -= 1) { - const mby10 = m / 10; - const mmod10 = (m: u32) - 10 * (mby10: u32); - buf[k] = '0' + mmod10: u8; - k -= 1; - m = mby10; - }; - buf[k] = '.'; - buf[0] = '0' + m: u8; - h = olen: i32 + 1; - }; - buf[h] = 'e'; - h += 1; - let ex = if (exp < 0) { - buf[h] = '-'; - h += 1; - yield -exp; - } else exp; - const elen = declen(ex: u64); - let k = h + elen: i32 - 1; - for (let a = elen: i32; a > 0; a -= 1) { - const eby10 = ex / 10; - const emod10 = (ex: u32) - 10 * (eby10: u32); - buf[k] = '0' + emod10: u8; - k -= 1; - ex = eby10; - }; - h += elen: i32; - return h: size; + encode_f_mp(&m, h, f, prec, flag)?; }; }; +// Converts any [[types::floating]] to a string in base 10. The return value +// must be freed. +// +// A precision of void yields the smallest number of digits that can be parsed +// into the exact same number. Otherwise, the meaning depends on f: +// - ffmt::F, ffmt::E: Number of digits after the decimal point. +// - ffmt::G: Number of significant digits. 0 is equivalent to 1 precision, and +// trailing zeros are removed. +export fn ftosf( + n: types::floating, + f: ffmt, + prec: (void | uint), + flag: fflags, +) str = { + let m = memio::dynamic(); + fftosf(&m, n, f, prec, flag)!; + return memio::string(&m)!; +}; + // Converts a f64 to a string in base 10. The return value is statically // allocated and will be overwritten on subsequent calls; see [[strings::dup]] -// to duplicate the result. +// to duplicate the result. The result is equivalent to [[ftosf]] with format G +// and precision void. export fn f64tos(n: f64) const str = { // The biggest string produced by a f64 number in base 10 would have the // negative sign, followed by a digit and decimal point, and then @@ -623,111 +337,23 @@ export fn f64tos(n: f64) const str = { // sign and the maximum of three digits for exponent. // (1 + 1 + 1 + 16 + 1 + 1 + 3) = 24 static let buf: [24]u8 = [0...]; - const bits = math::f64bits(n); - const sign = (bits >> (math::F64_MANTISSA_BITS + - math::F64_EXPONENT_BITS)): size; // bits >> 63 - const mantissa = bits & ((1u64 << math::F64_MANTISSA_BITS) - 1); - const exponent = ((bits >> math::F64_MANTISSA_BITS) & - (1u64 << math::F64_EXPONENT_BITS) - 1): u32; - if (mantissa == 0 && exponent == 0) { - return if (sign == 0) "0.0" else "-0.0"; - } else if (exponent == ((1 << math::F64_EXPONENT_BITS) - 1)) { - if (mantissa != 0) { - return "NaN"; - }; - return if (sign == 0) "Infinity" else "-Infinity"; - }; - const d = f64todecf64(mantissa, exponent); - if (sign != 0) { - buf[0] = '-'; - }; - let z = encode_base10(buf[sign..], d.mantissa, d.exponent, - F64_DECIMAL_DIGITS) + sign; - let s = types::string { data = &buf, length = z, ... }; - return *(&s: *str); + let m = memio::fixed(buf); + fftosf(&m, n, ffmt::G, void, 0)!; + return memio::string(&m)!; }; // Converts a f32 to a string in base 10. The return value is statically // allocated and will be overwritten on subsequent calls; see [[strings::dup]] -// to duplicate the result. +// to duplicate the result. The result is equivalent to [[ftosf]] with format G +// and precision void. export fn f32tos(n: f32) const str = { // The biggest string produced by a f32 number in base 10 would have the // negative sign, followed by a digit and decimal point, and then seven // more decimal digits, followed by 'e' and another negative sign and // the maximum of two digits for exponent. // (1 + 1 + 1 + 7 + 1 + 1 + 2) = 14 - static let buf: [16]u8 = [0...]; - const bits = math::f32bits(n); - const sign = (bits >> (math::F32_MANTISSA_BITS + - math::F32_EXPONENT_BITS)): size; // bits >> 31 - const mantissa = bits & ((1u32 << math::F32_MANTISSA_BITS) - 1): u32; - const exponent = (bits >> math::F32_MANTISSA_BITS): u32 & - ((1u32 << math::F32_EXPONENT_BITS) - 1): u32; - if (mantissa == 0 && exponent == 0) { - return if (sign == 0) "0.0" else "-0.0"; - } else if (exponent == ((1 << math::F32_EXPONENT_BITS) - 1): u32) { - if (mantissa != 0) { - return "NaN"; - }; - return if (sign == 0) "Infinity" else "-Infinity"; - }; - const d = f32todecf32(mantissa, exponent); - if (sign != 0) { - buf[0] = '-'; - }; - let z = encode_base10(buf[sign..], d.mantissa, d.exponent, - F32_DECIMAL_DIGITS) + sign; - const s = types::string { data = &buf, length = z, ... }; - return *(&s: *str); -}; - -@test fn f64tos() void = { - assert(f64tos(0.0) == "0.0"); - assert(f64tos(-0.0) == "-0.0"); - assert(f64tos(1.0 / 0.0) == "Infinity"); - assert(f64tos(-1.0 / 0.0) == "-Infinity"); - assert(f64tos(0.0 / 0.0) == "NaN"); - assert(f64tos(1.0) == "1.0"); - assert(f64tos(0.3) == "0.3"); - assert(f64tos(0.0031415) == "0.0031415"); - assert(f64tos(0.0000012345678) == "1.2345678e-6"); - assert(f64tos(1.414) == "1.414"); - assert(f64tos(1e234f64) == "1e234"); - assert(f64tos(1.2e-34) == "1.2e-34"); - assert(f64tos(-6345.9721) == "-6345.9721"); - assert(f64tos(1.23456789e67) == "1.23456789e67"); - assert(f64tos(11.2233445566778899e20) == "1.122334455667789e21"); - assert(f64tos(1000000.0e9) == "1000000000000000.0"); - assert(f64tos(9007199254740991.0) == "9007199254740991.0"); - assert(f64tos(90071992547409915.0) == "90071992547409920.0"); - assert(f64tos(90071992547409925.0) == "90071992547409920.0"); - assert(f64tos(5.0e-324) == "5e-324"); - assert(f64tos(2.2250738585072014e-308) == "2.2250738585072014e-308"); - assert(f64tos(1.7976931348623157e308) == "1.7976931348623157e308"); + static let buf: [14]u8 = [0...]; + let m = memio::fixed(buf); + fftosf(&m, n, ffmt::G, void, 0)!; + return memio::string(&m)!; }; - -@test fn f32tos() void = { - assert(f32tos(0.0) == "0.0"); - assert(f32tos(-0.0) == "-0.0"); - assert(f32tos(1.0 / 0.0) == "Infinity"); - assert(f32tos(-1.0 / 0.0) == "-Infinity"); - assert(f32tos(0.0 / 0.0) == "NaN"); - assert(f32tos(1.0) == "1.0"); - assert(f32tos(-8.0) == "-8.0"); - assert(f32tos(1.23) == "1.23"); - assert(f32tos(-0.618) == "-0.618"); - assert(f32tos(0.00456) == "0.00456"); - assert(f32tos(0.00000000000434655) == "4.34655e-12"); - assert(f32tos(123456.78) == "123456.78"); - assert(f32tos(-1.234567) == "-1.234567"); - assert(f32tos(12345.6789) == "12345.679"); - assert(f32tos(1.23e30) == "1.23e30"); - assert(f32tos(1.23e-30) == "1.23e-30"); - assert(f32tos(16777215.0) == "16777215.0"); - assert(f32tos(167772155.0) == "167772160.0"); - assert(f32tos(167772145.0) == "167772140.0"); - assert(f32tos(1.0e-45) == "1e-45"); - assert(f32tos(1.1754944e-38) == "1.1754944e-38"); - assert(f32tos(3.4028235e+38) == "3.4028235e38"); -}; - diff --git a/strconv/ftos_multiprecision.ha b/strconv/ftos_multiprecision.ha @@ -0,0 +1,305 @@ +// SPDX-License-Identifier: MPL-2.0 +// (c) Hare authors <https://harelang.org> + +// Multiprecision float to string based on golang's strconv/decimal.go. + +use ascii; +use math; +use strings; + +type mp = struct { + // Numbers 0-9, not ascii. Length for small numbers is + // log10(mantissa * 5^-exp). Subnormal doubles have min exp -1074 and + // max mantissa 4e16, giving at most 767 digits. + buf: [768]u8, + // Number of valid digits in buf. May be 0 if the number rounds to 0. + nd: uint, + // Decimal point index, may be negative. + // -1 means 0.0ddd... + // 0 means 0.ddd... + // 1 means d.dd... + // and so on + dp: int, +}; + +// These come from golang. The index into the table is amount of shift, up to +// 60. The number is the count of new digits that will be added by the shift, +// but one fewer if the number's prefix is smaller than the string prefix. +// +// For example, leftcheats[2] is (1, "25"). Any number left shifted by 2 will +// therefore be 1 digit longer, or zero digits longer if its first two digits +// are smaller than 25. +const leftcheats: [](size, str) = [ + (0, ""), + (1, "5"), + (1, "25"), + (1, "125"), + (2, "625"), + (2, "3125"), + (2, "15625"), + (3, "78125"), + (3, "390625"), + (3, "1953125"), + (4, "9765625"), + (4, "48828125"), + (4, "244140625"), + (4, "1220703125"), + (5, "6103515625"), + (5, "30517578125"), + (5, "152587890625"), + (6, "762939453125"), + (6, "3814697265625"), + (6, "19073486328125"), + (7, "95367431640625"), + (7, "476837158203125"), + (7, "2384185791015625"), + (7, "11920928955078125"), + (8, "59604644775390625"), + (8, "298023223876953125"), + (8, "1490116119384765625"), + (9, "7450580596923828125"), + (9, "37252902984619140625"), + (9, "186264514923095703125"), + (10, "931322574615478515625"), + (10, "4656612873077392578125"), + (10, "23283064365386962890625"), + (10, "116415321826934814453125"), + (11, "582076609134674072265625"), + (11, "2910383045673370361328125"), + (11, "14551915228366851806640625"), + (12, "72759576141834259033203125"), + (12, "363797880709171295166015625"), + (12, "1818989403545856475830078125"), + (13, "9094947017729282379150390625"), + (13, "45474735088646411895751953125"), + (13, "227373675443232059478759765625"), + (13, "1136868377216160297393798828125"), + (14, "5684341886080801486968994140625"), + (14, "28421709430404007434844970703125"), + (14, "142108547152020037174224853515625"), + (15, "710542735760100185871124267578125"), + (15, "3552713678800500929355621337890625"), + (15, "17763568394002504646778106689453125"), + (16, "88817841970012523233890533447265625"), + (16, "444089209850062616169452667236328125"), + (16, "2220446049250313080847263336181640625"), + (16, "11102230246251565404236316680908203125"), + (17, "55511151231257827021181583404541015625"), + (17, "277555756156289135105907917022705078125"), + (17, "1387778780781445675529539585113525390625"), + (18, "6938893903907228377647697925567626953125"), + (18, "34694469519536141888238489627838134765625"), + (18, "173472347597680709441192448139190673828125"), + (19, "867361737988403547205962240695953369140625"), +]; + +fn prefix_less_than_mp(m: *mp, s: str) bool = { + const u = strings::toutf8(s); + for (let i = 0z; i < len(s); i += 1) { + if (i >= m.nd) + return true; + if (m.buf[i] + '0': u8 != u[i]) + return m.buf[i] + '0': u8 < u[i]; + }; + return false; +}; + +// Shift left by k. +fn shl_mp(m: *mp, k: u64) void = { + let delta = leftcheats[k].0; + if (prefix_less_than_mp(m, leftcheats[k].1)) + delta -= 1; + let r = (m.nd - 1): int; + let w = m.nd + delta; + let n = 0u64; + for (r >= 0; r -= 1) { + n += m.buf[r]: u64 << k; + const quo = n / 10; + const rem = n - 10 * quo; + w -= 1; + m.buf[w] = rem: u8; + n = quo; + }; + for (n > 0) { + const quo = n / 10; + const rem = n - 10 * quo; + w -= 1; + m.buf[w] = rem: u8; + n = quo; + }; + m.nd += delta: uint; + m.dp += delta: int; +}; + +// Shift right by k. +fn shr_mp(m: *mp, k: u64) void = { + let r = 0z; + let w = 0z; + let n = 0u64; + const mask = (1 << k) - 1; + + for (n >> k == 0; r += 1) { + if (r >= m.nd) { + for (n >> k == 0) { + n *= 10; + r += 1; + }; + break; + }; + n = 10 * n + m.buf[r]; + }; + m.dp -= r: int - 1; + + for (r < m.nd; r += 1) { + const c = m.buf[r]; + const dig = n >> k; + n &= mask; + m.buf[w] = dig: u8; + w += 1; + n = n * 10 + c; + }; + + for (n > 0; w += 1) { + const dig = n >> k; + n &= mask; + m.buf[w] = dig: u8; + n = n * 10; + }; + m.nd = w: uint; +}; + +// Shift right (k < 0) or left (k > 0). We can only shift up to 60 at a time +// without losing bits, so break up big shifts. +fn shift_mp(m: *mp, k: int) void = { + if (k < 0) { + let nk = (-k): uint; + for (nk > 60) { + shr_mp(m, 60); + nk -= 60; + }; + shr_mp(m, nk); + } else if (k > 0) { + for (k > 60) { + shl_mp(m, 60); + k -= 60; + }; + shl_mp(m, k: uint); + }; +}; + +fn init_mp(m: *mp, mantissa: u64, exponent: u32, eb: u64, mb: u64) void = { + let e2 = (eb + mb): i32; + let m2: u64 = 0; + if (exponent == 0) { + e2 = 1 - e2; + m2 = mantissa; + } else { + e2 = (exponent: i32) - e2; + m2 = (1u64 << mb) | mantissa; + }; + + m.nd = declen(m2); + m.dp = m.nd: int; + for (let i = 0z; i < m.nd; i += 1) { + m.buf[m.nd - i - 1] = (m2 % 10): u8; + m2 /= 10; + }; + shift_mp(m, e2); +}; + +fn init_mp_dec(m: *mp, mantissa: u64, exponent: i32) void = { + const dl = declen(mantissa); + for (let i = 0u; i < dl; i += 1) { + m.buf[dl - i - 1] = (mantissa % 10): u8; + mantissa /= 10; + }; + m.nd = dl; + m.dp = dl: i32 + exponent; +}; + +fn round_up_mp(m: *mp) void = { + for (let i = 1z; i <= m.nd; i += 1) { + if (m.buf[m.nd - i] < 9) { + m.buf[m.nd - i] += 1; + return; + } else { + m.buf[m.nd - i] = 0; + }; + }; + // All high + m.buf[0] = 1; + m.nd = 1; + m.dp += 1; +}; + +// Compute the number of figs to round to for the given arguments. +fn compute_round_mp(m: *mp, f: ffmt, prec: (void | uint), flag: fflags) uint = { + // nd is the number of sig figs that we want to end up with + let nd: int = match (prec) { + case void => + // we should only get here if Ryu did not extend past the + // decimal point + assert(ffpoint(flag)); + yield m.nd: int + (if (m.dp > 0) m.dp else 0); + case let u: uint => + yield switch (f) { + case ffmt::E => + yield u: int + 1; + case ffmt::F => + yield u: int + m.dp; + case ffmt::G => + yield if (u == 0) 1 else u: int; + }; + }; + const nde = if (nd < 2) 2 else nd; + const ndf = if (m.dp >= 0 && nd < m.dp + 1) m.dp + 1 else nd; + if (ffpoint(flag)) { + nd = switch (f) { + case ffmt::E => + // need at least two digits, d.de0. + yield nde; + case ffmt::F => + // need enough to clear the decimal point by one. + yield ndf; + case ffmt::G => + // XXX: dup'd with the condition in ftosf_handle + if (m.dp < -1 || m.dp - m.nd: int > 2) + yield nde; + yield ndf; + }; + }; + if (nd <= 0) { + nd = 0; + }; + return if (nd: uint > m.nd) m.nd else nd: uint; +}; + +fn round_mp(m: *mp, nd: uint) void = { + assert(nd <= m.nd); + if (nd == m.nd) + return; + const oldnd = m.nd; + m.nd = nd; + if (m.buf[nd] > 5) { + round_up_mp(m); + } else if (m.buf[nd] == 5) { + let gt = false; + for (let i = m.nd + 1; i < oldnd; i += 1) { + if (m.buf[i] > 0) { + round_up_mp(m); + gt = true; + break; + }; + }; + if (!gt && nd > 0 && m.buf[nd - 1] & 1 > 0) { + round_up_mp(m); + }; + }; +}; + +// Remove trailing zeros. +fn trim_mp(m: *mp) void = { + for (m.nd > 1 && m.buf[m.nd - 1] == 0) { + m.nd -= 1; + }; +}; diff --git a/strconv/ftos_ryu.ha b/strconv/ftos_ryu.ha @@ -0,0 +1,503 @@ +// SPDX-License-Identifier: MPL-2.0 +// (c) Hare authors <https://harelang.org> + +// Using Ryū: fast float-to-string conversion algorithm by Ulf Adams. +// https://doi.org/10.1145/3192366.3192369 +// This Hare implementation is translated from the original +// C implementation here: https://github.com/ulfjack/ryu + +use ascii; +use math; +use types; + +type r128 = struct { + hi: u64, + lo: u64, +}; + +// TODO: use 128-bit integers when implemented +fn u128mul(a: u64, b: u64) r128 = { + const a0 = a: u32: u64, a1 = a >> 32; + const b0 = b: u32: u64, b1 = b >> 32; + const p00 = a0 * b0, p01 = a0 * b1, p10 = a1 * b0, p11 = a1 * b1; + const p00_lo = p00: u32: u64, p00_hi = p00 >> 32; + const mid1 = p10 + p00_hi; + const mid1_lo = mid1: u32: u64, mid1_hi = mid1 >> 32; + const mid2 = p01 + mid1_lo; + const mid2_lo = mid2: u32: u64, mid2_hi = mid2 >> 32; + const r_hi = p11 + mid1_hi + mid2_hi; + const r_lo = (mid2_lo << 32) | p00_lo; + return r128 { hi = r_hi, lo = r_lo }; +}; + +// TODO: Same as above +fn u128rshift(lo: u64, hi: u64, s: u32) u64 = { + assert(0 <= s); + assert(s <= 64); + return (hi << (64 - s)) | (lo >> s); +}; + +fn pow5fac(value: u64) u32 = { + const m_inv_5: u64 = 14757395258967641293; // 5 * m_inv_5 = 1 (mod 2^64) + const n_div_5: u64 = 3689348814741910323; + let count: u32 = 0; + for (true) { + assert(value != 0); + value *= m_inv_5; + if (value > n_div_5) break; + count += 1; + }; + return count; +}; + +fn pow5fac32(value: u32) u32 = { + let count: u32 = 0; + for (true) { + assert(value != 0); + const q = value / 5, r = value % 5; + if (r != 0) break; + value = q; + count += 1; + }; + return count; +}; + +fn ibool(b: bool) u8 = if (b) 1 else 0; + +fn pow5multiple(v: u64, p: u32) bool = pow5fac(v) >= p; +fn pow5multiple32(v: u32, p: u32) bool = pow5fac32(v) >= p; + +fn pow2multiple(v: u64, p: u32) bool = { + assert(v > 0); + assert(p < 64); + return (v & ((1u64 << p) - 1)) == 0; +}; + +fn pow2multiple32(v: u32, p: u32) bool = { + assert(v > 0); + assert(p < 32); + return (v & ((1u32 << p) - 1)) == 0; +}; + +fn mulshift64(m: u64, mul: (u64, u64), j: u32) u64 = { + // m is maximum 55 bits + let r0 = u128mul(m, mul.0), r1 = u128mul(m, mul.1); + const sum = r1.lo + r0.hi; + r1.hi += ibool(sum < r0.hi); + return u128rshift(sum, r1.hi, j - 64); +}; + +fn mulshiftall64( + m: u64, + mul: (u64, u64), + j: i32, + mm_shift: u32, +) (u64, u64, u64) = { + m <<= 1; + const r0 = u128mul(m, mul.0), r1 = u128mul(m, mul.1); + const lo = r0.lo, tmp = r0.hi, mid = tmp + r1.lo; + const hi = r1.hi + ibool(mid < tmp); + const lo2 = lo + mul.0; + const mid2 = mid + mul.1 + ibool(lo2 < lo); + const hi2 = hi + ibool(mid2 < mid); + const v_plus = u128rshift(mid2, hi2, (j - 64 - 1): u32); + const v_minus = if (mm_shift == 1) { + const lo3 = lo - mul.0; + const mid3 = mid - mul.1 - ibool(lo3 > lo); + const hi3 = hi - ibool(mid3 > mid); + yield u128rshift(mid3, hi3, (j - 64 - 1): u32); + } else { + const lo3 = lo + lo; + const mid3 = mid + mid + ibool(lo3 < lo); + const hi3 = hi + hi + ibool(mid3 < mid); + const lo4 = lo3 - mul.0; + const mid4 = mid3 - mul.1 - ibool(lo4 > lo3); + const hi4 = hi3 - ibool(mid4 > mid3); + yield u128rshift(mid4, hi4, (j - 64): u32); + }; + const v_rounded = u128rshift(mid, hi, (j - 64 - 1): u32); + return (v_plus, v_rounded, v_minus); +}; + +fn mulshift32(m: u32, a: u64, s: u32) u32 = { + assert(s > 32); + const a_lo = a: u32: u64, a_hi = a >> 32; + const b0 = m * a_lo, b1 = m * a_hi; + const sum = (b0 >> 32) + b1, ss = sum >> (s - 32); + assert(ss <= types::U32_MAX); + return ss: u32; +}; + +fn mulpow5inv_divpow2(m: u32, q: u32, j: i32) u32 = { + const pow5 = f64computeinvpow5(q); + return mulshift32(m, pow5.1 + 1, j: u32); +}; + +fn mulpow5_divpow2(m: u32, i: u32, j: i32) u32 = { + const pow5 = f64computepow5(i); + return mulshift32(m, pow5.1, j: u32); +}; + +fn log2pow5(e: u32) u32 = { + assert(e <= 3528); + return ((e * 1217359) >> 19); +}; + +fn ceil_log2pow5(e: u32) u32 = log2pow5(e) + 1; + +fn pow5bits(e: u32) u32 = ceil_log2pow5(e); + +fn log10pow2(e: u32) u32 = { + assert(e <= 1650); + return ((e * 78913) >> 18); +}; + +fn log10pow5(e: u32) u32 = { + assert(e <= 2620); + return ((e * 732923) >> 20); +}; + +def F64_POW5_INV_BITCOUNT: u8 = 125; +def F64_POW5_BITCOUNT: u8 = 125; + +def F32_POW5_INV_BITCOUNT: u8 = F64_POW5_INV_BITCOUNT - 64; +def F32_POW5_BITCOUNT: u8 = F64_POW5_BITCOUNT - 64; + +const F64_POW5_INV_SPLIT2: [15][2]u64 = [ + [1, 2305843009213693952], + [5955668970331000884, 1784059615882449851], + [8982663654677661702, 1380349269358112757], + [7286864317269821294, 2135987035920910082], + [7005857020398200553, 1652639921975621497], + [17965325103354776697, 1278668206209430417], + [8928596168509315048, 1978643211784836272], + [10075671573058298858, 1530901034580419511], + [597001226353042382, 1184477304306571148], + [1527430471115325346, 1832889850782397517], + [12533209867169019542, 1418129833677084982], + [5577825024675947042, 2194449627517475473], + [11006974540203867551, 1697873161311732311], + [10313493231639821582, 1313665730009899186], + [12701016819766672773, 2032799256770390445], +]; + +const POW5_INV_OFFSETS: [19]u32 = [ + 0x54544554, 0x04055545, 0x10041000, 0x00400414, 0x40010000, 0x41155555, + 0x00000454, 0x00010044, 0x40000000, 0x44000041, 0x50454450, 0x55550054, + 0x51655554, 0x40004000, 0x01000001, 0x00010500, 0x51515411, 0x05555554, + 0x00000000 +]; + +const F64_POW5_SPLIT2: [13][2]u64 = [ + [0, 1152921504606846976], + [0, 1490116119384765625], + [1032610780636961552, 1925929944387235853], + [7910200175544436838, 1244603055572228341], + [16941905809032713930, 1608611746708759036], + [13024893955298202172, 2079081953128979843], + [6607496772837067824, 1343575221513417750], + [17332926989895652603, 1736530273035216783], + [13037379183483547984, 2244412773384604712], + [1605989338741628675, 1450417759929778918], + [9630225068416591280, 1874621017369538693], + [665883850346957067, 1211445438634777304], + [14931890668723713708, 1565756531257009982] +]; + +const POW5_OFFSETS: [21]u32 = [ + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x40000000, 0x59695995, + 0x55545555, 0x56555515, 0x41150504, 0x40555410, 0x44555145, 0x44504540, + 0x45555550, 0x40004000, 0x96440440, 0x55565565, 0x54454045, 0x40154151, + 0x55559155, 0x51405555, 0x00000105 +]; + +def POW5_TABLE_SZ: u8 = 26; + +const POW5_TABLE: [POW5_TABLE_SZ]u64 = [ + 1u64, 5u64, 25u64, 125u64, 625u64, 3125u64, 15625u64, 78125u64, + 390625u64, 1953125u64, 9765625u64, 48828125u64, 244140625u64, + 1220703125u64, 6103515625u64, 30517578125u64, 152587890625u64, + 762939453125u64, 3814697265625u64, 19073486328125u64, 95367431640625u64, + 476837158203125u64, 2384185791015625u64, 11920928955078125u64, + 59604644775390625u64, 298023223876953125u64 //, 1490116119384765625u64 +]; + +fn f64computeinvpow5(i: u32) (u64, u64) = { + const base = ((i + POW5_TABLE_SZ - 1) / POW5_TABLE_SZ): u32; + const base2 = base * POW5_TABLE_SZ; + const mul = F64_POW5_INV_SPLIT2[base]; + const off = base2 - i; + if (off == 0) { + return (mul[0], mul[1]); + }; + const m = POW5_TABLE[off]; + const r1 = u128mul(m, mul[1]), r0 = u128mul(m, mul[0] - 1); + let high1 = r1.hi, low1 = r1.lo, high0 = r0.hi, low0 = r0.lo; + const sum = high0 + low1; + if (sum < high0) { + high1 += 1; + }; + const delta = pow5bits(base2) - pow5bits(i); + const res0 = u128rshift(low0, sum, delta) + 1 + + ((POW5_INV_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3); + const res1 = u128rshift(sum, high1, delta); + return (res0, res1); +}; + +fn f64computepow5(i: u32) (u64, u64) = { + const base = i / POW5_TABLE_SZ, base2 = base * POW5_TABLE_SZ; + const mul = F64_POW5_SPLIT2[base]; + const off = i - base2; + if (off == 0) { + return (mul[0], mul[1]); + }; + const m = POW5_TABLE[off]; + const r1 = u128mul(m, mul[1]), r0 = u128mul(m, mul[0]); + let high1 = r1.hi, low1 = r1.lo, high0 = r0.hi, low0 = r0.lo; + const sum = high0 + low1; + if (sum < high0) { + high1 += 1; + }; + const delta = pow5bits(i) - pow5bits(base2); + const res0 = u128rshift(low0, sum, delta) + + ((POW5_OFFSETS[i / 16] >> ((i % 16) << 1)) & 3); + const res1 = u128rshift(sum, high1, delta); + return (res0, res1); +}; + +type decf64 = struct { + mantissa: u64, + exponent: i32, +}; + +fn f64todecf64(mantissa: u64, exponent: u32) decf64 = { + let e2 = (math::F64_EXPONENT_BIAS + math::F64_MANTISSA_BITS + 2): i32; + let m2: u64 = 0; + if (exponent == 0) { + e2 = 1 - e2; + m2 = mantissa; + } else { + e2 = (exponent: i32) - e2; + m2 = (1u64 << math::F64_MANTISSA_BITS) | mantissa; + }; + const accept_bounds = (m2 & 1) == 0; + const mv = 4 * m2; + const mm_shift = ibool(mantissa != 0 || exponent <= 1); + let vp: u64 = 0, vr: u64 = 0, vm: u64 = 0; + let e10: i32 = 0; + let vm_trailing_zeros = false, vr_trailing_zeros = false; + if (e2 >= 0) { + const q = log10pow2(e2: u32) - ibool(e2 > 3); + e10 = q: i32; + const k = F64_POW5_INV_BITCOUNT + pow5bits(q) - 1; + const i = -e2 + (q + k): i32; + let pow5 = f64computeinvpow5(q); + const res = mulshiftall64(m2, pow5, i, mm_shift); + vp = res.0; vr = res.1; vm = res.2; + if (q <= 21) { + if ((mv - 5 * (mv / 5)) == 0) { + vr_trailing_zeros = pow5multiple(mv, q); + } else if (accept_bounds) { + vm_trailing_zeros = pow5multiple(mv - 1 - + mm_shift, q); + } else { + vp -= ibool(pow5multiple(mv + 2, q)); + }; + }; + } else { + const q = log10pow5((-e2): u32) - ibool(-e2 > 1); + e10 = e2 + (q: i32); + const i = -e2 - (q: i32); + const k = pow5bits(i: u32): i32 - F64_POW5_BITCOUNT: i32; + const j = (q: i32) - k; + let pow5 = f64computepow5(i: u32); + const res = mulshiftall64(m2, pow5, j, mm_shift); + vp = res.0; vr = res.1; vm = res.2; + if (q <= 1) { + vr_trailing_zeros = true; + if (accept_bounds) { + vm_trailing_zeros = mm_shift == 1; + } else { + vp -= 1; + }; + } else if (q < 63) { + vr_trailing_zeros = pow2multiple(mv, q); + }; + }; + let removed: i32 = 0, last_removed_digit: u8 = 0; + let output: u64 = 0; + if (vm_trailing_zeros || vr_trailing_zeros) { + for (true) { + const vpby10 = vp / 10, vmby10 = vm / 10; + if (vpby10 <= vmby10) break; + const vmmod10 = (vm: u32) - 10 * (vmby10: u32); + const vrby10 = vr / 10; + const vrmod10 = (vr: u32) - 10 * (vrby10: u32); + vm_trailing_zeros &&= vmmod10 == 0; + vr_trailing_zeros &&= last_removed_digit == 0; + last_removed_digit = vrmod10: u8; + vr = vrby10; vp = vpby10; vm = vmby10; + removed += 1; + }; + if (vm_trailing_zeros) { + for (true) { + const vmby10 = vm / 10; + const vmmod10 = (vm: u32) - 10 * (vmby10: u32); + if (vmmod10 != 0) break; + const vpby10 = vp / 10, vrby10 = vr / 10; + const vrmod10 = (vr: u32) - 10 * (vrby10: u32); + vr_trailing_zeros &&= last_removed_digit == 0; + last_removed_digit = vrmod10: u8; + vr = vrby10; vp = vpby10; vm = vmby10; + removed += 1; + }; + }; + if (vr_trailing_zeros && last_removed_digit == 5 && + (vr & 1 == 0)) { + // round to even + last_removed_digit = 4; + }; + output = vr + ibool((vr == vm && + (!accept_bounds || !vm_trailing_zeros)) || + last_removed_digit >= 5); + } else { + let round_up = false; + const vpby100 = vp / 100, vmby100 = vm / 100; + if (vpby100 > vmby100) { + const vrby100 = vr / 100; + const vrmod100 = (vr: u32) - 100 * (vrby100: u32); + round_up = vrmod100 >= 50; + vr = vrby100; vp = vpby100; vm = vmby100; + removed += 2; + }; + for (true) { + const vmby10 = vm / 10, vpby10 = vp / 10; + if (vpby10 <= vmby10) break; + const vrby10 = vr / 10; + const vrmod10 = (vr: u32) - 10 * (vrby10: u32); + round_up = vrmod10 >= 5; + vr = vrby10; vp = vpby10; vm = vmby10; + removed += 1; + }; + output = vr + ibool(vr == vm || round_up); + }; + const exp = e10 + removed; + return decf64 { exponent = exp, mantissa = output }; +}; + +type decf32 = struct { + mantissa: u32, + exponent: i32, +}; + +fn f32todecf32(mantissa: u32, exponent: u32) decf32 = { + let e2 = (math::F32_EXPONENT_BIAS + math::F32_MANTISSA_BITS + 2): i32; + let m2: u32 = 0; + if (exponent == 0) { + e2 = 1 - e2; + m2 = mantissa; + } else { + e2 = (exponent: i32) - e2; + m2 = (1u32 << math::F32_MANTISSA_BITS: u32) | mantissa; + }; + const accept_bounds = (m2 & 1) == 0; + const mv = 4 * m2, mp = mv + 2; + const mm_shift = ibool(mantissa != 0 || exponent <= 1); + const mm = mv - 1 - mm_shift; + let vr: u32 = 0, vp: u32 = 0, vm: u32 = 0; + let e10: i32 = 0; + let vm_trailing_zeroes = false, vr_trailing_zeroes = false; + let last_removed_digit: u8 = 0; + if (e2 >= 0) { + const q = log10pow2(e2: u32); + e10 = q: i32; + const k = F32_POW5_INV_BITCOUNT + pow5bits(q) - 1; + const i = -e2 + (q + k): i32; + vr = mulpow5inv_divpow2(mv, q, i); + vp = mulpow5inv_divpow2(mp, q, i); + vm = mulpow5inv_divpow2(mm, q, i); + if (q != 0 && (vp - 1) / 10 <= vm / 10) { + const l = F32_POW5_INV_BITCOUNT + pow5bits(q - 1) - 1; + last_removed_digit = (mulpow5inv_divpow2(mv, q - 1, + -e2 + ((q + l): i32) - 1) % 10): u8; + }; + if (q <= 9) { + if (mv % 5 == 0) { + vr_trailing_zeroes = pow5multiple32(mv, q); + } else if (accept_bounds) { + vm_trailing_zeroes = pow5multiple32(mm, q); + } else { + vp -= ibool(pow5multiple32(mp, q)); + }; + }; + } else { + const q = log10pow5((-e2): u32); + e10 = (q: i32) + e2; + const i = (-e2 - (q: i32)): u32; + const k = pow5bits(i) - F32_POW5_BITCOUNT; + let j = (q: i32) - k: i32; + vr = mulpow5_divpow2(mv, i, j); + vp = mulpow5_divpow2(mp, i, j); + vm = mulpow5_divpow2(mm, i, j); + if (q != 0 && (vp - 1) / 10 <= vm / 10) { + j = (q: i32) - 1 - (pow5bits(i + 1): i32 - + F32_POW5_BITCOUNT: i32); + last_removed_digit = (mulpow5_divpow2(mv, + (i + 1), j) % 10): u8; + }; + if (q <= 1) { + vr_trailing_zeroes = true; + if (accept_bounds) { + vm_trailing_zeroes = mm_shift == 1; + } else { + vp -= 1; + }; + } else if (q < 31) { + vr_trailing_zeroes = pow2multiple32(mv, q - 1); + }; + }; + let removed: i32 = 0, output: u32 = 0; + if (vm_trailing_zeroes || vr_trailing_zeroes) { + for (vp / 10 > vm / 10) { + vm_trailing_zeroes &&= (vm - (vm / 10) * 10) == 0; + vr_trailing_zeroes &&= last_removed_digit == 0; + last_removed_digit = (vr % 10): u8; + vr /= 10; + vp /= 10; + vm /= 10; + removed += 1; + }; + if (vm_trailing_zeroes) { + for (vm % 10 == 0) { + vr_trailing_zeroes &&= last_removed_digit == 0; + last_removed_digit = (vr % 10): u8; + vr /= 10; + vp /= 10; + vm /= 10; + removed += 1; + }; + }; + if (vr_trailing_zeroes && last_removed_digit == 5 && + vr % 2 == 0) { + // round to even + last_removed_digit = 4; + }; + output = vr + ibool((vr == vm && + (!accept_bounds || !vm_trailing_zeroes)) || + last_removed_digit >= 5); + } else { + for (vp / 10 > vm / 10) { + last_removed_digit = (vr % 10): u8; + vr /= 10; + vp /= 10; + vm /= 10; + removed += 1; + }; + output = vr + ibool(vr == vm || last_removed_digit >= 5); + }; + const exp = e10 + removed; + return decf32 { mantissa = output, exponent = exp }; +}; + +def F32_DECIMAL_DIGITS: i32 = 9; +def F64_DECIMAL_DIGITS: i32 = 17; diff --git a/strconv/ftostof+test.ha b/strconv/ftostof+test.ha @@ -1,33 +0,0 @@ -// SPDX-License-Identifier: MPL-2.0 -// (c) Hare authors <https://harelang.org> - -use math; - -// test round-trip accuracy of ftos followed by stof. -@test fn ftostof() void = { - const tcs: []f64 = [ - 1.0, - 0.0, - 0.1, - -0.0, - 1f64 / 3f64, - 4f64 / 3f64, - - 1f64 / 0f64, - -1f64 / 0f64, - - 0f64 / 0f64, - ]; - for (let i = 0z; i < len(tcs); i += 1) { - const res64 = stof64(f64tos(tcs[i]))!; - const res32 = stof32(f32tos(tcs[i]: f32))!; - // there can be multiple NaNs. only test isnan. - if (math::isnan(tcs[i])) { - assert(math::isnan(res64)); - assert(math::isnan(res32)); - } else { - assert(tcs[i] == res64); - assert(tcs[i]: f32 == res32); - }; - }; -};