strconv: replace Ryu ftoa with Dragonbox

Dragonbox is a faster ftoa algorithm that provides the same guarantees
as Ryu: round-trip conversion, shortest length, and correct rounding.
Dragonbox only supports shortest-precision conversion, so we continue to
use Ryu-printf for fixed precision.

The new implementation has been fuzz-tested against the current
Ryu implementation in addition to the existing test suite.
Benchmarks show at least ~15-20% performance improvement.

The following shows the relevant output from benchstat. Full benchmark
results and plots are available at:
https://github.com/taichimaeda/dragonbox-bench/

goos: darwin
goarch: arm64
pkg: strconv
cpu: Apple M1
                                    │   old.txt    │               new.txt                │
                                    │    sec/op    │    sec/op     vs base                │
FormatFloat/Decimal-8                 32.71n ± 14%   31.89n ± 12%        ~ (p=0.165 n=10)
FormatFloat/Float-8                   45.54n ±  1%   42.48n ±  0%   -6.70% (p=0.000 n=10)
FormatFloat/Exp-8                     50.06n ±  0%   32.27n ±  1%  -35.54% (p=0.000 n=10)
FormatFloat/NegExp-8                  47.15n ±  1%   31.33n ±  0%  -33.56% (p=0.000 n=10)
FormatFloat/LongExp-8                 46.15n ±  1%   43.66n ±  0%   -5.38% (p=0.000 n=10)
FormatFloat/Big-8                     50.02n ±  0%   39.36n ±  0%  -21.31% (p=0.000 n=10)
FormatFloat/BinaryExp-8               27.89n ±  0%   27.88n ±  1%        ~ (p=0.798 n=10)
FormatFloat/32Integer-8               31.41n ±  0%   23.00n ±  3%  -26.79% (p=0.000 n=10)
FormatFloat/32ExactFraction-8         44.93n ±  1%   29.91n ±  0%  -33.43% (p=0.000 n=10)
FormatFloat/32Point-8                 43.22n ±  1%   33.82n ±  0%  -21.74% (p=0.000 n=10)
FormatFloat/32Exp-8                   45.91n ±  0%   25.48n ±  0%  -44.50% (p=0.000 n=10)
FormatFloat/32NegExp-8                44.66n ±  0%   25.12n ±  0%  -43.76% (p=0.000 n=10)
FormatFloat/32Shortest-8              37.96n ±  0%   27.83n ±  1%  -26.68% (p=0.000 n=10)
FormatFloat/Slowpath64-8              47.74n ±  2%   45.85n ±  0%   -3.96% (p=0.000 n=10)
FormatFloat/SlowpathDenormal64-8      42.78n ±  1%   41.46n ±  0%   -3.07% (p=0.000 n=10)
FormatFloat/ShorterIntervalCase32-8                  25.49n ±  2%
FormatFloat/ShorterIntervalCase64-8                  27.72n ±  1%
geomean                               41.95n         31.89n        -22.11%

Fixes #74886

Co-authored-by: Junekey Jeon <j6jeon@ucsd.edu>

Change-Id: I923f7259c9cecd0896b2340a43d1041cc2ed7787
GitHub-Last-Rev: fd735db
GitHub-Pull-Request: #75195
Reviewed-on: https://go-review.googlesource.com/c/go/+/700075
Reviewed-by: Russ Cox <rsc@golang.org>
Reviewed-by: Alan Donovan <adonovan@google.com>
TryBot-Bypass: Russ Cox <rsc@golang.org>

Author: taichimaeda

src/internal/strconv/ftoa.go

@@ -86,6 +86,7 @@ func genericFtoa(dst []byte, val float64, fmt byte, prec, bitSize int) []byte {
 	neg := bits>>(flt.expbits+flt.mantbits) != 0
 	exp := int(bits>>flt.mantbits) & (1<<flt.expbits - 1)
 	mant := bits & (uint64(1)<<flt.mantbits - 1)
+	denorm := false
 
 	switch exp {
 	case 1<<flt.expbits - 1:
@@ -104,6 +105,7 @@ func genericFtoa(dst []byte, val float64, fmt byte, prec, bitSize int) []byte {
 	case 0:
 		// denormalized
 		exp++
+		denorm = true
 
 	default:
 		// add implicit top bit
@@ -130,10 +132,10 @@ func genericFtoa(dst []byte, val float64, fmt byte, prec, bitSize int) []byte {
 		return formatDigits(dst, shortest, neg, digs, prec, fmt)
 	}
 	if shortest {
-		// Use Ryu algorithm.
+		// Use the Dragonbox algorithm.
 		var buf [32]byte
 		digs.d = buf[:]
-		ryuFtoaShortest(&digs, mant, exp-int(flt.mantbits), flt)
+		dboxFtoa(&digs, mant, exp-int(flt.mantbits), denorm, bitSize)
 		// Precision for shortest representation mode.
 		switch fmt {
 		case 'e', 'E':

src/internal/strconv/ftoa_test.go

@@ -351,6 +351,10 @@ var ftoaBenches = []struct {
 	// 622666234635.321497e-320 ~= 622666234635.3215e-320
 	// making it hard to find the 3rd digit
 	{"SlowpathDenormal64", 622666234635.3213e-320, 'e', -1, 64},
+
+	// Trigger the shorter interval case (3.90625e-3 = 1/256).
+	{"ShorterIntervalCase32", 3.90625e-3, 'e', -1, 32},
+	{"ShorterIntervalCase64", 3.90625e-3, 'e', -1, 64},
 }
 
 func BenchmarkFormatFloat(b *testing.B) {

src/internal/strconv/ftoadbox.go

@@ -0,0 +1,349 @@
+// Copyright 2025 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package strconv
+
+// Binary to decimal conversion using the Dragonbox algorithm by Junekey Jeon.
+//
+// Fixed precision format is not supported by the Dragonbox algorithm
+// so we continue to use Ryū-printf for this purpose.
+// See https://github.com/jk-jeon/dragonbox/issues/38 for more details.
+//
+// For binary to decimal rounding, uses round to nearest, tie to even.
+// For decimal to binary rounding, assumes round to nearest, tie to even.
+//
+// The original paper by Junekey Jeon can be found at:
+// https://github.com/jk-jeon/dragonbox/blob/d5dc40ae6a3f1a4559cda816738df2d6255b4e24/other_files/Dragonbox.pdf
+//
+// The reference implementation in C++ by Junekey Jeon can be found at:
+// https://github.com/jk-jeon/dragonbox/blob/6c7c925b571d54486b9ffae8d9d18a822801cbda/subproject/simple/include/simple_dragonbox.h
+
+// dragonboxFtoa computes the decimal significand and exponent
+// from the binary significand and exponent using the Dragonbox algorithm
+// and formats the decimal floating point number in d.
+func dboxFtoa(d *decimalSlice, mant uint64, exp int, denorm bool, bitSize int) {
+	if bitSize == 32 {
+		dboxFtoa32(d, uint32(mant), exp, denorm)
+		return
+	}
+	dboxFtoa64(d, mant, exp, denorm)
+}
+
+func dboxFtoa64(d *decimalSlice, mant uint64, exp int, denorm bool) {
+	if mant == 1<<float64MantBits && !denorm {
+		// Algorithm 5.6 (page 24).
+		k0 := -mulLog10_2MinusLog10_4Over3(exp)
+		φ, β := dboxPow64(k0, exp)
+		xi, zi := dboxRange64(φ, β)
+		if exp != 2 && exp != 3 {
+			xi++
+		}
+		q := zi / 10
+		if xi <= q*10 {
+			q, zeros := trimZeros(q)
+			dboxDigits(d, q, -k0+1+zeros)
+			return
+		}
+		yru := dboxRoundUp64(φ, β)
+		if exp == -77 && yru%2 != 0 {
+			yru--
+		} else if yru < xi {
+			yru++
+		}
+		dboxDigits(d, yru, -k0)
+		return
+	}
+
+	// κ = 2 for float64 (section 5.1.3)
+	const (
+		κ     = 2
+		p10κ  = 100       // 10**κ
+		p10κ1 = p10κ * 10 // 10**(κ+1)
+	)
+
+	// Algorithm 5.2 (page 15).
+	k0 := -mulLog10_2(exp)
+	φ, β := dboxPow64(κ+k0, exp)
+	zi, exact := dboxMulPow64(uint64(mant*2+1)<<β, φ)
+	s, r := zi/p10κ1, uint32(zi%p10κ1)
+	δi := dboxDelta64(φ, β)
+
+	if r < δi {
+		if r != 0 || !exact || mant%2 == 0 {
+			s, zeros := trimZeros(s)
+			dboxDigits(d, s, -k0+1+zeros)
+			return
+		}
+		s--
+		r = p10κ * 10
+	} else if r == δi {
+		parity, exact := dboxParity64(uint64(mant*2-1), φ, β)
+		if parity || (exact && mant%2 == 0) {
+			s, zeros := trimZeros(s)
+			dboxDigits(d, s, -k0+1+zeros)
+			return
+		}
+	}
+
+	// Algorithm 5.4 (page 18).
+	D := r + p10κ/2 - δi/2
+	t, ρ := D/p10κ, D%p10κ
+	yru := 10*s + uint64(t)
+	if ρ == 0 {
+		parity, exact := dboxParity64(mant*2, φ, β)
+		if parity != ((D-p10κ/2)%2 != 0) || exact && yru%2 != 0 {
+			yru--
+		}
+	}
+	dboxDigits(d, yru, -k0)
+}
+
+// Almost identical to dragonboxFtoa64.
+// This is kept as a separate copy to minimize runtime overhead.
+func dboxFtoa32(d *decimalSlice, mant uint32, exp int, denorm bool) {
+	if mant == 1<<float32MantBits && !denorm {
+		// Algorithm 5.6 (page 24).
+		k0 := -mulLog10_2MinusLog10_4Over3(exp)
+		φ, β := dboxPow32(k0, exp)
+		xi, zi := dboxRange32(φ, β)
+		if exp != 2 && exp != 3 {
+			xi++
+		}
+		q := zi / 10
+		if xi <= q*10 {
+			q, zeros := trimZeros(uint64(q))
+			dboxDigits(d, q, -k0+1+zeros)
+			return
+		}
+		yru := dboxRoundUp32(φ, β)
+		if exp == -77 && yru%2 != 0 {
+			yru--
+		} else if yru < xi {
+			yru++
+		}
+		dboxDigits(d, uint64(yru), -k0)
+		return
+	}
+
+	// κ = 1 for float32 (section 5.1.3)
+	const (
+		κ     = 1
+		p10κ  = 10
+		p10κ1 = p10κ * 10
+	)
+
+	// Algorithm 5.2 (page 15).
+	k0 := -mulLog10_2(exp)
+	φ, β := dboxPow32(κ+k0, exp)
+	zi, exact := dboxMulPow32(uint32(mant*2+1)<<β, φ)
+	s, r := zi/p10κ1, uint32(zi%p10κ1)
+	δi := dboxDelta32(φ, β)
+
+	if r < δi {
+		if r != 0 || !exact || mant%2 == 0 {
+			s, zeros := trimZeros(uint64(s))
+			dboxDigits(d, s, -k0+1+zeros)
+			return
+		}
+		s--
+		r = p10κ * 10
+	} else if r == δi {
+		parity, exact := dboxParity32(uint32(mant*2-1), φ, β)
+		if parity || (exact && mant%2 == 0) {
+			s, zeros := trimZeros(uint64(s))
+			dboxDigits(d, s, -k0+1+zeros)
+			return
+		}
+	}
+
+	// Algorithm 5.4 (page 18).
+	D := r + p10κ/2 - δi/2
+	t, ρ := D/p10κ, D%p10κ
+	yru := 10*s + uint32(t)
+	if ρ == 0 {
+		parity, exact := dboxParity32(mant*2, φ, β)
+		if parity != ((D-p10κ/2)%2 != 0) || exact && yru%2 != 0 {
+			yru--
+		}
+	}
+	dboxDigits(d, uint64(yru), -k0)
+}
+
+// dboxDigits emits decimal digits of mant in d for float64
+// and adjusts the decimal point based on exp.
+func dboxDigits(d *decimalSlice, mant uint64, exp int) {
+	i := formatBase10(d.d, mant)
+	d.d = d.d[i:]
+	d.nd = len(d.d)
+	d.dp = d.nd + exp
+}
+
+// uadd128 returns the full 128 bits of u + n.
+func uadd128(u uint128, n uint64) uint128 {
+	sum := uint64(u.Lo + n)
+	// Check if lo is wrapped around.
+	if sum < u.Lo {
+		u.Hi++
+	}
+	u.Lo = sum
+	return u
+}
+
+// umul64 returns the full 64 bits of x * y.
+func umul64(x, y uint32) uint64 {
+	return uint64(x) * uint64(y)
+}
+
+// umul96Upper64 returns the upper 64 bits (out of 96 bits) of x * y.
+func umul96Upper64(x uint32, y uint64) uint64 {
+	yh := uint32(y >> 32)
+	yl := uint32(y)
+
+	xyh := umul64(x, yh)
+	xyl := umul64(x, yl)
+
+	return xyh + (xyl >> 32)
+}
+
+// umul96Lower64 returns the lower 64 bits (out of 96 bits) of x * y.
+func umul96Lower64(x uint32, y uint64) uint64 {
+	return uint64(uint64(x) * y)
+}
+
+// umul128Upper64 returns the upper 64 bits (out of 128 bits) of x * y.
+func umul128Upper64(x, y uint64) uint64 {
+	a := uint32(x >> 32)
+	b := uint32(x)
+	c := uint32(y >> 32)
+	d := uint32(y)
+
+	ac := umul64(a, c)
+	bc := umul64(b, c)
+	ad := umul64(a, d)
+	bd := umul64(b, d)
+
+	intermediate := (bd >> 32) + uint64(uint32(ad)) + uint64(uint32(bc))
+
+	return ac + (intermediate >> 32) + (ad >> 32) + (bc >> 32)
+}
+
+// umul192Upper128 returns the upper 128 bits (out of 192 bits) of x * y.
+func umul192Upper128(x uint64, y uint128) uint128 {
+	r := umul128(x, y.Hi)
+	t := umul128Upper64(x, y.Lo)
+	return uadd128(r, t)
+}
+
+// umul192Lower128 returns the lower 128 bits (out of 192 bits) of x * y.
+func umul192Lower128(x uint64, y uint128) uint128 {
+	high := x * y.Hi
+	highLow := umul128(x, y.Lo)
+	return uint128{uint64(high + highLow.Hi), highLow.Lo}
+}
+
+// dboxMulPow64 computes x^(i), y^(i), z^(i)
+// from the precomputed value of φ̃k for float64
+// and also checks if x^(f), y^(f), z^(f) == 0 (section 5.2.1).
+func dboxMulPow64(u uint64, phi uint128) (intPart uint64, isInt bool) {
+	r := umul192Upper128(u, phi)
+	intPart = r.Hi
+	isInt = r.Lo == 0
+	return
+}
+
+// dboxMulPow32 computes x^(i), y^(i), z^(i)
+// from the precomputed value of φ̃k for float32
+// and also checks if x^(f), y^(f), z^(f) == 0 (section 5.2.1).
+func dboxMulPow32(u uint32, phi uint64) (intPart uint32, isInt bool) {
+	r := umul96Upper64(u, phi)
+	intPart = uint32(r >> 32)
+	isInt = uint32(r) == 0
+	return
+}
+
+// dboxParity64 computes only the parity of x^(i), y^(i), z^(i)
+// from the precomputed value of φ̃k for float64
+// and also checks if x^(f), y^(f), z^(f) = 0 (section 5.2.1).
+func dboxParity64(mant2 uint64, phi uint128, beta int) (parity bool, isInt bool) {
+	r := umul192Lower128(mant2, phi)
+	parity = ((r.Hi >> (64 - beta)) & 1) != 0
+	isInt = ((uint64(r.Hi << beta)) | (r.Lo >> (64 - beta))) == 0
+	return
+}
+
+// dboxParity32 computes only the parity of x^(i), y^(i), z^(i)
+// from the precomputed value of φ̃k for float32
+// and also checks if x^(f), y^(f), z^(f) = 0 (section 5.2.1).
+func dboxParity32(mant2 uint32, phi uint64, beta int) (parity bool, isInt bool) {
+	r := umul96Lower64(mant2, phi)
+	parity = ((r >> (64 - beta)) & 1) != 0
+	isInt = uint32(r>>(32-beta)) == 0
+	return
+}
+
+// dboxDelta64 returns δ^(i) from the precomputed value of φ̃k for float64.
+func dboxDelta64(φ uint128, β int) uint32 {
+	return uint32(φ.Hi >> (64 - 1 - β))
+}
+
+// dboxDelta32 returns δ^(i) from the precomputed value of φ̃k for float32.
+func dboxDelta32(φ uint64, β int) uint32 {
+	return uint32(φ >> (64 - 1 - β))
+}
+
+// mulLog10_2MinusLog10_4Over3 computes
+// ⌊e*log10(2)-log10(4/3)⌋ = ⌊log10(2^e)-log10(4/3)⌋ (section 6.3).
+func mulLog10_2MinusLog10_4Over3(e int) int {
+	// e should be in the range [-2985, 2936].
+	return (e*631305 - 261663) >> 21
+}
+
+const (
+	floatMantBits64 = 52 // p = 52 for float64.
+	floatMantBits32 = 23 // p = 23 for float32.
+)
+
+// dboxRange64 returns the left and right float64 endpoints.
+func dboxRange64(φ uint128, β int) (left, right uint64) {
+	left = (φ.Hi - (φ.Hi >> (float64MantBits + 2))) >> (64 - float64MantBits - 1 - β)
+	right = (φ.Hi + (φ.Hi >> (float64MantBits + 1))) >> (64 - float64MantBits - 1 - β)
+	return left, right
+}
+
+// dboxRange32 returns the left and right float32 endpoints.
+func dboxRange32(φ uint64, β int) (left, right uint32) {
+	left = uint32((φ - (φ >> (floatMantBits32 + 2))) >> (64 - floatMantBits32 - 1 - β))
+	right = uint32((φ + (φ >> (floatMantBits32 + 1))) >> (64 - floatMantBits32 - 1 - β))
+	return left, right
+}
+
+// dboxRoundUp64 computes the round up of y (i.e., y^(ru)).
+func dboxRoundUp64(phi uint128, beta int) uint64 {
+	return (phi.Hi>>(128/2-floatMantBits64-2-beta) + 1) / 2
+}
+
+// dboxRoundUp32 computes the round up of y (i.e., y^(ru)).
+func dboxRoundUp32(phi uint64, beta int) uint32 {
+	return uint32(phi>>(64-floatMantBits32-2-beta)+1) / 2
+}
+
+// dboxPow64 gets the precomputed value of φ̃̃k for float64.
+func dboxPow64(k, e int) (φ uint128, β int) {
+	φ, e1, _ := pow10(k)
+	if k < 0 || k > 55 {
+		φ.Lo++
+	}
+	β = e + e1 - 1
+	return φ, β
+}
+
+// dboxPow32 gets the precomputed value of φ̃̃k for float32.
+func dboxPow32(k, e int) (mant uint64, exp int) {
+	m, e1, _ := pow10(k)
+	if k < 0 || k > 27 {
+		m.Hi++
+	}
+	exp = e + e1 - 1
+	return m.Hi, exp
+}