Source File
	ftoadbox.go

Belonging Package
	internal/strconv

// 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
}

The pages are generated with Golds v0.8.3-preview. (GOOS=linux GOARCH=amd64)
Golds is a Go 101 project developed by Tapir Liu.
PR and bug reports are welcome and can be submitted to the issue list.
Please follow @zigo_101 (reachable from the left QR code) to get the latest news of Golds.