arithmetic.gno
2.59 Kb ยท 146 lines
1package int256
2
3import "math/bits"
4
5func udivrem(quot, u []uint64, d *Int) (rem Int) {
6 var dLen int
7 for i := len(d) - 1; i >= 0; i-- {
8 if d[i] != 0 {
9 dLen = i + 1
10 break
11 }
12 }
13
14 shift := uint(bits.LeadingZeros64(d[dLen-1]))
15
16 var dnStorage Int
17 dn := dnStorage[:dLen]
18 for i := dLen - 1; i > 0; i-- {
19 dn[i] = (d[i] << shift) | (d[i-1] >> (64 - shift))
20 }
21 dn[0] = d[0] << shift
22
23 var uLen int
24 for i := len(u) - 1; i >= 0; i-- {
25 if u[i] != 0 {
26 uLen = i + 1
27 break
28 }
29 }
30
31 if uLen < dLen {
32 copy(rem[:], u)
33 return rem
34 }
35
36 var unStorage [9]uint64
37 un := unStorage[:uLen+1]
38 un[uLen] = u[uLen-1] >> (64 - shift)
39 for i := uLen - 1; i > 0; i-- {
40 un[i] = (u[i] << shift) | (u[i-1] >> (64 - shift))
41 }
42 un[0] = u[0] << shift
43
44 if dLen == 1 {
45 r := udivremBy1(quot, un, dn[0])
46 rem.SetUint64(r >> shift)
47 return rem
48 }
49
50 udivremKnuth(quot, un, dn)
51
52 for i := 0; i < dLen-1; i++ {
53 rem[i] = (un[i] >> shift) | (un[i+1] << (64 - shift))
54 }
55 rem[dLen-1] = un[dLen-1] >> shift
56
57 return rem
58}
59
60func udivremBy1(quot, u []uint64, d uint64) (rem uint64) {
61 reciprocal := reciprocal2by1(d)
62 rem = u[len(u)-1]
63 for j := len(u) - 2; j >= 0; j-- {
64 quot[j], rem = udivrem2by1(rem, u[j], d, reciprocal)
65 }
66 return rem
67}
68
69func udivremKnuth(quot, u, d []uint64) {
70 dLen := len(d)
71 dh := d[dLen-1]
72 dl := d[dLen-2]
73 reciprocal := reciprocal2by1(dh)
74
75 for j := len(u) - dLen - 1; j >= 0; j-- {
76 u2 := u[j+dLen]
77 u1 := u[j+dLen-1]
78 u0 := u[j+dLen-2]
79
80 var qhat, rhat uint64
81 if u2 >= dh {
82 qhat = MAX_UINT64
83 } else {
84 qhat, rhat = udivrem2by1(u2, u1, dh, reciprocal)
85 ph, pl := bits.Mul64(qhat, dl)
86 if ph > rhat || (ph == rhat && pl > u0) {
87 qhat--
88 }
89 }
90
91 borrow := subMulTo(u[j:], d, qhat)
92 u[j+dLen] = u2 - borrow
93 if u2 < borrow {
94 qhat--
95 u[j+dLen] += addTo(u[j:], d)
96 }
97
98 quot[j] = qhat
99 }
100}
101
102func reciprocal2by1(d uint64) uint64 {
103 reciprocal, _ := bits.Div64(^d, MAX_UINT64, d)
104 return reciprocal
105}
106
107func udivrem2by1(uh, ul, d, reciprocal uint64) (quot, rem uint64) {
108 qh, ql := bits.Mul64(reciprocal, uh)
109 ql, carry := bits.Add64(ql, ul, 0)
110 qh += uh + carry
111 qh++
112
113 r := ul - qh*d
114
115 if r > ql {
116 qh--
117 r += d
118 }
119
120 if r >= d {
121 qh++
122 r -= d
123 }
124
125 return qh, r
126}
127
128func subMulTo(x, y []uint64, multiplier uint64) uint64 {
129 var borrow uint64
130 for i := 0; i < len(y); i++ {
131 s, carry1 := bits.Sub64(x[i], borrow, 0)
132 ph, pl := bits.Mul64(y[i], multiplier)
133 t, carry2 := bits.Sub64(s, pl, 0)
134 x[i] = t
135 borrow = ph + carry1 + carry2
136 }
137 return borrow
138}
139
140func addTo(x, y []uint64) uint64 {
141 var carry uint64
142 for i := 0; i < len(y); i++ {
143 x[i], carry = bits.Add64(x[i], y[i], carry)
144 }
145 return carry
146}