arithmetic.gno
10.16 Kb · 452 lines
1// arithmetic provides arithmetic operations for Uint objects.
2// This includes basic binary operations such as addition, subtraction, multiplication, division, and modulo operations
3// as well as overflow checks, and negation. These functions are essential for numeric
4// calculations using 256-bit unsigned integers.
5package uint256
6
7import (
8 "math/bits"
9)
10
11// Add sets z to the sum x+y and returns z.
12func (z *Uint) Add(x, y *Uint) *Uint {
13 var carry uint64
14 z[0], carry = bits.Add64(x[0], y[0], 0)
15 z[1], carry = bits.Add64(x[1], y[1], carry)
16 z[2], carry = bits.Add64(x[2], y[2], carry)
17 z[3], _ = bits.Add64(x[3], y[3], carry)
18 return z
19}
20
21// AddOverflow sets z to the sum x+y and returns z and true if overflow occurred.
22func (z *Uint) AddOverflow(x, y *Uint) (*Uint, bool) {
23 var carry uint64
24 z[0], carry = bits.Add64(x[0], y[0], 0)
25 z[1], carry = bits.Add64(x[1], y[1], carry)
26 z[2], carry = bits.Add64(x[2], y[2], carry)
27 z[3], carry = bits.Add64(x[3], y[3], carry)
28 return z, carry != 0
29}
30
31// Sub sets z to the difference x-y and returns z.
32func (z *Uint) Sub(x, y *Uint) *Uint {
33 var carry uint64
34 z[0], carry = bits.Sub64(x[0], y[0], 0)
35 z[1], carry = bits.Sub64(x[1], y[1], carry)
36 z[2], carry = bits.Sub64(x[2], y[2], carry)
37 z[3], _ = bits.Sub64(x[3], y[3], carry)
38 return z
39}
40
41// SubOverflow sets z to the difference x-y and returns z and true if underflow occurred.
42func (z *Uint) SubOverflow(x, y *Uint) (*Uint, bool) {
43 var carry uint64
44 z[0], carry = bits.Sub64(x[0], y[0], 0)
45 z[1], carry = bits.Sub64(x[1], y[1], carry)
46 z[2], carry = bits.Sub64(x[2], y[2], carry)
47 z[3], carry = bits.Sub64(x[3], y[3], carry)
48 return z, carry != 0
49}
50
51// Neg returns -x mod 2^256.
52func (z *Uint) Neg(x *Uint) *Uint {
53 return z.Sub(Zero(), x)
54}
55
56// Mul sets z to the product x*y and returns z.
57func (z *Uint) Mul(x, y *Uint) *Uint {
58 var (
59 res Uint
60 carry uint64
61 res1, res2, res3 uint64
62 )
63
64 carry, res[0] = bits.Mul64(x[0], y[0])
65 carry, res1 = umulHop(carry, x[1], y[0])
66 carry, res2 = umulHop(carry, x[2], y[0])
67 res3 = x[3]*y[0] + carry
68
69 carry, res[1] = umulHop(res1, x[0], y[1])
70 carry, res2 = umulStep(res2, x[1], y[1], carry)
71 res3 = res3 + x[2]*y[1] + carry
72
73 carry, res[2] = umulHop(res2, x[0], y[2])
74 res3 = res3 + x[1]*y[2] + carry
75
76 res[3] = res3 + x[0]*y[3]
77
78 return z.Set(&res)
79}
80
81// MulOverflow sets z to the product x*y and returns z and true if overflow occurred.
82func (z *Uint) MulOverflow(x, y *Uint) (*Uint, bool) {
83 p := umul(x, y)
84 copy(z[:], p[:4])
85 return z, (p[4] | p[5] | p[6] | p[7]) != 0
86}
87
88// Div sets z to the quotient x/y and returns z.
89// It panics if y == 0.
90func (z *Uint) Div(x, y *Uint) *Uint {
91 if y.IsZero() {
92 panic("division by zero")
93 }
94 if y.Gt(x) {
95 return z.Clear()
96 }
97 if x.Eq(y) {
98 return z.SetOne()
99 }
100 // Shortcut some cases
101 if x.IsUint64() {
102 return z.SetUint64(x.Uint64() / y.Uint64())
103 }
104
105 // At this point, we know
106 // x/y ; x > y > 0
107
108 var quot Uint
109 udivrem(quot[:], x[:], y)
110 return z.Set(")
111}
112
113// Mod sets z to the modulus x%y and returns z.
114// It panics if y == 0.
115func (z *Uint) Mod(x, y *Uint) *Uint {
116 if y.IsZero() {
117 panic("modulo by zero")
118 }
119 if x.IsZero() {
120 return z.Clear()
121 }
122 switch x.Cmp(y) {
123 case -1:
124 // x < y
125 copy(z[:], x[:])
126 return z
127 case 0:
128 // x == y
129 return z.Clear() // They are equal
130 }
131
132 // At this point:
133 // x != 0
134 // y != 0
135 // x > y
136
137 // Shortcut trivial case
138 if x.IsUint64() {
139 return z.SetUint64(x.Uint64() % y.Uint64())
140 }
141
142 var quot Uint
143 *z = udivrem(quot[:], x[:], y)
144 return z
145}
146
147// MulMod sets z to (x * y) mod m and returns z.
148// It panics if m == 0.
149func (z *Uint) MulMod(x, y, m *Uint) *Uint {
150 if m.IsZero() {
151 panic("modulo by zero")
152 }
153 if x.IsZero() || y.IsZero() {
154 return z.Clear()
155 }
156 p := umul(x, y)
157
158 if m[3] != 0 {
159 mu := Reciprocal(m)
160 r := reduce4(p, m, mu)
161 return z.Set(&r)
162 }
163
164 var (
165 pl Uint
166 ph Uint
167 )
168
169 pl[0], pl[1], pl[2], pl[3] = p[0], p[1], p[2], p[3]
170 ph[0], ph[1], ph[2], ph[3] = p[4], p[5], p[6], p[7]
171
172 // If the multiplication is within 256 bits use Mod().
173 if ph.IsZero() {
174 return z.Mod(&pl, m)
175 }
176
177 var quot [8]uint64
178 rem := udivrem(quot[:], p[:], m)
179 return z.Set(&rem)
180}
181
182// DivMod sets z to the quotient x/y and m to the modulus x%y, returning the pair (z, m).
183// It panics if y == 0.
184func (z *Uint) DivMod(x, y, m *Uint) (*Uint, *Uint) {
185 if y.IsZero() {
186 panic("division by zero")
187 }
188
189 switch x.Cmp(y) {
190 case -1:
191 // x < y
192 return z.Clear(), m.Set(x)
193 case 0:
194 // x == y
195 return z.SetOne(), m.Clear()
196 }
197
198 // At this point:
199 // x != 0
200 // y != 0
201 // x > y
202
203 // Shortcut trivial case
204 if x.IsUint64() {
205 x0, y0 := x.Uint64(), y.Uint64()
206 return z.SetUint64(x0 / y0), m.SetUint64(x0 % y0)
207 }
208
209 var quot Uint
210 *m = udivrem(quot[:], x[:], y)
211 *z = quot
212 return z, m
213}
214
215// udivrem divides u by d and produces both quotient and remainder.
216// The quotient is stored in provided quot - len(u)-len(d)+1 words.
217// It loosely follows the Knuth's division algorithm (sometimes referenced as "schoolbook" division) using 64-bit words.
218// See Knuth, Volume 2, section 4.3.1, Algorithm D.
219func udivrem(quot, u []uint64, d *Uint) (rem Uint) {
220 var dLen int
221 for i := len(d) - 1; i >= 0; i-- {
222 if d[i] != 0 {
223 dLen = i + 1
224 break
225 }
226 }
227
228 shift := uint(bits.LeadingZeros64(d[dLen-1]))
229
230 var dnStorage Uint
231 dn := dnStorage[:dLen]
232 for i := dLen - 1; i > 0; i-- {
233 dn[i] = (d[i] << shift) | (d[i-1] >> (64 - shift))
234 }
235 dn[0] = d[0] << shift
236
237 var uLen int
238 for i := len(u) - 1; i >= 0; i-- {
239 if u[i] != 0 {
240 uLen = i + 1
241 break
242 }
243 }
244
245 if uLen < dLen {
246 copy(rem[:], u)
247 return rem
248 }
249
250 var unStorage [9]uint64
251 un := unStorage[:uLen+1]
252 un[uLen] = u[uLen-1] >> (64 - shift)
253 for i := uLen - 1; i > 0; i-- {
254 un[i] = (u[i] << shift) | (u[i-1] >> (64 - shift))
255 }
256 un[0] = u[0] << shift
257
258 if dLen == 1 {
259 r := udivremBy1(quot, un, dn[0])
260 rem.SetUint64(r >> shift)
261 return rem
262 }
263
264 udivremKnuth(quot, un, dn)
265
266 for i := 0; i < dLen-1; i++ {
267 rem[i] = (un[i] >> shift) | (un[i+1] << (64 - shift))
268 }
269 rem[dLen-1] = un[dLen-1] >> shift
270
271 return rem
272}
273
274// umul computes full 256 x 256 -> 512 multiplication.
275func umul(x, y *Uint) [8]uint64 {
276 var res [8]uint64
277
278 topX := highestNonZeroWord(x)
279 topY := highestNonZeroWord(y)
280
281 if topX < 0 || topY < 0 {
282 return res
283 }
284
285 lenX := topX + 1
286 lenY := topY + 1
287
288 for i := 0; i < lenX; i++ {
289 xi := x[i]
290 if xi == 0 {
291 continue
292 }
293 var carry uint64
294 k := i
295 for j := 0; j < lenY; j++ {
296 hi, lo := bits.Mul64(xi, y[j])
297 lo, c := bits.Add64(lo, res[k], 0)
298 hi += c
299 lo, c = bits.Add64(lo, carry, 0)
300 hi += c
301 res[k] = lo
302 carry = hi
303 k++
304 }
305 res[i+lenY] = carry
306 }
307
308 return res
309}
310
311// highestNonZeroWord returns the highest index with non-zero value or -1 if the Uint is zero.
312func highestNonZeroWord(u *Uint) int {
313 for i := 3; i >= 0; i-- {
314 if u[i] != 0 {
315 return i
316 }
317 }
318 return -1
319}
320
321// umulStep computes (hi * 2^64 + lo) = z + (x * y) + carry.
322func umulStep(z, x, y, carry uint64) (hi, lo uint64) {
323 hi, lo = bits.Mul64(x, y)
324 lo, carry = bits.Add64(lo, carry, 0)
325 hi += carry
326 lo, carry = bits.Add64(lo, z, 0)
327 hi += carry
328 return hi, lo
329}
330
331// umulHop computes (hi * 2^64 + lo) = z + (x * y)
332func umulHop(z, x, y uint64) (hi, lo uint64) {
333 hi, lo = bits.Mul64(x, y)
334 lo, carry := bits.Add64(lo, z, 0)
335 hi += carry
336 return hi, lo
337}
338
339// udivremBy1 divides u by single normalized word d and produces both quotient and remainder.
340// The quotient is stored in provided quot.
341func udivremBy1(quot, u []uint64, d uint64) (rem uint64) {
342 reciprocal := reciprocal2by1(d)
343 rem = u[len(u)-1] // Set the top word as remainder.
344 for j := len(u) - 2; j >= 0; j-- {
345 quot[j], rem = udivrem2by1(rem, u[j], d, reciprocal)
346 }
347 return rem
348}
349
350// udivremKnuth implements the division of u by normalized multiple word d from the Knuth's division algorithm.
351// The quotient is stored in provided quot - len(u)-len(d) words.
352// Updates u to contain the remainder - len(d) words.
353func udivremKnuth(quot, u, d []uint64) {
354 dLen := len(d)
355 dh := d[dLen-1]
356 dl := d[dLen-2]
357 reciprocal := reciprocal2by1(dh)
358
359 for j := len(u) - dLen - 1; j >= 0; j-- {
360 u2 := u[j+dLen]
361 u1 := u[j+dLen-1]
362 u0 := u[j+dLen-2]
363
364 var qhat, rhat uint64
365 if u2 >= dh { // Division overflows.
366 qhat = 18446744073709551615 // max uint64
367 // NOTE: Add "qhat one to big" adjustment (not needed for correctness, but helps avoiding "add back" case).
368 } else {
369 qhat, rhat = udivrem2by1(u2, u1, dh, reciprocal)
370 ph, pl := bits.Mul64(qhat, dl)
371 if ph > rhat || (ph == rhat && pl > u0) {
372 qhat--
373 // NOTE: Add "qhat one to big" adjustment (not needed for correctness, but helps avoiding "add back" case).
374 }
375 }
376
377 // Multiply and subtract.
378 borrow := subMulTo(u[j:], d, qhat)
379 u[j+dLen] = u2 - borrow
380 if u2 < borrow { // Too much subtracted, add back.
381 qhat--
382 u[j+dLen] += addTo(u[j:], d)
383 }
384
385 quot[j] = qhat // Store quotient digit.
386 }
387}
388
389// isBitSet returns true if bit n-th is set, where n = 0 is LSB.
390// The n must be <= 255.
391func (z *Uint) isBitSet(n uint) bool {
392 return (z[n/64] & (1 << (n % 64))) != 0
393}
394
395func (z *Uint) IsOverflow() bool {
396 return z.isBitSet(255)
397}
398
399// addTo computes x += y.
400// Requires len(x) >= len(y).
401func addTo(x, y []uint64) uint64 {
402 var carry uint64
403 for i := 0; i < len(y); i++ {
404 x[i], carry = bits.Add64(x[i], y[i], carry)
405 }
406 return carry
407}
408
409// subMulTo computes x -= y * multiplier.
410// Requires len(x) >= len(y).
411func subMulTo(x, y []uint64, multiplier uint64) uint64 {
412 var borrow uint64
413 for i := 0; i < len(y); i++ {
414 s, carry1 := bits.Sub64(x[i], borrow, 0)
415 ph, pl := bits.Mul64(y[i], multiplier)
416 t, carry2 := bits.Sub64(s, pl, 0)
417 x[i] = t
418 borrow = ph + carry1 + carry2
419 }
420 return borrow
421}
422
423// reciprocal2by1 computes <^d, ^0> / d.
424func reciprocal2by1(d uint64) uint64 {
425 reciprocal, _ := bits.Div64(^d, 18446744073709551615, d)
426 return reciprocal
427}
428
429// udivrem2by1 divides <uh, ul> / d and produces both quotient and remainder.
430// It uses the provided d's reciprocal.
431// Implementation ported from https://github.com/chfast/intx and is based on
432// "Improved division by invariant integers", Algorithm 4.
433func udivrem2by1(uh, ul, d, reciprocal uint64) (quot, rem uint64) {
434 qh, ql := bits.Mul64(reciprocal, uh)
435 ql, carry := bits.Add64(ql, ul, 0)
436 qh, _ = bits.Add64(qh, uh, carry)
437 qh++
438
439 r := ul - qh*d
440
441 if r > ql {
442 qh--
443 r += d
444 }
445
446 if r >= d {
447 qh++
448 r -= d
449 }
450
451 return qh, r
452}
