math.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  // A uint128 is a 128-bit uint.
    10  // The fields are exported to make them visible to package strconv_test.
    11  type uint128 struct {
    12  	Hi uint64
    13  	Lo uint64
    14  }
    15  
    16  // umul128 returns the 128-bit product x*y.
    17  func umul128(x, y uint64) uint128 {
    18  	hi, lo := bits.Mul64(x, y)
    19  	return uint128{hi, lo}
    20  }
    21  
    22  // umul192 returns the 192-bit product x*y in three uint64s.
    23  func umul192(x uint64, y uint128) (hi, mid, lo uint64) {
    24  	mid1, lo := bits.Mul64(x, y.Lo)
    25  	hi, mid2 := bits.Mul64(x, y.Hi)
    26  	mid, carry := bits.Add64(mid1, mid2, 0)
    27  	return hi + carry, mid, lo
    28  }
    29  
    30  // pow10 returns the 128-bit mantissa and binary exponent of 10**e.
    31  // That is, 10^e = mant/2^128 * 2**exp.
    32  // If e is out of range, pow10 returns ok=false.
    33  func pow10(e int) (mant uint128, exp int, ok bool) {
    34  	if e < pow10Min || e > pow10Max {
    35  		return
    36  	}
    37  	return pow10Tab[e-pow10Min], 1 + mulLog2_10(e), true
    38  }
    39  
    40  // mulLog10_2 returns math.Floor(x * log(2)/log(10)) for an integer x in
    41  // the range -1600 <= x && x <= +1600.
    42  //
    43  // The range restriction lets us work in faster integer arithmetic instead of
    44  // slower floating point arithmetic. Correctness is verified by unit tests.
    45  func mulLog10_2(x int) int {
    46  	// log(2)/log(10) ≈ 0.30102999566 ≈ 78913 / 2^18
    47  	return (x * 78913) >> 18
    48  }
    49  
    50  // mulLog2_10 returns math.Floor(x * log(10)/log(2)) for an integer x in
    51  // the range -500 <= x && x <= +500.
    52  //
    53  // The range restriction lets us work in faster integer arithmetic instead of
    54  // slower floating point arithmetic. Correctness is verified by unit tests.
    55  func mulLog2_10(x int) int {
    56  	// log(10)/log(2) ≈ 3.32192809489 ≈ 108853 / 2^15
    57  	return (x * 108853) >> 15
    58  }
    59  
    60  func bool2uint(b bool) uint {
    61  	if b {
    62  		return 1
    63  	}
    64  	return 0
    65  }
    66  
    67  // Exact Division and Remainder Checking
    68  //
    69  // An exact division x/c (exact means x%c == 0)
    70  // can be implemented by x*m where m is the multiplicative inverse of c (m*c == 1).
    71  //
    72  // Since c is also the multiplicative inverse of m, x*m is lossless,
    73  // and all the exact multiples of c map to all of [0, maxUint64/c].
    74  // The non-multiples are forced to map to larger values.
    75  // This also gives a quick test for whether x is an exact multiple of c:
    76  // compute the exact division and check whether it's at most maxUint64/c:
    77  //	x%c == 0 => x*m <= maxUint64/c.
    78  //
    79  // Only odd c have multiplicative inverses mod powers of two.
    80  // To do an exact divide x / (c<<s) we can use (x/c)>>s instead.
    81  // And to check for remainder, we need to check that those low s
    82  // bits are all zero before we shift them away. We can merge that
    83  // with the <= for the exact odd remainder check by rotating the
    84  // shifted bits into the high part instead:
    85  // 	x%(c<<s) == 0 => bits.RotateLeft64(x*m, -s) <= maxUint64/c.
    86  //
    87  // The compiler does this transformation automatically in general,
    88  // but we apply it here by hand in a few ways that the compiler can't help with.
    89  //
    90  // For a more detailed explanation, see
    91  // Henry S. Warren, Jr., Hacker's Delight, 2nd ed., sections 10-16 and 10-17.
    92  
    93  // divisiblePow5 reports whether x is divisible by 5^p.
    94  // It returns false for p not in [1, 22],
    95  // because we only care about float64 mantissas, and 5^23 > 2^53.
    96  func divisiblePow5(x uint64, p int) bool {
    97  	return 1 <= p && p <= 22 && x*div5Tab[p-1][0] <= div5Tab[p-1][1]
    98  }
    99  
   100  const maxUint64 = 1<<64 - 1
   101  
   102  // div5Tab[p-1] is the multiplicative inverse of 5^p and maxUint64/5^p.
   103  var div5Tab = [22][2]uint64{
   104  	{0xcccccccccccccccd, maxUint64 / 5},
   105  	{0x8f5c28f5c28f5c29, maxUint64 / 5 / 5},
   106  	{0x1cac083126e978d5, maxUint64 / 5 / 5 / 5},
   107  	{0xd288ce703afb7e91, maxUint64 / 5 / 5 / 5 / 5},
   108  	{0x5d4e8fb00bcbe61d, maxUint64 / 5 / 5 / 5 / 5 / 5},
   109  	{0x790fb65668c26139, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5},
   110  	{0xe5032477ae8d46a5, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   111  	{0xc767074b22e90e21, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   112  	{0x8e47ce423a2e9c6d, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   113  	{0x4fa7f60d3ed61f49, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   114  	{0x0fee64690c913975, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   115  	{0x3662e0e1cf503eb1, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   116  	{0xa47a2cf9f6433fbd, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   117  	{0x54186f653140a659, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   118  	{0x7738164770402145, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   119  	{0xe4a4d1417cd9a041, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   120  	{0xc75429d9e5c5200d, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   121  	{0xc1773b91fac10669, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   122  	{0x26b172506559ce15, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   123  	{0xd489e3a9addec2d1, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   124  	{0x90e860bb892c8d5d, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   125  	{0x502e79bf1b6f4f79, maxUint64 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5 / 5},
   126  }
   127  
   128  // trimZeros trims trailing zeros from x.
   129  // It finds the largest p such that x % 10^p == 0
   130  // and then returns x / 10^p, p.
   131  //
   132  // This is here for reference and tested, because it is an optimization
   133  // used by other ftoa algorithms, but in our implementations it has
   134  // never been benchmarked to be faster than trimming zeros after
   135  // formatting into decimal bytes.
   136  func trimZeros(x uint64) (uint64, int) {
   137  	const (
   138  		div1e8m  = 0xc767074b22e90e21
   139  		div1e8le = maxUint64 / 100000000
   140  
   141  		div1e4m  = 0xd288ce703afb7e91
   142  		div1e4le = maxUint64 / 10000
   143  
   144  		div1e2m  = 0x8f5c28f5c28f5c29
   145  		div1e2le = maxUint64 / 100
   146  
   147  		div1e1m  = 0xcccccccccccccccd
   148  		div1e1le = maxUint64 / 10
   149  	)
   150  
   151  	// _ = assert[x - y] asserts at compile time that x == y.
   152  	// Assert that the multiplicative inverses are correct
   153  	// by checking that (div1eNm * 5^N) % 1<<64 == 1.
   154  	var assert [1]struct{}
   155  	_ = assert[(div1e8m*5*5*5*5*5*5*5*5)%(1<<64)-1]
   156  	_ = assert[(div1e4m*5*5*5*5)%(1<<64)-1]
   157  	_ = assert[(div1e2m*5*5)%(1<<64)-1]
   158  	_ = assert[(div1e1m*5)%(1<<64)-1]
   159  
   160  	// Cut 8 zeros, then 4, then 2, then 1.
   161  	p := 0
   162  	for d := bits.RotateLeft64(x*div1e8m, -8); d <= div1e8le; d = bits.RotateLeft64(x*div1e8m, -8) {
   163  		x = d
   164  		p += 8
   165  	}
   166  	if d := bits.RotateLeft64(x*div1e4m, -4); d <= div1e4le {
   167  		x = d
   168  		p += 4
   169  	}
   170  	if d := bits.RotateLeft64(x*div1e2m, -2); d <= div1e2le {
   171  		x = d
   172  		p += 2
   173  	}
   174  	if d := bits.RotateLeft64(x*div1e1m, -1); d <= div1e1le {
   175  		x = d
   176  		p += 1
   177  	}
   178  	return x, p
   179  }
   180
View as plain text