Branch data Line data Source code
1 : : // Copyright (c) 2017-2022 The Bitcoin Core developers
2 : : // Distributed under the MIT software license, see the accompanying
3 : : // file COPYING or http://www.opensource.org/licenses/mit-license.php.
4 : :
5 : : #include <crypto/muhash.h>
6 : :
7 : : #include <crypto/chacha20.h>
8 : : #include <crypto/common.h>
9 : : #include <hash.h>
10 : :
11 : : #include <cassert>
12 : : #include <cstdio>
13 : : #include <limits>
14 : :
15 : : namespace {
16 : :
17 : : using limb_t = Num3072::limb_t;
18 : : using double_limb_t = Num3072::double_limb_t;
19 : : constexpr int LIMB_SIZE = Num3072::LIMB_SIZE;
20 : : /** 2^3072 - 1103717, the largest 3072-bit safe prime number, is used as the modulus. */
21 : : constexpr limb_t MAX_PRIME_DIFF = 1103717;
22 : :
23 : : /** Extract the lowest limb of [c0,c1,c2] into n, and left shift the number by 1 limb. */
24 : 1591958640 : inline void extract3(limb_t& c0, limb_t& c1, limb_t& c2, limb_t& n)
25 : : {
26 : 1591958640 : n = c0;
27 : 1591958640 : c0 = c1;
28 : 1591958640 : c1 = c2;
29 : 1591958640 : c2 = 0;
30 : : }
31 : :
32 : : /** [c0,c1] = a * b */
33 : 13665250 : inline void mul(limb_t& c0, limb_t& c1, const limb_t& a, const limb_t& b)
34 : : {
35 : 13665250 : double_limb_t t = (double_limb_t)a * b;
36 : 13665250 : c1 = t >> LIMB_SIZE;
37 : 13665250 : c0 = t;
38 : : }
39 : :
40 : : /* [c0,c1,c2] += n * [d0,d1,d2]. c2 is 0 initially */
41 : 1558792835 : inline void mulnadd3(limb_t& c0, limb_t& c1, limb_t& c2, limb_t& d0, limb_t& d1, limb_t& d2, const limb_t& n)
42 : : {
43 : 1558792835 : double_limb_t t = (double_limb_t)d0 * n + c0;
44 : 1558792835 : c0 = t;
45 : 1558792835 : t >>= LIMB_SIZE;
46 : 1558792835 : t += (double_limb_t)d1 * n + c1;
47 : 1558792835 : c1 = t;
48 : 1558792835 : t >>= LIMB_SIZE;
49 : 1558792835 : c2 = t + d2 * n;
50 : 1558792835 : }
51 : :
52 : : /* [c0,c1] *= n */
53 : 33165805 : inline void muln2(limb_t& c0, limb_t& c1, const limb_t& n)
54 : : {
55 : 33165805 : double_limb_t t = (double_limb_t)c0 * n;
56 : 33165805 : c0 = t;
57 : 33165805 : t >>= LIMB_SIZE;
58 : 33165805 : t += (double_limb_t)c1 * n;
59 : 33165805 : c1 = t;
60 : : }
61 : :
62 : : /** [c0,c1,c2] += a * b */
63 : 2234225390 : inline void muladd3(limb_t& c0, limb_t& c1, limb_t& c2, const limb_t& a, const limb_t& b)
64 : : {
65 : 2234225390 : double_limb_t t = (double_limb_t)a * b;
66 : 2234225390 : limb_t th = t >> LIMB_SIZE;
67 : 2234225390 : limb_t tl = t;
68 : :
69 : 2234225390 : c0 += tl;
70 [ + + ]: 2234225390 : th += (c0 < tl) ? 1 : 0;
71 : 2234225390 : c1 += th;
72 [ + + ]: 2234225390 : c2 += (c1 < th) ? 1 : 0;
73 : 2234225390 : }
74 : :
75 : : /** [c0,c1,c2] += 2 * a * b */
76 : 37083062040 : inline void muldbladd3(limb_t& c0, limb_t& c1, limb_t& c2, const limb_t& a, const limb_t& b)
77 : : {
78 : 37083062040 : double_limb_t t = (double_limb_t)a * b;
79 : 37083062040 : limb_t th = t >> LIMB_SIZE;
80 : 37083062040 : limb_t tl = t;
81 : :
82 : 37083062040 : c0 += tl;
83 [ + + ]: 37083062040 : limb_t tt = th + ((c0 < tl) ? 1 : 0);
84 : 37083062040 : c1 += tt;
85 [ + + ]: 37083062040 : c2 += (c1 < tt) ? 1 : 0;
86 : 37083062040 : c0 += tl;
87 [ + + ]: 37083062040 : th += (c0 < tl) ? 1 : 0;
88 : 37083062040 : c1 += th;
89 [ + + ]: 37083062040 : c2 += (c1 < th) ? 1 : 0;
90 : 37083062040 : }
91 : :
92 : : /**
93 : : * Add limb a to [c0,c1]: [c0,c1] += a. Then extract the lowest
94 : : * limb of [c0,c1] into n, and left shift the number by 1 limb.
95 : : * */
96 : 1591959360 : inline void addnextract2(limb_t& c0, limb_t& c1, const limb_t& a, limb_t& n)
97 : : {
98 : 1591959360 : limb_t c2 = 0;
99 : :
100 : : // add
101 : 1591959360 : c0 += a;
102 [ + + ]: 1591959360 : if (c0 < a) {
103 : 459742 : c1 += 1;
104 : :
105 : : // Handle case when c1 has overflown
106 [ - + ]: 459742 : if (c1 == 0)
107 : 0 : c2 = 1;
108 : : }
109 : :
110 : : // extract
111 : 1591959360 : n = c0;
112 : 1591959360 : c0 = c1;
113 : 1591959360 : c1 = c2;
114 : 1591959360 : }
115 : :
116 : : /** in_out = in_out^(2^sq) * mul */
117 : 149870 : inline void square_n_mul(Num3072& in_out, const int sq, const Num3072& mul)
118 : : {
119 [ + + ]: 11111790 : for (int j = 0; j < sq; ++j) in_out.Square();
120 : 149870 : in_out.Multiply(mul);
121 : 149870 : }
122 : :
123 : : } // namespace
124 : :
125 : : /** Indicates whether d is larger than the modulus. */
126 : 33197920 : bool Num3072::IsOverflow() const
127 : : {
128 [ + + ]: 33197920 : if (this->limbs[0] <= std::numeric_limits<limb_t>::max() - MAX_PRIME_DIFF) return false;
129 [ + + ]: 720 : for (int i = 1; i < LIMBS; ++i) {
130 [ + - ]: 705 : if (this->limbs[i] != std::numeric_limits<limb_t>::max()) return false;
131 : : }
132 : : return true;
133 : : }
134 : :
135 : 15 : void Num3072::FullReduce()
136 : : {
137 : 15 : limb_t c0 = MAX_PRIME_DIFF;
138 : 15 : limb_t c1 = 0;
139 [ + + ]: 735 : for (int i = 0; i < LIMBS; ++i) {
140 : 720 : addnextract2(c0, c1, this->limbs[i], this->limbs[i]);
141 : : }
142 : 15 : }
143 : :
144 : 10705 : Num3072 Num3072::GetInverse() const
145 : : {
146 : : // For fast exponentiation a sliding window exponentiation with repunit
147 : : // precomputation is utilized. See "Fast Point Decompression for Standard
148 : : // Elliptic Curves" (Brumley, Järvinen, 2008).
149 : :
150 [ + + ]: 139165 : Num3072 p[12]; // p[i] = a^(2^(2^i)-1)
151 : 10705 : Num3072 out;
152 : :
153 : 10705 : p[0] = *this;
154 : :
155 [ + + ]: 128460 : for (int i = 0; i < 11; ++i) {
156 : 117755 : p[i + 1] = p[i];
157 [ + + ]: 22030890 : for (int j = 0; j < (1 << i); ++j) p[i + 1].Square();
158 : 117755 : p[i + 1].Multiply(p[i]);
159 : : }
160 : :
161 : 10705 : out = p[11];
162 : :
163 : 10705 : square_n_mul(out, 512, p[9]);
164 : 10705 : square_n_mul(out, 256, p[8]);
165 : 10705 : square_n_mul(out, 128, p[7]);
166 : 10705 : square_n_mul(out, 64, p[6]);
167 : 10705 : square_n_mul(out, 32, p[5]);
168 : 10705 : square_n_mul(out, 8, p[3]);
169 : 10705 : square_n_mul(out, 2, p[1]);
170 : 10705 : square_n_mul(out, 1, p[0]);
171 : 10705 : square_n_mul(out, 5, p[2]);
172 : 10705 : square_n_mul(out, 3, p[0]);
173 : 10705 : square_n_mul(out, 2, p[0]);
174 : 10705 : square_n_mul(out, 4, p[0]);
175 : 10705 : square_n_mul(out, 4, p[1]);
176 : 10705 : square_n_mul(out, 3, p[0]);
177 : :
178 : 10705 : return out;
179 : : }
180 : :
181 : 290750 : void Num3072::Multiply(const Num3072& a)
182 : : {
183 : 290750 : limb_t c0 = 0, c1 = 0, c2 = 0;
184 : 290750 : Num3072 tmp;
185 : :
186 : : /* Compute limbs 0..N-2 of this*a into tmp, including one reduction. */
187 [ + + ]: 13956000 : for (int j = 0; j < LIMBS - 1; ++j) {
188 : 13665250 : limb_t d0 = 0, d1 = 0, d2 = 0;
189 : 13665250 : mul(d0, d1, this->limbs[1 + j], a.limbs[LIMBS + j - (1 + j)]);
190 [ + + ]: 327966000 : for (int i = 2 + j; i < LIMBS; ++i) muladd3(d0, d1, d2, this->limbs[i], a.limbs[LIMBS + j - i]);
191 : 13665250 : mulnadd3(c0, c1, c2, d0, d1, d2, MAX_PRIME_DIFF);
192 [ + + ]: 341631250 : for (int i = 0; i < j + 1; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[j - i]);
193 : 13665250 : extract3(c0, c1, c2, tmp.limbs[j]);
194 : : }
195 : :
196 : : /* Compute limb N-1 of a*b into tmp. */
197 [ + - ]: 290750 : assert(c2 == 0);
198 [ + + ]: 14246750 : for (int i = 0; i < LIMBS; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[LIMBS - 1 - i]);
199 : 290750 : extract3(c0, c1, c2, tmp.limbs[LIMBS - 1]);
200 : :
201 : : /* Perform a second reduction. */
202 : 290750 : muln2(c0, c1, MAX_PRIME_DIFF);
203 [ + + ]: 14246750 : for (int j = 0; j < LIMBS; ++j) {
204 : 13956000 : addnextract2(c0, c1, tmp.limbs[j], this->limbs[j]);
205 : : }
206 : :
207 [ - + ]: 290750 : assert(c1 == 0);
208 [ - + ]: 290750 : assert(c0 == 0 || c0 == 1);
209 : :
210 : : /* Perform up to two more reductions if the internal state has already
211 : : * overflown the MAX of Num3072 or if it is larger than the modulus or
212 : : * if both are the case.
213 : : * */
214 [ + + ]: 290750 : if (this->IsOverflow()) this->FullReduce();
215 [ - + ]: 290750 : if (c0) this->FullReduce();
216 : 290750 : }
217 : :
218 : 32875055 : void Num3072::Square()
219 : : {
220 : 32875055 : limb_t c0 = 0, c1 = 0, c2 = 0;
221 : 32875055 : Num3072 tmp;
222 : :
223 : : /* Compute limbs 0..N-2 of this*this into tmp, including one reduction. */
224 [ + + ]: 1578002640 : for (int j = 0; j < LIMBS - 1; ++j) {
225 : 1545127585 : limb_t d0 = 0, d1 = 0, d2 = 0;
226 [ + + ]: 19692157945 : for (int i = 0; i < (LIMBS - 1 - j) / 2; ++i) muldbladd3(d0, d1, d2, this->limbs[i + j + 1], this->limbs[LIMBS - 1 - i]);
227 [ + + ]: 1545127585 : if ((j + 1) & 1) muladd3(d0, d1, d2, this->limbs[(LIMBS - 1 - j) / 2 + j + 1], this->limbs[LIMBS - 1 - (LIMBS - 1 - j) / 2]);
228 : 1545127585 : mulnadd3(c0, c1, c2, d0, d1, d2, MAX_PRIME_DIFF);
229 [ + + ]: 19692157945 : for (int i = 0; i < (j + 1) / 2; ++i) muldbladd3(c0, c1, c2, this->limbs[i], this->limbs[j - i]);
230 [ + + ]: 1545127585 : if ((j + 1) & 1) muladd3(c0, c1, c2, this->limbs[(j + 1) / 2], this->limbs[j - (j + 1) / 2]);
231 : 1545127585 : extract3(c0, c1, c2, tmp.limbs[j]);
232 : : }
233 : :
234 [ + - ]: 32875055 : assert(c2 == 0);
235 [ + + ]: 821876375 : for (int i = 0; i < LIMBS / 2; ++i) muldbladd3(c0, c1, c2, this->limbs[i], this->limbs[LIMBS - 1 - i]);
236 : 32875055 : extract3(c0, c1, c2, tmp.limbs[LIMBS - 1]);
237 : :
238 : : /* Perform a second reduction. */
239 : 32875055 : muln2(c0, c1, MAX_PRIME_DIFF);
240 [ + + ]: 1610877695 : for (int j = 0; j < LIMBS; ++j) {
241 : 1578002640 : addnextract2(c0, c1, tmp.limbs[j], this->limbs[j]);
242 : : }
243 : :
244 [ - + ]: 32875055 : assert(c1 == 0);
245 [ - + ]: 32875055 : assert(c0 == 0 || c0 == 1);
246 : :
247 : : /* Perform up to two more reductions if the internal state has already
248 : : * overflown the MAX of Num3072 or if it is larger than the modulus or
249 : : * if both are the case.
250 : : * */
251 [ - + ]: 32875055 : if (this->IsOverflow()) this->FullReduce();
252 [ - + ]: 32875055 : if (c0) this->FullReduce();
253 : 32875055 : }
254 : :
255 : 33326972 : void Num3072::SetToOne()
256 : : {
257 : 33326972 : this->limbs[0] = 1;
258 [ + + ]: 1599694656 : for (int i = 1; i < LIMBS; ++i) this->limbs[i] = 0;
259 : 33326972 : }
260 : :
261 : 10705 : void Num3072::Divide(const Num3072& a)
262 : : {
263 [ + + ]: 10705 : if (this->IsOverflow()) this->FullReduce();
264 : :
265 : 10705 : Num3072 inv{};
266 [ - + ]: 10705 : if (a.IsOverflow()) {
267 : 0 : Num3072 b = a;
268 : 0 : b.FullReduce();
269 : 0 : inv = b.GetInverse();
270 : : } else {
271 : 10705 : inv = a.GetInverse();
272 : : }
273 : :
274 : 10705 : this->Multiply(inv);
275 [ - + ]: 10705 : if (this->IsOverflow()) this->FullReduce();
276 : 10705 : }
277 : :
278 : 12200 : Num3072::Num3072(const unsigned char (&data)[BYTE_SIZE]) {
279 [ + + ]: 597800 : for (int i = 0; i < LIMBS; ++i) {
280 : 585600 : if (sizeof(limb_t) == 4) {
281 : : this->limbs[i] = ReadLE32(data + 4 * i);
282 : 585600 : } else if (sizeof(limb_t) == 8) {
283 : 585600 : this->limbs[i] = ReadLE64(data + 8 * i);
284 : : }
285 : : }
286 : 12200 : }
287 : :
288 : 10705 : void Num3072::ToBytes(unsigned char (&out)[BYTE_SIZE]) {
289 [ + + ]: 524545 : for (int i = 0; i < LIMBS; ++i) {
290 : 513840 : if (sizeof(limb_t) == 4) {
291 : : WriteLE32(out + i * 4, this->limbs[i]);
292 : 513840 : } else if (sizeof(limb_t) == 8) {
293 : 513840 : WriteLE64(out + i * 8, this->limbs[i]);
294 : : }
295 : : }
296 : 10705 : }
297 : :
298 : 12200 : Num3072 MuHash3072::ToNum3072(Span<const unsigned char> in) {
299 : 12200 : unsigned char tmp[Num3072::BYTE_SIZE];
300 : :
301 : 12200 : uint256 hashed_in{(HashWriter{} << in).GetSHA256()};
302 : 12200 : static_assert(sizeof(tmp) % ChaCha20Aligned::BLOCKLEN == 0);
303 : 12200 : ChaCha20Aligned{MakeByteSpan(hashed_in)}.Keystream(MakeWritableByteSpan(tmp));
304 : 12200 : Num3072 out{tmp};
305 : :
306 : 12200 : return out;
307 : : }
308 : :
309 : 186 : MuHash3072::MuHash3072(Span<const unsigned char> in) noexcept
310 : : {
311 : 186 : m_numerator = ToNum3072(in);
312 : 186 : }
313 : :
314 : 10705 : void MuHash3072::Finalize(uint256& out) noexcept
315 : : {
316 : 10705 : m_numerator.Divide(m_denominator);
317 : 10705 : m_denominator.SetToOne(); // Needed to keep the MuHash object valid
318 : :
319 : 10705 : unsigned char data[Num3072::BYTE_SIZE];
320 : 10705 : m_numerator.ToBytes(data);
321 : :
322 : 10705 : out = (HashWriter{} << data).GetSHA256();
323 : 10705 : }
324 : :
325 : 92 : MuHash3072& MuHash3072::operator*=(const MuHash3072& mul) noexcept
326 : : {
327 : 92 : m_numerator.Multiply(mul.m_numerator);
328 : 92 : m_denominator.Multiply(mul.m_denominator);
329 : 92 : return *this;
330 : : }
331 : :
332 : 111 : MuHash3072& MuHash3072::operator/=(const MuHash3072& div) noexcept
333 : : {
334 : 111 : m_numerator.Multiply(div.m_denominator);
335 : 111 : m_denominator.Multiply(div.m_numerator);
336 : 111 : return *this;
337 : : }
338 : :
339 : 11766 : MuHash3072& MuHash3072::Insert(Span<const unsigned char> in) noexcept {
340 : 11766 : m_numerator.Multiply(ToNum3072(in));
341 : 11766 : return *this;
342 : : }
343 : :
344 : 248 : MuHash3072& MuHash3072::Remove(Span<const unsigned char> in) noexcept {
345 : 248 : m_denominator.Multiply(ToNum3072(in));
346 : 248 : return *this;
347 : : }
|
