Black Lives Matter. Support the Equal Justice Initiative.

Source file src/cmd/compile/internal/ssa/magic.go

Documentation: cmd/compile/internal/ssa

     1  // Copyright 2016 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 ssa
     6  
     7  import (
     8  	"math/big"
     9  	"math/bits"
    10  )
    11  
    12  // So you want to compute x / c for some constant c?
    13  // Machine division instructions are slow, so we try to
    14  // compute this division with a multiplication + a few
    15  // other cheap instructions instead.
    16  // (We assume here that c != 0, +/- 1, or +/- 2^i.  Those
    17  // cases are easy to handle in different ways).
    18  
    19  // Technique from https://gmplib.org/~tege/divcnst-pldi94.pdf
    20  
    21  // First consider unsigned division.
    22  // Our strategy is to precompute 1/c then do
    23  //   ⎣x / c⎦ = ⎣x * (1/c)⎦.
    24  // 1/c is less than 1, so we can't compute it directly in
    25  // integer arithmetic.  Let's instead compute 2^e/c
    26  // for a value of e TBD (^ = exponentiation).  Then
    27  //   ⎣x / c⎦ = ⎣x * (2^e/c) / 2^e⎦.
    28  // Dividing by 2^e is easy.  2^e/c isn't an integer, unfortunately.
    29  // So we must approximate it.  Let's call its approximation m.
    30  // We'll then compute
    31  //   ⎣x * m / 2^e⎦
    32  // Which we want to be equal to ⎣x / c⎦ for 0 <= x < 2^n-1
    33  // where n is the word size.
    34  // Setting x = c gives us c * m >= 2^e.
    35  // We'll chose m = ⎡2^e/c⎤ to satisfy that equation.
    36  // What remains is to choose e.
    37  // Let m = 2^e/c + delta, 0 <= delta < 1
    38  //   ⎣x * (2^e/c + delta) / 2^e⎦
    39  //   ⎣x / c + x * delta / 2^e⎦
    40  // We must have x * delta / 2^e < 1/c so that this
    41  // additional term never rounds differently than ⎣x / c⎦ does.
    42  // Rearranging,
    43  //   2^e > x * delta * c
    44  // x can be at most 2^n-1 and delta can be at most 1.
    45  // So it is sufficient to have 2^e >= 2^n*c.
    46  // So we'll choose e = n + s, with s = ⎡log2(c)⎤.
    47  //
    48  // An additional complication arises because m has n+1 bits in it.
    49  // Hardware restricts us to n bit by n bit multiplies.
    50  // We divide into 3 cases:
    51  //
    52  // Case 1: m is even.
    53  //   ⎣x / c⎦ = ⎣x * m / 2^(n+s)⎦
    54  //   ⎣x / c⎦ = ⎣x * (m/2) / 2^(n+s-1)⎦
    55  //   ⎣x / c⎦ = ⎣x * (m/2) / 2^n / 2^(s-1)⎦
    56  //   ⎣x / c⎦ = ⎣⎣x * (m/2) / 2^n⎦ / 2^(s-1)⎦
    57  //   multiply + shift
    58  //
    59  // Case 2: c is even.
    60  //   ⎣x / c⎦ = ⎣(x/2) / (c/2)⎦
    61  //   ⎣x / c⎦ = ⎣⎣x/2⎦ / (c/2)⎦
    62  //     This is just the original problem, with x' = ⎣x/2⎦, c' = c/2, n' = n-1.
    63  //       s' = s-1
    64  //       m' = ⎡2^(n'+s')/c'⎤
    65  //          = ⎡2^(n+s-1)/c⎤
    66  //          = ⎡m/2⎤
    67  //   ⎣x / c⎦ = ⎣x' * m' / 2^(n'+s')⎦
    68  //   ⎣x / c⎦ = ⎣⎣x/2⎦ * ⎡m/2⎤ / 2^(n+s-2)⎦
    69  //   ⎣x / c⎦ = ⎣⎣⎣x/2⎦ * ⎡m/2⎤ / 2^n⎦ / 2^(s-2)⎦
    70  //   shift + multiply + shift
    71  //
    72  // Case 3: everything else
    73  //   let k = m - 2^n. k fits in n bits.
    74  //   ⎣x / c⎦ = ⎣x * m / 2^(n+s)⎦
    75  //   ⎣x / c⎦ = ⎣x * (2^n + k) / 2^(n+s)⎦
    76  //   ⎣x / c⎦ = ⎣(x + x * k / 2^n) / 2^s⎦
    77  //   ⎣x / c⎦ = ⎣(x + ⎣x * k / 2^n⎦) / 2^s⎦
    78  //   ⎣x / c⎦ = ⎣(x + ⎣x * k / 2^n⎦) / 2^s⎦
    79  //   ⎣x / c⎦ = ⎣⎣(x + ⎣x * k / 2^n⎦) / 2⎦ / 2^(s-1)⎦
    80  //   multiply + avg + shift
    81  //
    82  // These can be implemented in hardware using:
    83  //  ⎣a * b / 2^n⎦ - aka high n bits of an n-bit by n-bit multiply.
    84  //  ⎣(a+b) / 2⎦   - aka "average" of two n-bit numbers.
    85  //                  (Not just a regular add & shift because the intermediate result
    86  //                   a+b has n+1 bits in it.  Nevertheless, can be done
    87  //                   in 2 instructions on x86.)
    88  
    89  // umagicOK reports whether we should strength reduce a n-bit divide by c.
    90  func umagicOK(n uint, c int64) bool {
    91  	// Convert from ConstX auxint values to the real uint64 constant they represent.
    92  	d := uint64(c) << (64 - n) >> (64 - n)
    93  
    94  	// Doesn't work for 0.
    95  	// Don't use for powers of 2.
    96  	return d&(d-1) != 0
    97  }
    98  
    99  // umagicOKn reports whether we should strength reduce an unsigned n-bit divide by c.
   100  // We can strength reduce when c != 0 and c is not a power of two.
   101  func umagicOK8(c int8) bool   { return c&(c-1) != 0 }
   102  func umagicOK16(c int16) bool { return c&(c-1) != 0 }
   103  func umagicOK32(c int32) bool { return c&(c-1) != 0 }
   104  func umagicOK64(c int64) bool { return c&(c-1) != 0 }
   105  
   106  type umagicData struct {
   107  	s int64  // ⎡log2(c)⎤
   108  	m uint64 // ⎡2^(n+s)/c⎤ - 2^n
   109  }
   110  
   111  // umagic computes the constants needed to strength reduce unsigned n-bit divides by the constant uint64(c).
   112  // The return values satisfy for all 0 <= x < 2^n
   113  //  floor(x / uint64(c)) = x * (m + 2^n) >> (n+s)
   114  func umagic(n uint, c int64) umagicData {
   115  	// Convert from ConstX auxint values to the real uint64 constant they represent.
   116  	d := uint64(c) << (64 - n) >> (64 - n)
   117  
   118  	C := new(big.Int).SetUint64(d)
   119  	s := C.BitLen()
   120  	M := big.NewInt(1)
   121  	M.Lsh(M, n+uint(s))     // 2^(n+s)
   122  	M.Add(M, C)             // 2^(n+s)+c
   123  	M.Sub(M, big.NewInt(1)) // 2^(n+s)+c-1
   124  	M.Div(M, C)             // ⎡2^(n+s)/c⎤
   125  	if M.Bit(int(n)) != 1 {
   126  		panic("n+1st bit isn't set")
   127  	}
   128  	M.SetBit(M, int(n), 0)
   129  	m := M.Uint64()
   130  	return umagicData{s: int64(s), m: m}
   131  }
   132  
   133  func umagic8(c int8) umagicData   { return umagic(8, int64(c)) }
   134  func umagic16(c int16) umagicData { return umagic(16, int64(c)) }
   135  func umagic32(c int32) umagicData { return umagic(32, int64(c)) }
   136  func umagic64(c int64) umagicData { return umagic(64, c) }
   137  
   138  // For signed division, we use a similar strategy.
   139  // First, we enforce a positive c.
   140  //   x / c = -(x / (-c))
   141  // This will require an additional Neg op for c<0.
   142  //
   143  // If x is positive we're in a very similar state
   144  // to the unsigned case above.  We define:
   145  //   s = ⎡log2(c)⎤-1
   146  //   m = ⎡2^(n+s)/c⎤
   147  // Then
   148  //   ⎣x / c⎦ = ⎣x * m / 2^(n+s)⎦
   149  // If x is negative we have
   150  //   ⎡x / c⎤ = ⎣x * m / 2^(n+s)⎦ + 1
   151  // (TODO: derivation?)
   152  //
   153  // The multiply is a bit odd, as it is a signed n-bit value
   154  // times an unsigned n-bit value.  For n smaller than the
   155  // word size, we can extend x and m appropriately and use the
   156  // signed multiply instruction.  For n == word size,
   157  // we must use the signed multiply high and correct
   158  // the result by adding x*2^n.
   159  //
   160  // Adding 1 if x<0 is done by subtracting x>>(n-1).
   161  
   162  func smagicOK(n uint, c int64) bool {
   163  	if c < 0 {
   164  		// Doesn't work for negative c.
   165  		return false
   166  	}
   167  	// Doesn't work for 0.
   168  	// Don't use it for powers of 2.
   169  	return c&(c-1) != 0
   170  }
   171  
   172  // smagicOKn reports whether we should strength reduce an signed n-bit divide by c.
   173  func smagicOK8(c int8) bool   { return smagicOK(8, int64(c)) }
   174  func smagicOK16(c int16) bool { return smagicOK(16, int64(c)) }
   175  func smagicOK32(c int32) bool { return smagicOK(32, int64(c)) }
   176  func smagicOK64(c int64) bool { return smagicOK(64, c) }
   177  
   178  type smagicData struct {
   179  	s int64  // ⎡log2(c)⎤-1
   180  	m uint64 // ⎡2^(n+s)/c⎤
   181  }
   182  
   183  // magic computes the constants needed to strength reduce signed n-bit divides by the constant c.
   184  // Must have c>0.
   185  // The return values satisfy for all -2^(n-1) <= x < 2^(n-1)
   186  //  trunc(x / c) = x * m >> (n+s) + (x < 0 ? 1 : 0)
   187  func smagic(n uint, c int64) smagicData {
   188  	C := new(big.Int).SetInt64(c)
   189  	s := C.BitLen() - 1
   190  	M := big.NewInt(1)
   191  	M.Lsh(M, n+uint(s))     // 2^(n+s)
   192  	M.Add(M, C)             // 2^(n+s)+c
   193  	M.Sub(M, big.NewInt(1)) // 2^(n+s)+c-1
   194  	M.Div(M, C)             // ⎡2^(n+s)/c⎤
   195  	if M.Bit(int(n)) != 0 {
   196  		panic("n+1st bit is set")
   197  	}
   198  	if M.Bit(int(n-1)) == 0 {
   199  		panic("nth bit is not set")
   200  	}
   201  	m := M.Uint64()
   202  	return smagicData{s: int64(s), m: m}
   203  }
   204  
   205  func smagic8(c int8) smagicData   { return smagic(8, int64(c)) }
   206  func smagic16(c int16) smagicData { return smagic(16, int64(c)) }
   207  func smagic32(c int32) smagicData { return smagic(32, int64(c)) }
   208  func smagic64(c int64) smagicData { return smagic(64, c) }
   209  
   210  // Divisibility x%c == 0 can be checked more efficiently than directly computing
   211  // the modulus x%c and comparing against 0.
   212  //
   213  // The same "Division by invariant integers using multiplication" paper
   214  // by Granlund and Montgomery referenced above briefly mentions this method
   215  // and it is further elaborated in "Hacker's Delight" by Warren Section 10-17
   216  //
   217  // The first thing to note is that for odd integers, exact division can be computed
   218  // by using the modular inverse with respect to the word size 2^n.
   219  //
   220  // Given c, compute m such that (c * m) mod 2^n == 1
   221  // Then if c divides x (x%c ==0), the quotient is given by q = x/c == x*m mod 2^n
   222  //
   223  // x can range from 0, c, 2c, 3c, ... ⎣(2^n - 1)/c⎦ * c the maximum multiple
   224  // Thus, x*m mod 2^n is 0, 1, 2, 3, ... ⎣(2^n - 1)/c⎦
   225  // i.e. the quotient takes all values from zero up to max = ⎣(2^n - 1)/c⎦
   226  //
   227  // If x is not divisible by c, then x*m mod 2^n must take some larger value than max.
   228  //
   229  // This gives x*m mod 2^n <= ⎣(2^n - 1)/c⎦ as a test for divisibility
   230  // involving one multiplication and compare.
   231  //
   232  // To extend this to even integers, consider c = d0 * 2^k where d0 is odd.
   233  // We can test whether x is divisible by both d0 and 2^k.
   234  // For d0, the test is the same as above.  Let m be such that m*d0 mod 2^n == 1
   235  // Then x*m mod 2^n <= ⎣(2^n - 1)/d0⎦ is the first test.
   236  // The test for divisibility by 2^k is a check for k trailing zeroes.
   237  // Note that since d0 is odd, m is odd and thus x*m will have the same number of
   238  // trailing zeroes as x.  So the two tests are,
   239  //
   240  // x*m mod 2^n <= ⎣(2^n - 1)/d0⎦
   241  // and x*m ends in k zero bits
   242  //
   243  // These can be combined into a single comparison by the following
   244  // (theorem ZRU in Hacker's Delight) for unsigned integers.
   245  //
   246  // x <= a and x ends in k zero bits if and only if RotRight(x ,k) <= ⎣a/(2^k)⎦
   247  // Where RotRight(x ,k) is right rotation of x by k bits.
   248  //
   249  // To prove the first direction, x <= a -> ⎣x/(2^k)⎦ <= ⎣a/(2^k)⎦
   250  // But since x ends in k zeroes all the rotated bits would be zero too.
   251  // So RotRight(x, k) == ⎣x/(2^k)⎦ <= ⎣a/(2^k)⎦
   252  //
   253  // If x does not end in k zero bits, then RotRight(x, k)
   254  // has some non-zero bits in the k highest bits.
   255  // ⎣x/(2^k)⎦ has all zeroes in the k highest bits,
   256  // so RotRight(x, k) > ⎣x/(2^k)⎦
   257  //
   258  // Finally, if x > a and has k trailing zero bits, then RotRight(x, k) == ⎣x/(2^k)⎦
   259  // and ⎣x/(2^k)⎦ must be greater than ⎣a/(2^k)⎦, that is the top n-k bits of x must
   260  // be greater than the top n-k bits of a because the rest of x bits are zero.
   261  //
   262  // So the two conditions about can be replaced with the single test
   263  //
   264  // RotRight(x*m mod 2^n, k) <= ⎣(2^n - 1)/c⎦
   265  //
   266  // Where d0*2^k was replaced by c on the right hand side.
   267  
   268  // udivisibleOK reports whether we should strength reduce an unsigned n-bit divisibilty check by c.
   269  func udivisibleOK(n uint, c int64) bool {
   270  	// Convert from ConstX auxint values to the real uint64 constant they represent.
   271  	d := uint64(c) << (64 - n) >> (64 - n)
   272  
   273  	// Doesn't work for 0.
   274  	// Don't use for powers of 2.
   275  	return d&(d-1) != 0
   276  }
   277  
   278  func udivisibleOK8(c int8) bool   { return udivisibleOK(8, int64(c)) }
   279  func udivisibleOK16(c int16) bool { return udivisibleOK(16, int64(c)) }
   280  func udivisibleOK32(c int32) bool { return udivisibleOK(32, int64(c)) }
   281  func udivisibleOK64(c int64) bool { return udivisibleOK(64, c) }
   282  
   283  type udivisibleData struct {
   284  	k   int64  // trailingZeros(c)
   285  	m   uint64 // m * (c>>k) mod 2^n == 1 multiplicative inverse of odd portion modulo 2^n
   286  	max uint64 // ⎣(2^n - 1)/ c⎦ max value to for divisibility
   287  }
   288  
   289  func udivisible(n uint, c int64) udivisibleData {
   290  	// Convert from ConstX auxint values to the real uint64 constant they represent.
   291  	d := uint64(c) << (64 - n) >> (64 - n)
   292  
   293  	k := bits.TrailingZeros64(d)
   294  	d0 := d >> uint(k) // the odd portion of the divisor
   295  
   296  	mask := ^uint64(0) >> (64 - n)
   297  
   298  	// Calculate the multiplicative inverse via Newton's method.
   299  	// Quadratic convergence doubles the number of correct bits per iteration.
   300  	m := d0            // initial guess correct to 3-bits d0*d0 mod 8 == 1
   301  	m = m * (2 - m*d0) // 6-bits
   302  	m = m * (2 - m*d0) // 12-bits
   303  	m = m * (2 - m*d0) // 24-bits
   304  	m = m * (2 - m*d0) // 48-bits
   305  	m = m * (2 - m*d0) // 96-bits >= 64-bits
   306  	m = m & mask
   307  
   308  	max := mask / d
   309  
   310  	return udivisibleData{
   311  		k:   int64(k),
   312  		m:   m,
   313  		max: max,
   314  	}
   315  }
   316  
   317  func udivisible8(c int8) udivisibleData   { return udivisible(8, int64(c)) }
   318  func udivisible16(c int16) udivisibleData { return udivisible(16, int64(c)) }
   319  func udivisible32(c int32) udivisibleData { return udivisible(32, int64(c)) }
   320  func udivisible64(c int64) udivisibleData { return udivisible(64, c) }
   321  
   322  // For signed integers, a similar method follows.
   323  //
   324  // Given c > 1 and odd, compute m such that (c * m) mod 2^n == 1
   325  // Then if c divides x (x%c ==0), the quotient is given by q = x/c == x*m mod 2^n
   326  //
   327  // x can range from ⎡-2^(n-1)/c⎤ * c, ... -c, 0, c, ...  ⎣(2^(n-1) - 1)/c⎦ * c
   328  // Thus, x*m mod 2^n is ⎡-2^(n-1)/c⎤, ... -2, -1, 0, 1, 2, ... ⎣(2^(n-1) - 1)/c⎦
   329  //
   330  // So, x is a multiple of c if and only if:
   331  // ⎡-2^(n-1)/c⎤ <= x*m mod 2^n <= ⎣(2^(n-1) - 1)/c⎦
   332  //
   333  // Since c > 1 and odd, this can be simplified by
   334  // ⎡-2^(n-1)/c⎤ == ⎡(-2^(n-1) + 1)/c⎤ == -⎣(2^(n-1) - 1)/c⎦
   335  //
   336  // -⎣(2^(n-1) - 1)/c⎦ <= x*m mod 2^n <= ⎣(2^(n-1) - 1)/c⎦
   337  //
   338  // To extend this to even integers, consider c = d0 * 2^k where d0 is odd.
   339  // We can test whether x is divisible by both d0 and 2^k.
   340  //
   341  // Let m be such that (d0 * m) mod 2^n == 1.
   342  // Let q = x*m mod 2^n. Then c divides x if:
   343  //
   344  // -⎣(2^(n-1) - 1)/d0⎦ <= q <= ⎣(2^(n-1) - 1)/d0⎦ and q ends in at least k 0-bits
   345  //
   346  // To transform this to a single comparison, we use the following theorem (ZRS in Hacker's Delight).
   347  //
   348  // For a >= 0 the following conditions are equivalent:
   349  // 1) -a <= x <= a and x ends in at least k 0-bits
   350  // 2) RotRight(x+a', k) <= ⎣2a'/2^k⎦
   351  //
   352  // Where a' = a & -2^k (a with its right k bits set to zero)
   353  //
   354  // To see that 1 & 2 are equivalent, note that -a <= x <= a is equivalent to
   355  // -a' <= x <= a' if and only if x ends in at least k 0-bits.  Adding -a' to each side gives,
   356  // 0 <= x + a' <= 2a' and x + a' ends in at least k 0-bits if and only if x does since a' has
   357  // k 0-bits by definition.  We can use theorem ZRU above with x -> x + a' and a -> 2a' giving 1) == 2).
   358  //
   359  // Let m be such that (d0 * m) mod 2^n == 1.
   360  // Let q = x*m mod 2^n.
   361  // Let a' = ⎣(2^(n-1) - 1)/d0⎦ & -2^k
   362  //
   363  // Then the divisibility test is:
   364  //
   365  // RotRight(q+a', k) <= ⎣2a'/2^k⎦
   366  //
   367  // Note that the calculation is performed using unsigned integers.
   368  // Since a' can have n-1 bits, 2a' may have n bits and there is no risk of overflow.
   369  
   370  // sdivisibleOK reports whether we should strength reduce a signed n-bit divisibilty check by c.
   371  func sdivisibleOK(n uint, c int64) bool {
   372  	if c < 0 {
   373  		// Doesn't work for negative c.
   374  		return false
   375  	}
   376  	// Doesn't work for 0.
   377  	// Don't use it for powers of 2.
   378  	return c&(c-1) != 0
   379  }
   380  
   381  func sdivisibleOK8(c int8) bool   { return sdivisibleOK(8, int64(c)) }
   382  func sdivisibleOK16(c int16) bool { return sdivisibleOK(16, int64(c)) }
   383  func sdivisibleOK32(c int32) bool { return sdivisibleOK(32, int64(c)) }
   384  func sdivisibleOK64(c int64) bool { return sdivisibleOK(64, c) }
   385  
   386  type sdivisibleData struct {
   387  	k   int64  // trailingZeros(c)
   388  	m   uint64 // m * (c>>k) mod 2^n == 1 multiplicative inverse of odd portion modulo 2^n
   389  	a   uint64 // ⎣(2^(n-1) - 1)/ (c>>k)⎦ & -(1<<k) additive constant
   390  	max uint64 // ⎣(2 a) / (1<<k)⎦ max value to for divisibility
   391  }
   392  
   393  func sdivisible(n uint, c int64) sdivisibleData {
   394  	d := uint64(c)
   395  	k := bits.TrailingZeros64(d)
   396  	d0 := d >> uint(k) // the odd portion of the divisor
   397  
   398  	mask := ^uint64(0) >> (64 - n)
   399  
   400  	// Calculate the multiplicative inverse via Newton's method.
   401  	// Quadratic convergence doubles the number of correct bits per iteration.
   402  	m := d0            // initial guess correct to 3-bits d0*d0 mod 8 == 1
   403  	m = m * (2 - m*d0) // 6-bits
   404  	m = m * (2 - m*d0) // 12-bits
   405  	m = m * (2 - m*d0) // 24-bits
   406  	m = m * (2 - m*d0) // 48-bits
   407  	m = m * (2 - m*d0) // 96-bits >= 64-bits
   408  	m = m & mask
   409  
   410  	a := ((mask >> 1) / d0) & -(1 << uint(k))
   411  	max := (2 * a) >> uint(k)
   412  
   413  	return sdivisibleData{
   414  		k:   int64(k),
   415  		m:   m,
   416  		a:   a,
   417  		max: max,
   418  	}
   419  }
   420  
   421  func sdivisible8(c int8) sdivisibleData   { return sdivisible(8, int64(c)) }
   422  func sdivisible16(c int16) sdivisibleData { return sdivisible(16, int64(c)) }
   423  func sdivisible32(c int32) sdivisibleData { return sdivisible(32, int64(c)) }
   424  func sdivisible64(c int64) sdivisibleData { return sdivisible(64, c) }
   425  

View as plain text