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