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
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)
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)
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