ftoafixed.go

     1  // Copyright 2025 The Go Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package strconv
     6  
     7  import "math/bits"
     8  
     9  var uint64pow10 = [...]uint64{
    10  	1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
    11  	1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
    12  }
    13  
    14  // fixedFtoa formats a number of decimal digits of mant*(2^exp) into d,
    15  // where mant > 0 and 1 ≤ digits ≤ 18.
    16  // If fmt == 'f', digits is a conservative overestimate, and the final
    17  // number of digits is prec past the decimal point.
    18  func fixedFtoa(d *decimalSlice, mant uint64, exp, digits, prec int, fmt byte) {
    19  	// The strategy here is to multiply (mant * 2^exp) by a power of 10
    20  	// to make the resulting integer be the number of digits we want.
    21  	//
    22  	// Adams proved in the Ryu paper that 128-bit precision in the
    23  	// power-of-10 constant is sufficient to produce correctly
    24  	// rounded output for all float64s, up to 18 digits.
    25  	// https://dl.acm.org/doi/10.1145/3192366.3192369
    26  	//
    27  	// TODO(rsc): The paper is not focused on, nor terribly clear about,
    28  	// this fact in this context, and the proof seems too complicated.
    29  	// Post a shorter, more direct proof and link to it here.
    30  
    31  	if digits > 18 {
    32  		panic("fixedFtoa called with digits > 18")
    33  	}
    34  
    35  	// Shift mantissa to have 64 bits,
    36  	// so that the 192-bit product below will
    37  	// have at least 63 bits in its top word.
    38  	b := 64 - bits.Len64(mant)
    39  	mant <<= b
    40  	exp -= b
    41  
    42  	// We have f = mant * 2^exp ≥ 2^(63+exp)
    43  	// and we want to multiply it by some 10^p
    44  	// to make it have the number of digits plus one rounding bit:
    45  	//
    46  	//	2 * 10^(digits-1) ≤ f * 10^p < ~2 * 10^digits
    47  	//
    48  	// The lower bound is required, but the upper bound is approximate:
    49  	// we must not have too few digits, but we can round away extra ones.
    50  	//
    51  	//	f * 10^p ≥ 2 * 10^(digits-1)
    52  	//	10^p ≥ 2 * 10^(digits-1) / f                         [dividing by f]
    53  	//	p ≥ (log₁₀ 2) + (digits-1) - log₁₀ f                 [taking log₁₀]
    54  	//	p ≥ (log₁₀ 2) + (digits-1) - log₁₀ (mant * 2^exp)    [expanding f]
    55  	//	p ≥ (log₁₀ 2) + (digits-1) - (log₁₀ 2) * (64 + exp)  [mant < 2⁶⁴]
    56  	//	p ≥ (digits - 1) - (log₁₀ 2) * (63 + exp)            [refactoring]
    57  	//
    58  	// Once we have p, we can compute the scaled value:
    59  	//
    60  	//	dm * 2^de = mant * 2^exp * 10^p
    61  	//	          = mant * 2^exp * pow/2^128 * 2^exp2.
    62  	//	          = (mant * pow/2^128) * 2^(exp+exp2).
    63  	p := (digits - 1) - mulLog10_2(63+exp)
    64  	pow, exp2, ok := pow10(p)
    65  	if !ok {
    66  		// This never happens due to the range of float32/float64 exponent
    67  		panic("fixedFtoa: pow10 out of range")
    68  	}
    69  	if -22 <= p && p < 0 {
    70  		// Special case: Let q=-p. q is in [1,22]. We are dividing by 10^q
    71  		// and the mantissa may be a multiple of 5^q (5^22 < 2^53),
    72  		// in which case the division must be computed exactly and
    73  		// recorded as exact for correct rounding. Our normal computation is:
    74  		//
    75  		//	dm = floor(mant * floor(10^p * 2^s))
    76  		//
    77  		// for some scaling shift s. To make this an exact division,
    78  		// it suffices to change the inner floor to a ceil:
    79  		//
    80  		//	dm = floor(mant * ceil(10^p * 2^s))
    81  		//
    82  		// In the range of values we are using, the floor and ceil
    83  		// cancel each other out and the high 64 bits of the product
    84  		// come out exactly right.
    85  		// (This is the same trick compilers use for division by constants.
    86  		// See Hacker's Delight, 2nd ed., Chapter 10.)
    87  		pow.Lo++
    88  	}
    89  	dm, lo1, lo0 := umul192(mant, pow)
    90  	de := exp + exp2
    91  
    92  	// Check whether any bits have been truncated from dm.
    93  	// If so, set dt != 0. If not, leave dt == 0 (meaning dm is exact).
    94  	var dt uint
    95  	switch {
    96  	default:
    97  		// Most powers of 10 use a truncated constant,
    98  		// meaning the result is also truncated.
    99  		dt = 1
   100  	case 0 <= p && p <= 55:
   101  		// Small positive powers of 10 (up to 10⁵⁵) can be represented
   102  		// precisely in a 128-bit mantissa (5⁵⁵ ≤ 2¹²⁸), so the only truncation
   103  		// comes from discarding the low bits of the 192-bit product.
   104  		//
   105  		// TODO(rsc): The new proof mentioned above should also
   106  		// prove that we can't have lo1 == 0 and lo0 != 0.
   107  		// After proving that, drop computation and use of lo0 here.
   108  		dt = bool2uint(lo1|lo0 != 0)
   109  	case -22 <= p && p < 0 && divisiblePow5(mant, -p):
   110  		// If the original mantissa was a multiple of 5^p,
   111  		// the result is exact. (See comment above for pow.Lo++.)
   112  		dt = 0
   113  	}
   114  
   115  	// The value we want to format is dm * 2^de, where de < 0.
   116  	// Multply by 2^de by shifting, but leave one extra bit for rounding.
   117  	// After the shift, the "integer part" of dm is dm>>1,
   118  	// the "rounding bit" (the first fractional bit) is dm&1,
   119  	// and the "truncated bit" (have any bits been discarded?) is dt.
   120  	shift := -de - 1
   121  	dt |= bool2uint(dm&(1<<shift-1) != 0)
   122  	dm >>= shift
   123  
   124  	// Set decimal point in eventual formatted digits,
   125  	// so we can update it as we adjust the digits.
   126  	d.dp = digits - p
   127  
   128  	// Trim excess digit if any, updating truncation and decimal point.
   129  	// The << 1 is leaving room for the rounding bit.
   130  	max := uint64pow10[digits] << 1
   131  	if dm >= max {
   132  		var r uint
   133  		dm, r = dm/10, uint(dm%10)
   134  		dt |= bool2uint(r != 0)
   135  		d.dp++
   136  	}
   137  
   138  	// If this is %.*f we may have overestimated the digits needed.
   139  	// Now that we know where the decimal point is,
   140  	// trim to the actual number of digits, which is d.dp+prec.
   141  	if fmt == 'f' && digits != d.dp+prec {
   142  		for digits > d.dp+prec {
   143  			var r uint
   144  			dm, r = dm/10, uint(dm%10)
   145  			dt |= bool2uint(r != 0)
   146  			digits--
   147  		}
   148  
   149  		// Dropping those digits can create a new leftmost
   150  		// non-zero digit, like if we are formatting %.1f and
   151  		// convert 0.09 -> 0.1. Detect and adjust for that.
   152  		if digits <= 0 {
   153  			digits = 1
   154  			d.dp++
   155  		}
   156  
   157  		max = uint64pow10[digits] << 1
   158  	}
   159  
   160  	// Round and shift away rounding bit.
   161  	// We want to round up when
   162  	// (a) the fractional part is > 0.5 (dm&1 != 0 and dt == 1)
   163  	// (b) or the fractional part is ≥ 0.5 and the integer part is odd
   164  	//     (dm&1 != 0 and dm&2 != 0).
   165  	// The bitwise expression encodes that logic.
   166  	dm += uint64(uint(dm) & (dt | uint(dm)>>1) & 1)
   167  	dm >>= 1
   168  	if dm == max>>1 {
   169  		// 999... rolled over to 1000...
   170  		dm = uint64pow10[digits-1]
   171  		d.dp++
   172  	}
   173  
   174  	// Format digits into d.
   175  	if dm != 0 {
   176  		if formatBase10(d.d[:digits], dm) != 0 {
   177  			panic("formatBase10")
   178  		}
   179  		d.nd = digits
   180  		for d.d[d.nd-1] == '0' {
   181  			d.nd--
   182  		}
   183  	}
   184  }
   185
View as plain text