Black Lives Matter. Support the Equal Justice Initiative.

# Source file src/math/big/arith.go

## Documentation: math/big

```     1  // Copyright 2009 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
4
5  // This file provides Go implementations of elementary multi-precision
6  // arithmetic operations on word vectors. These have the suffix _g.
7  // These are needed for platforms without assembly implementations of these routines.
8  // This file also contains elementary operations that can be implemented
9  // sufficiently efficiently in Go.
10
11  package big
12
13  import "math/bits"
14
15  // A Word represents a single digit of a multi-precision unsigned integer.
16  type Word uint
17
18  const (
19  	_S = _W / 8 // word size in bytes
20
21  	_W = bits.UintSize // word size in bits
22  	_B = 1 << _W       // digit base
23  	_M = _B - 1        // digit mask
24  )
25
26  // Many of the loops in this file are of the form
27  //   for i := 0; i < len(z) && i < len(x) && i < len(y); i++
28  // i < len(z) is the real condition.
29  // However, checking i < len(x) && i < len(y) as well is faster than
30  // having the compiler do a bounds check in the body of the loop;
31  // remarkably it is even faster than hoisting the bounds check
32  // out of the loop, by doing something like
33  //   _, _ = x[len(z)-1], y[len(z)-1]
34  // There are other ways to hoist the bounds check out of the loop,
35  // but the compiler's BCE isn't powerful enough for them (yet?).
36  // See the discussion in CL 164966.
37
38  // ----------------------------------------------------------------------------
39  // Elementary operations on words
40  //
41  // These operations are used by the vector operations below.
42
43  // z1<<_W + z0 = x*y
44  func mulWW_g(x, y Word) (z1, z0 Word) {
45  	hi, lo := bits.Mul(uint(x), uint(y))
46  	return Word(hi), Word(lo)
47  }
48
49  // z1<<_W + z0 = x*y + c
50  func mulAddWWW_g(x, y, c Word) (z1, z0 Word) {
51  	hi, lo := bits.Mul(uint(x), uint(y))
52  	var cc uint
53  	lo, cc = bits.Add(lo, uint(c), 0)
54  	return Word(hi + cc), Word(lo)
55  }
56
57  // nlz returns the number of leading zeros in x.
58  // Wraps bits.LeadingZeros call for convenience.
59  func nlz(x Word) uint {
61  }
62
63  // The resulting carry c is either 0 or 1.
64  func addVV_g(z, x, y []Word) (c Word) {
65  	// The comment near the top of this file discusses this for loop condition.
66  	for i := 0; i < len(z) && i < len(x) && i < len(y); i++ {
67  		zi, cc := bits.Add(uint(x[i]), uint(y[i]), uint(c))
68  		z[i] = Word(zi)
69  		c = Word(cc)
70  	}
71  	return
72  }
73
74  // The resulting carry c is either 0 or 1.
75  func subVV_g(z, x, y []Word) (c Word) {
76  	// The comment near the top of this file discusses this for loop condition.
77  	for i := 0; i < len(z) && i < len(x) && i < len(y); i++ {
78  		zi, cc := bits.Sub(uint(x[i]), uint(y[i]), uint(c))
79  		z[i] = Word(zi)
80  		c = Word(cc)
81  	}
82  	return
83  }
84
85  // The resulting carry c is either 0 or 1.
86  func addVW_g(z, x []Word, y Word) (c Word) {
87  	c = y
88  	// The comment near the top of this file discusses this for loop condition.
89  	for i := 0; i < len(z) && i < len(x); i++ {
90  		zi, cc := bits.Add(uint(x[i]), uint(c), 0)
91  		z[i] = Word(zi)
92  		c = Word(cc)
93  	}
94  	return
95  }
96
98  // The only difference is that we check on every iteration
99  // whether we are done with carries,
100  // and if so, switch to a much faster copy instead.
101  // This is only a good idea for large z,
102  // because the overhead of the check and the function call
103  // outweigh the benefits when z is small.
104  func addVWlarge(z, x []Word, y Word) (c Word) {
105  	c = y
106  	// The comment near the top of this file discusses this for loop condition.
107  	for i := 0; i < len(z) && i < len(x); i++ {
108  		if c == 0 {
109  			copy(z[i:], x[i:])
110  			return
111  		}
112  		zi, cc := bits.Add(uint(x[i]), uint(c), 0)
113  		z[i] = Word(zi)
114  		c = Word(cc)
115  	}
116  	return
117  }
118
119  func subVW_g(z, x []Word, y Word) (c Word) {
120  	c = y
121  	// The comment near the top of this file discusses this for loop condition.
122  	for i := 0; i < len(z) && i < len(x); i++ {
123  		zi, cc := bits.Sub(uint(x[i]), uint(c), 0)
124  		z[i] = Word(zi)
125  		c = Word(cc)
126  	}
127  	return
128  }
129
131  func subVWlarge(z, x []Word, y Word) (c Word) {
132  	c = y
133  	// The comment near the top of this file discusses this for loop condition.
134  	for i := 0; i < len(z) && i < len(x); i++ {
135  		if c == 0 {
136  			copy(z[i:], x[i:])
137  			return
138  		}
139  		zi, cc := bits.Sub(uint(x[i]), uint(c), 0)
140  		z[i] = Word(zi)
141  		c = Word(cc)
142  	}
143  	return
144  }
145
146  func shlVU_g(z, x []Word, s uint) (c Word) {
147  	if s == 0 {
148  		copy(z, x)
149  		return
150  	}
151  	if len(z) == 0 {
152  		return
153  	}
154  	s &= _W - 1 // hint to the compiler that shifts by s don't need guard code
155  	ŝ := _W - s
156  	ŝ &= _W - 1 // ditto
157  	c = x[len(z)-1] >> ŝ
158  	for i := len(z) - 1; i > 0; i-- {
159  		z[i] = x[i]<<s | x[i-1]>>ŝ
160  	}
161  	z[0] = x[0] << s
162  	return
163  }
164
165  func shrVU_g(z, x []Word, s uint) (c Word) {
166  	if s == 0 {
167  		copy(z, x)
168  		return
169  	}
170  	if len(z) == 0 {
171  		return
172  	}
173  	if len(x) != len(z) {
174  		// This is an invariant guaranteed by the caller.
175  		panic("len(x) != len(z)")
176  	}
177  	s &= _W - 1 // hint to the compiler that shifts by s don't need guard code
178  	ŝ := _W - s
179  	ŝ &= _W - 1 // ditto
180  	c = x[0] << ŝ
181  	for i := 1; i < len(z); i++ {
182  		z[i-1] = x[i-1]>>s | x[i]<<ŝ
183  	}
184  	z[len(z)-1] = x[len(z)-1] >> s
185  	return
186  }
187
188  func mulAddVWW_g(z, x []Word, y, r Word) (c Word) {
189  	c = r
190  	// The comment near the top of this file discusses this for loop condition.
191  	for i := 0; i < len(z) && i < len(x); i++ {
192  		c, z[i] = mulAddWWW_g(x[i], y, c)
193  	}
194  	return
195  }
196
197  func addMulVVW_g(z, x []Word, y Word) (c Word) {
198  	// The comment near the top of this file discusses this for loop condition.
199  	for i := 0; i < len(z) && i < len(x); i++ {
200  		z1, z0 := mulAddWWW_g(x[i], y, z[i])
201  		lo, cc := bits.Add(uint(z0), uint(c), 0)
202  		c, z[i] = Word(cc), Word(lo)
203  		c += z1
204  	}
205  	return
206  }
207
208  // q = ( x1 << _W + x0 - r)/y. m = floor(( _B^2 - 1 ) / d - _B). Requiring x1<y.
209  // An approximate reciprocal with a reference to "Improved Division by Invariant Integers
210  // (IEEE Transactions on Computers, 11 Jun. 2010)"
211  func divWW(x1, x0, y, m Word) (q, r Word) {
212  	s := nlz(y)
213  	if s != 0 {
214  		x1 = x1<<s | x0>>(_W-s)
215  		x0 <<= s
216  		y <<= s
217  	}
218  	d := uint(y)
219  	// We know that
220  	//   m = ⎣(B^2-1)/d⎦-B
221  	//   ⎣(B^2-1)/d⎦ = m+B
222  	//   (B^2-1)/d = m+B+delta1    0 <= delta1 <= (d-1)/d
223  	//   B^2/d = m+B+delta2        0 <= delta2 <= 1
224  	// The quotient we're trying to compute is
225  	//   quotient = ⎣(x1*B+x0)/d⎦
226  	//            = ⎣(x1*B*(B^2/d)+x0*(B^2/d))/B^2⎦
227  	//            = ⎣(x1*B*(m+B+delta2)+x0*(m+B+delta2))/B^2⎦
228  	//            = ⎣(x1*m+x1*B+x0)/B + x0*m/B^2 + delta2*(x1*B+x0)/B^2⎦
229  	// The latter two terms of this three-term sum are between 0 and 1.
230  	// So we can compute just the first term, and we will be low by at most 2.
231  	t1, t0 := bits.Mul(uint(m), uint(x1))
232  	_, c := bits.Add(t0, uint(x0), 0)
233  	t1, _ = bits.Add(t1, uint(x1), c)
234  	// The quotient is either t1, t1+1, or t1+2.
235  	// We'll try t1 and adjust if needed.
236  	qq := t1
237  	// compute remainder r=x-d*q.
238  	dq1, dq0 := bits.Mul(d, qq)
239  	r0, b := bits.Sub(uint(x0), dq0, 0)
240  	r1, _ := bits.Sub(uint(x1), dq1, b)
241  	// The remainder we just computed is bounded above by B+d:
242  	// r = x1*B + x0 - d*q.
243  	//   = x1*B + x0 - d*⎣(x1*m+x1*B+x0)/B⎦
244  	//   = x1*B + x0 - d*((x1*m+x1*B+x0)/B-alpha)                                   0 <= alpha < 1
245  	//   = x1*B + x0 - x1*d/B*m                         - x1*d - x0*d/B + d*alpha
246  	//   = x1*B + x0 - x1*d/B*⎣(B^2-1)/d-B⎦             - x1*d - x0*d/B + d*alpha
247  	//   = x1*B + x0 - x1*d/B*⎣(B^2-1)/d-B⎦             - x1*d - x0*d/B + d*alpha
248  	//   = x1*B + x0 - x1*d/B*((B^2-1)/d-B-beta)        - x1*d - x0*d/B + d*alpha   0 <= beta < 1
249  	//   = x1*B + x0 - x1*B + x1/B + x1*d + x1*d/B*beta - x1*d - x0*d/B + d*alpha
250  	//   =        x0        + x1/B        + x1*d/B*beta        - x0*d/B + d*alpha
251  	//   = x0*(1-d/B) + x1*(1+d*beta)/B + d*alpha
252  	//   <  B*(1-d/B) +  d*B/B          + d          because x0<B (and 1-d/B>0), x1<d, 1+d*beta<=B, alpha<1
253  	//   =  B - d     +  d              + d
254  	//   = B+d
255  	// So r1 can only be 0 or 1. If r1 is 1, then we know q was too small.
256  	// Add 1 to q and subtract d from r. That guarantees that r is <B, so
257  	// we no longer need to keep track of r1.
258  	if r1 != 0 {
259  		qq++
260  		r0 -= d
261  	}
262  	// If the remainder is still too large, increment q one more time.
263  	if r0 >= d {
264  		qq++
265  		r0 -= d
266  	}
267  	return Word(qq), Word(r0 >> s)
268  }
269
270  // reciprocalWord return the reciprocal of the divisor. rec = floor(( _B^2 - 1 ) / u - _B). u = d1 << nlz(d1).
271  func reciprocalWord(d1 Word) Word {
272  	u := uint(d1 << nlz(d1))
273  	x1 := ^u
274  	x0 := uint(_M)
275  	rec, _ := bits.Div(x1, x0, u) // (_B^2-1)/U-_B = (_B*(_M-C)+_M)/U
276  	return Word(rec)
277  }
278
```

View as plain text