Source file src/math/big/arith.go
1 // Copyright 2009 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 // 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(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 { 60 return uint(bits.LeadingZeros(uint(x))) 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 97 // addVWlarge is addVW, but intended for large z. 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 130 // subVWlarge is to subVW as addVWlarge is to addVW. 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