# Source file src/strconv/ftoaryu.go

```     1  // Copyright 2021 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
4
5  package strconv
6
7  import (
8  	"math/bits"
9  )
10
11  // binary to decimal conversion using the Ryū algorithm.
12  //
13  // See Ulf Adams, "Ryū: Fast Float-to-String Conversion" (doi:10.1145/3192366.3192369)
14  //
15  // Fixed precision formatting is a variant of the original paper's
16  // algorithm, where a single multiplication by 10^k is required,
17  // sharing the same rounding guarantees.
18
19  // ryuFtoaFixed32 formats mant*(2^exp) with prec decimal digits.
20  func ryuFtoaFixed32(d *decimalSlice, mant uint32, exp int, prec int) {
21  	if prec < 0 {
22  		panic("ryuFtoaFixed32 called with negative prec")
23  	}
24  	if prec > 9 {
25  		panic("ryuFtoaFixed32 called with prec > 9")
26  	}
27  	// Zero input.
28  	if mant == 0 {
29  		d.nd, d.dp = 0, 0
30  		return
31  	}
32  	// Renormalize to a 25-bit mantissa.
33  	e2 := exp
34  	if b := bits.Len32(mant); b < 25 {
35  		mant <<= uint(25 - b)
36  		e2 += b - 25
37  	}
38  	// Choose an exponent such that rounded mant*(2^e2)*(10^q) has
39  	// at least prec decimal digits, i.e
40  	//     mant*(2^e2)*(10^q) >= 10^(prec-1)
41  	// Because mant >= 2^24, it is enough to choose:
42  	//     2^(e2+24) >= 10^(-q+prec-1)
43  	// or q = -mulByLog2Log10(e2+24) + prec - 1
44  	q := -mulByLog2Log10(e2+24) + prec - 1
45
46  	// Now compute mant*(2^e2)*(10^q).
47  	// Is it an exact computation?
48  	// Only small positive powers of 10 are exact (5^28 has 66 bits).
49  	exact := q <= 27 && q >= 0
50
51  	di, dexp2, d0 := mult64bitPow10(mant, e2, q)
52  	if dexp2 >= 0 {
53  		panic("not enough significant bits after mult64bitPow10")
54  	}
55  	// As a special case, computation might still be exact, if exponent
56  	// was negative and if it amounts to computing an exact division.
57  	// In that case, we ignore all lower bits.
58  	// Note that division by 10^11 cannot be exact as 5^11 has 26 bits.
59  	if q < 0 && q >= -10 && divisibleByPower5(uint64(mant), -q) {
60  		exact = true
61  		d0 = true
62  	}
63  	// Remove extra lower bits and keep rounding info.
64  	extra := uint(-dexp2)
65  	extraMask := uint32(1<<extra - 1)
66
67  	di, dfrac := di>>extra, di&extraMask
68  	roundUp := false
69  	if exact {
70  		// If we computed an exact product, d + 1/2
71  		// should round to d+1 if 'd' is odd.
72  		roundUp = dfrac > 1<<(extra-1) ||
73  			(dfrac == 1<<(extra-1) && !d0) ||
74  			(dfrac == 1<<(extra-1) && d0 && di&1 == 1)
75  	} else {
76  		// otherwise, d+1/2 always rounds up because
77  		// we truncated below.
78  		roundUp = dfrac>>(extra-1) == 1
79  	}
80  	if dfrac != 0 {
81  		d0 = false
82  	}
83  	// Proceed to the requested number of digits
84  	formatDecimal(d, uint64(di), !d0, roundUp, prec)
86  	d.dp -= q
87  }
88
89  // ryuFtoaFixed64 formats mant*(2^exp) with prec decimal digits.
90  func ryuFtoaFixed64(d *decimalSlice, mant uint64, exp int, prec int) {
91  	if prec > 18 {
92  		panic("ryuFtoaFixed64 called with prec > 18")
93  	}
94  	// Zero input.
95  	if mant == 0 {
96  		d.nd, d.dp = 0, 0
97  		return
98  	}
99  	// Renormalize to a 55-bit mantissa.
100  	e2 := exp
101  	if b := bits.Len64(mant); b < 55 {
102  		mant = mant << uint(55-b)
103  		e2 += b - 55
104  	}
105  	// Choose an exponent such that rounded mant*(2^e2)*(10^q) has
106  	// at least prec decimal digits, i.e
107  	//     mant*(2^e2)*(10^q) >= 10^(prec-1)
108  	// Because mant >= 2^54, it is enough to choose:
109  	//     2^(e2+54) >= 10^(-q+prec-1)
110  	// or q = -mulByLog2Log10(e2+54) + prec - 1
111  	//
112  	// The minimal required exponent is -mulByLog2Log10(1025)+18 = -291
113  	// The maximal required exponent is mulByLog2Log10(1074)+18 = 342
114  	q := -mulByLog2Log10(e2+54) + prec - 1
115
116  	// Now compute mant*(2^e2)*(10^q).
117  	// Is it an exact computation?
118  	// Only small positive powers of 10 are exact (5^55 has 128 bits).
119  	exact := q <= 55 && q >= 0
120
121  	di, dexp2, d0 := mult128bitPow10(mant, e2, q)
122  	if dexp2 >= 0 {
123  		panic("not enough significant bits after mult128bitPow10")
124  	}
125  	// As a special case, computation might still be exact, if exponent
126  	// was negative and if it amounts to computing an exact division.
127  	// In that case, we ignore all lower bits.
128  	// Note that division by 10^23 cannot be exact as 5^23 has 54 bits.
129  	if q < 0 && q >= -22 && divisibleByPower5(mant, -q) {
130  		exact = true
131  		d0 = true
132  	}
133  	// Remove extra lower bits and keep rounding info.
134  	extra := uint(-dexp2)
135  	extraMask := uint64(1<<extra - 1)
136
137  	di, dfrac := di>>extra, di&extraMask
138  	roundUp := false
139  	if exact {
140  		// If we computed an exact product, d + 1/2
141  		// should round to d+1 if 'd' is odd.
142  		roundUp = dfrac > 1<<(extra-1) ||
143  			(dfrac == 1<<(extra-1) && !d0) ||
144  			(dfrac == 1<<(extra-1) && d0 && di&1 == 1)
145  	} else {
146  		// otherwise, d+1/2 always rounds up because
147  		// we truncated below.
148  		roundUp = dfrac>>(extra-1) == 1
149  	}
150  	if dfrac != 0 {
151  		d0 = false
152  	}
153  	// Proceed to the requested number of digits
154  	formatDecimal(d, di, !d0, roundUp, prec)
156  	d.dp -= q
157  }
158
159  var uint64pow10 = [...]uint64{
160  	1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
161  	1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
162  }
163
164  // formatDecimal fills d with at most prec decimal digits
165  // of mantissa m. The boolean trunc indicates whether m
166  // is truncated compared to the original number being formatted.
167  func formatDecimal(d *decimalSlice, m uint64, trunc bool, roundUp bool, prec int) {
168  	max := uint64pow10[prec]
169  	trimmed := 0
170  	for m >= max {
171  		a, b := m/10, m%10
172  		m = a
173  		trimmed++
174  		if b > 5 {
175  			roundUp = true
176  		} else if b < 5 {
177  			roundUp = false
178  		} else { // b == 5
179  			// round up if there are trailing digits,
180  			// or if the new value of m is odd (round-to-even convention)
181  			roundUp = trunc || m&1 == 1
182  		}
183  		if b != 0 {
184  			trunc = true
185  		}
186  	}
187  	if roundUp {
188  		m++
189  	}
190  	if m >= max {
191  		// Happens if di was originally 99999....xx
192  		m /= 10
193  		trimmed++
194  	}
195  	// render digits (similar to formatBits)
196  	n := uint(prec)
197  	d.nd = prec
198  	v := m
199  	for v >= 100 {
200  		var v1, v2 uint64
201  		if v>>32 == 0 {
202  			v1, v2 = uint64(uint32(v)/100), uint64(uint32(v)%100)
203  		} else {
204  			v1, v2 = v/100, v%100
205  		}
206  		n -= 2
207  		d.d[n+1] = smallsString[2*v2+1]
208  		d.d[n+0] = smallsString[2*v2+0]
209  		v = v1
210  	}
211  	if v > 0 {
212  		n--
213  		d.d[n] = smallsString[2*v+1]
214  	}
215  	if v >= 10 {
216  		n--
217  		d.d[n] = smallsString[2*v]
218  	}
219  	for d.d[d.nd-1] == '0' {
220  		d.nd--
221  		trimmed++
222  	}
223  	d.dp = d.nd + trimmed
224  }
225
226  // ryuFtoaShortest formats mant*2^exp with prec decimal digits.
227  func ryuFtoaShortest(d *decimalSlice, mant uint64, exp int, flt *floatInfo) {
228  	if mant == 0 {
229  		d.nd, d.dp = 0, 0
230  		return
231  	}
232  	// If input is an exact integer with fewer bits than the mantissa,
233  	// the previous and next integer are not admissible representations.
234  	if exp <= 0 && bits.TrailingZeros64(mant) >= -exp {
235  		mant >>= uint(-exp)
236  		ryuDigits(d, mant, mant, mant, true, false)
237  		return
238  	}
239  	ml, mc, mu, e2 := computeBounds(mant, exp, flt)
240  	if e2 == 0 {
241  		ryuDigits(d, ml, mc, mu, true, false)
242  		return
243  	}
244  	// Find 10^q *larger* than 2^-e2
245  	q := mulByLog2Log10(-e2) + 1
246
247  	// We are going to multiply by 10^q using 128-bit arithmetic.
248  	// The exponent is the same for all 3 numbers.
249  	var dl, dc, du uint64
250  	var dl0, dc0, du0 bool
251  	if flt == &float32info {
252  		var dl32, dc32, du32 uint32
253  		dl32, _, dl0 = mult64bitPow10(uint32(ml), e2, q)
254  		dc32, _, dc0 = mult64bitPow10(uint32(mc), e2, q)
255  		du32, e2, du0 = mult64bitPow10(uint32(mu), e2, q)
256  		dl, dc, du = uint64(dl32), uint64(dc32), uint64(du32)
257  	} else {
258  		dl, _, dl0 = mult128bitPow10(ml, e2, q)
259  		dc, _, dc0 = mult128bitPow10(mc, e2, q)
260  		du, e2, du0 = mult128bitPow10(mu, e2, q)
261  	}
262  	if e2 >= 0 {
263  		panic("not enough significant bits after mult128bitPow10")
264  	}
265  	// Is it an exact computation?
266  	if q > 55 {
267  		// Large positive powers of ten are not exact
268  		dl0, dc0, du0 = false, false, false
269  	}
270  	if q < 0 && q >= -24 {
271  		// Division by a power of ten may be exact.
272  		// (note that 5^25 is a 59-bit number so division by 5^25 is never exact).
273  		if divisibleByPower5(ml, -q) {
274  			dl0 = true
275  		}
276  		if divisibleByPower5(mc, -q) {
277  			dc0 = true
278  		}
279  		if divisibleByPower5(mu, -q) {
280  			du0 = true
281  		}
282  	}
283  	// Express the results (dl, dc, du)*2^e2 as integers.
284  	// Extra bits must be removed and rounding hints computed.
285  	extra := uint(-e2)
286  	extraMask := uint64(1<<extra - 1)
287  	// Now compute the floored, integral base 10 mantissas.
288  	dl, fracl := dl>>extra, dl&extraMask
289  	dc, fracc := dc>>extra, dc&extraMask
290  	du, fracu := du>>extra, du&extraMask
291  	// Is it allowed to use 'du' as a result?
292  	// It is always allowed when it is truncated, but also
293  	// if it is exact and the original binary mantissa is even
294  	// When disallowed, we can subtract 1.
295  	uok := !du0 || fracu > 0
296  	if du0 && fracu == 0 {
297  		uok = mant&1 == 0
298  	}
299  	if !uok {
300  		du--
301  	}
302  	// Is 'dc' the correctly rounded base 10 mantissa?
303  	// The correct rounding might be dc+1
304  	cup := false // don't round up.
305  	if dc0 {
306  		// If we computed an exact product, the half integer
307  		// should round to next (even) integer if 'dc' is odd.
308  		cup = fracc > 1<<(extra-1) ||
309  			(fracc == 1<<(extra-1) && dc&1 == 1)
310  	} else {
311  		// otherwise, the result is a lower truncation of the ideal
312  		// result.
313  		cup = fracc>>(extra-1) == 1
314  	}
315  	// Is 'dl' an allowed representation?
316  	// Only if it is an exact value, and if the original binary mantissa
317  	// was even.
318  	lok := dl0 && fracl == 0 && (mant&1 == 0)
319  	if !lok {
320  		dl++
321  	}
322  	// We need to remember whether the trimmed digits of 'dc' are zero.
323  	c0 := dc0 && fracc == 0
324  	// render digits
325  	ryuDigits(d, dl, dc, du, c0, cup)
326  	d.dp -= q
327  }
328
329  // mulByLog2Log10 returns math.Floor(x * log(2)/log(10)) for an integer x in
330  // the range -1600 <= x && x <= +1600.
331  //
332  // The range restriction lets us work in faster integer arithmetic instead of
333  // slower floating point arithmetic. Correctness is verified by unit tests.
334  func mulByLog2Log10(x int) int {
335  	// log(2)/log(10) ≈ 0.30102999566 ≈ 78913 / 2^18
336  	return (x * 78913) >> 18
337  }
338
339  // mulByLog10Log2 returns math.Floor(x * log(10)/log(2)) for an integer x in
340  // the range -500 <= x && x <= +500.
341  //
342  // The range restriction lets us work in faster integer arithmetic instead of
343  // slower floating point arithmetic. Correctness is verified by unit tests.
344  func mulByLog10Log2(x int) int {
345  	// log(10)/log(2) ≈ 3.32192809489 ≈ 108853 / 2^15
346  	return (x * 108853) >> 15
347  }
348
349  // computeBounds returns a floating-point vector (l, c, u)×2^e2
350  // where the mantissas are 55-bit (or 26-bit) integers, describing the interval
351  // represented by the input float64 or float32.
352  func computeBounds(mant uint64, exp int, flt *floatInfo) (lower, central, upper uint64, e2 int) {
353  	if mant != 1<<flt.mantbits || exp == flt.bias+1-int(flt.mantbits) {
354  		// regular case (or denormals)
355  		lower, central, upper = 2*mant-1, 2*mant, 2*mant+1
356  		e2 = exp - 1
357  		return
358  	} else {
359  		// border of an exponent
360  		lower, central, upper = 4*mant-1, 4*mant, 4*mant+2
361  		e2 = exp - 2
362  		return
363  	}
364  }
365
366  func ryuDigits(d *decimalSlice, lower, central, upper uint64,
367  	c0, cup bool) {
368  	lhi, llo := divmod1e9(lower)
369  	chi, clo := divmod1e9(central)
370  	uhi, ulo := divmod1e9(upper)
371  	if uhi == 0 {
372  		// only low digits (for denormals)
373  		ryuDigits32(d, llo, clo, ulo, c0, cup, 8)
374  	} else if lhi < uhi {
375  		// truncate 9 digits at once.
376  		if llo != 0 {
377  			lhi++
378  		}
379  		c0 = c0 && clo == 0
380  		cup = (clo > 5e8) || (clo == 5e8 && cup)
381  		ryuDigits32(d, lhi, chi, uhi, c0, cup, 8)
382  		d.dp += 9
383  	} else {
384  		d.nd = 0
385  		// emit high part
386  		n := uint(9)
387  		for v := chi; v > 0; {
388  			v1, v2 := v/10, v%10
389  			v = v1
390  			n--
391  			d.d[n] = byte(v2 + '0')
392  		}
393  		d.d = d.d[n:]
394  		d.nd = int(9 - n)
395  		// emit low part
396  		ryuDigits32(d, llo, clo, ulo,
397  			c0, cup, d.nd+8)
398  	}
399  	// trim trailing zeros
400  	for d.nd > 0 && d.d[d.nd-1] == '0' {
401  		d.nd--
402  	}
403  	// trim initial zeros
404  	for d.nd > 0 && d.d[0] == '0' {
405  		d.nd--
406  		d.dp--
407  		d.d = d.d[1:]
408  	}
409  }
410
411  // ryuDigits32 emits decimal digits for a number less than 1e9.
412  func ryuDigits32(d *decimalSlice, lower, central, upper uint32,
413  	c0, cup bool, endindex int) {
414  	if upper == 0 {
415  		d.dp = endindex + 1
416  		return
417  	}
418  	trimmed := 0
419  	// Remember last trimmed digit to check for round-up.
420  	// c0 will be used to remember zeroness of following digits.
421  	cNextDigit := 0
422  	for upper > 0 {
423  		// Repeatedly compute:
424  		// l = Ceil(lower / 10^k)
425  		// c = Round(central / 10^k)
426  		// u = Floor(upper / 10^k)
427  		// and stop when c goes out of the (l, u) interval.
428  		l := (lower + 9) / 10
429  		c, cdigit := central/10, central%10
430  		u := upper / 10
431  		if l > u {
432  			// don't trim the last digit as it is forbidden to go below l
433  			// other, trim and exit now.
434  			break
435  		}
436  		// Check that we didn't cross the lower boundary.
437  		// The case where l < u but c == l-1 is essentially impossible,
438  		// but may happen if:
439  		//    lower   = ..11
440  		//    central = ..19
441  		//    upper   = ..31
442  		// and means that 'central' is very close but less than
443  		// an integer ending with many zeros, and usually
444  		// the "round-up" logic hides the problem.
445  		if l == c+1 && c < u {
446  			c++
447  			cdigit = 0
448  			cup = false
449  		}
450  		trimmed++
451  		// Remember trimmed digits of c
452  		c0 = c0 && cNextDigit == 0
453  		cNextDigit = int(cdigit)
454  		lower, central, upper = l, c, u
455  	}
456  	// should we round up?
457  	if trimmed > 0 {
458  		cup = cNextDigit > 5 ||
459  			(cNextDigit == 5 && !c0) ||
460  			(cNextDigit == 5 && c0 && central&1 == 1)
461  	}
462  	if central < upper && cup {
463  		central++
464  	}
465  	// We know where the number ends, fill directly
466  	endindex -= trimmed
467  	v := central
468  	n := endindex
469  	for n > d.nd {
470  		v1, v2 := v/100, v%100
471  		d.d[n] = smallsString[2*v2+1]
472  		d.d[n-1] = smallsString[2*v2+0]
473  		n -= 2
474  		v = v1
475  	}
476  	if n == d.nd {
477  		d.d[n] = byte(v + '0')
478  	}
479  	d.nd = endindex + 1
480  	d.dp = d.nd + trimmed
481  }
482
483  // mult64bitPow10 takes a floating-point input with a 25-bit
484  // mantissa and multiplies it with 10^q. The resulting mantissa
485  // is m*P >> 57 where P is a 64-bit element of the detailedPowersOfTen tables.
486  // It is typically 31 or 32-bit wide.
487  // The returned boolean is true if all trimmed bits were zero.
488  //
489  // That is:
490  //
491  //	m*2^e2 * round(10^q) = resM * 2^resE + ε
492  //	exact = ε == 0
493  func mult64bitPow10(m uint32, e2, q int) (resM uint32, resE int, exact bool) {
494  	if q == 0 {
495  		// P == 1<<63
496  		return m << 6, e2 - 6, true
497  	}
498  	if q < detailedPowersOfTenMinExp10 || detailedPowersOfTenMaxExp10 < q {
499  		// This never happens due to the range of float32/float64 exponent
500  		panic("mult64bitPow10: power of 10 is out of range")
501  	}
502  	pow := detailedPowersOfTen[q-detailedPowersOfTenMinExp10][1]
503  	if q < 0 {
504  		// Inverse powers of ten must be rounded up.
505  		pow += 1
506  	}
507  	hi, lo := bits.Mul64(uint64(m), pow)
508  	e2 += mulByLog10Log2(q) - 63 + 57
509  	return uint32(hi<<7 | lo>>57), e2, lo<<7 == 0
510  }
511
512  // mult128bitPow10 takes a floating-point input with a 55-bit
513  // mantissa and multiplies it with 10^q. The resulting mantissa
514  // is m*P >> 119 where P is a 128-bit element of the detailedPowersOfTen tables.
515  // It is typically 63 or 64-bit wide.
516  // The returned boolean is true is all trimmed bits were zero.
517  //
518  // That is:
519  //
520  //	m*2^e2 * round(10^q) = resM * 2^resE + ε
521  //	exact = ε == 0
522  func mult128bitPow10(m uint64, e2, q int) (resM uint64, resE int, exact bool) {
523  	if q == 0 {
524  		// P == 1<<127
525  		return m << 8, e2 - 8, true
526  	}
527  	if q < detailedPowersOfTenMinExp10 || detailedPowersOfTenMaxExp10 < q {
528  		// This never happens due to the range of float32/float64 exponent
529  		panic("mult128bitPow10: power of 10 is out of range")
530  	}
531  	pow := detailedPowersOfTen[q-detailedPowersOfTenMinExp10]
532  	if q < 0 {
533  		// Inverse powers of ten must be rounded up.
534  		pow[0] += 1
535  	}
536  	e2 += mulByLog10Log2(q) - 127 + 119
537
538  	// long multiplication
539  	l1, l0 := bits.Mul64(m, pow[0])
540  	h1, h0 := bits.Mul64(m, pow[1])
541  	mid, carry := bits.Add64(l1, h0, 0)
542  	h1 += carry
543  	return h1<<9 | mid>>55, e2, mid<<9 == 0 && l0 == 0
544  }
545
546  func divisibleByPower5(m uint64, k int) bool {
547  	if m == 0 {
548  		return true
549  	}
550  	for i := 0; i < k; i++ {
551  		if m%5 != 0 {
552  			return false
553  		}
554  		m /= 5
555  	}
556  	return true
557  }
558
559  // divmod1e9 computes quotient and remainder of division by 1e9,
560  // avoiding runtime uint64 division on 32-bit platforms.
561  func divmod1e9(x uint64) (uint32, uint32) {
562  	if !host32bit {
563  		return uint32(x / 1e9), uint32(x % 1e9)
564  	}
565  	// Use the same sequence of operations as the amd64 compiler.
566  	hi, _ := bits.Mul64(x>>1, 0x89705f4136b4a598) // binary digits of 1e-9
567  	q := hi >> 28
568  	return uint32(q), uint32(x - q*1e9)
569  }
570
```

View as plain text