LCOV - code coverage report
Current view: top level - src/crypto - muhash.cpp (source / functions) Coverage Total Hit
Test: test_bitcoin_coverage.info Lines: 99.0 % 298 295
Test Date: 2025-02-16 04:47:53 Functions: 100.0 % 25 25
Branches: 89.0 % 82 73

             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                 :             : #include <util/check.h>
      11                 :             : 
      12                 :             : #include <bit>
      13                 :             : #include <cassert>
      14                 :             : #include <cstdio>
      15                 :             : #include <limits>
      16                 :             : 
      17                 :             : namespace {
      18                 :             : 
      19                 :             : using limb_t = Num3072::limb_t;
      20                 :             : using signed_limb_t = Num3072::signed_limb_t;
      21                 :             : using double_limb_t = Num3072::double_limb_t;
      22                 :             : using signed_double_limb_t = Num3072::signed_double_limb_t;
      23                 :             : constexpr int LIMB_SIZE = Num3072::LIMB_SIZE;
      24                 :             : constexpr int SIGNED_LIMB_SIZE = Num3072::SIGNED_LIMB_SIZE;
      25                 :             : constexpr int LIMBS = Num3072::LIMBS;
      26                 :             : constexpr int SIGNED_LIMBS = Num3072::SIGNED_LIMBS;
      27                 :             : constexpr int FINAL_LIMB_POSITION = 3072 / SIGNED_LIMB_SIZE;
      28                 :             : constexpr int FINAL_LIMB_MODULUS_BITS = 3072 % SIGNED_LIMB_SIZE;
      29                 :             : constexpr limb_t MAX_LIMB = (limb_t)(-1);
      30                 :             : constexpr limb_t MAX_SIGNED_LIMB = (((limb_t)1) << SIGNED_LIMB_SIZE) - 1;
      31                 :             : /** 2^3072 - 1103717, the largest 3072-bit safe prime number, is used as the modulus. */
      32                 :             : constexpr limb_t MAX_PRIME_DIFF = 1103717;
      33                 :             : /** The modular inverse of (2**3072 - MAX_PRIME_DIFF) mod (MAX_SIGNED_LIMB + 1). */
      34                 :             : constexpr limb_t MODULUS_INVERSE = limb_t(0x70a1421da087d93);
      35                 :             : 
      36                 :             : 
      37                 :             : /** Extract the lowest limb of [c0,c1,c2] into n, and left shift the number by 1 limb. */
      38                 :       42240 : inline void extract3(limb_t& c0, limb_t& c1, limb_t& c2, limb_t& n)
      39                 :             : {
      40                 :       42240 :     n = c0;
      41                 :       42240 :     c0 = c1;
      42                 :       42240 :     c1 = c2;
      43                 :       42240 :     c2 = 0;
      44                 :             : }
      45                 :             : 
      46                 :             : /** [c0,c1] = a * b */
      47                 :       41360 : inline void mul(limb_t& c0, limb_t& c1, const limb_t& a, const limb_t& b)
      48                 :             : {
      49                 :       41360 :     double_limb_t t = (double_limb_t)a * b;
      50                 :       41360 :     c1 = t >> LIMB_SIZE;
      51                 :       41360 :     c0 = t;
      52                 :             : }
      53                 :             : 
      54                 :             : /* [c0,c1,c2] += n * [d0,d1,d2]. c2 is 0 initially */
      55                 :       41360 : 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)
      56                 :             : {
      57                 :       41360 :     double_limb_t t = (double_limb_t)d0 * n + c0;
      58                 :       41360 :     c0 = t;
      59                 :       41360 :     t >>= LIMB_SIZE;
      60                 :       41360 :     t += (double_limb_t)d1 * n + c1;
      61                 :       41360 :     c1 = t;
      62                 :       41360 :     t >>= LIMB_SIZE;
      63                 :       41360 :     c2 = t + d2 * n;
      64                 :       41360 : }
      65                 :             : 
      66                 :             : /* [c0,c1] *= n */
      67                 :         880 : inline void muln2(limb_t& c0, limb_t& c1, const limb_t& n)
      68                 :             : {
      69                 :         880 :     double_limb_t t = (double_limb_t)c0 * n;
      70                 :         880 :     c0 = t;
      71                 :         880 :     t >>= LIMB_SIZE;
      72                 :         880 :     t += (double_limb_t)c1 * n;
      73                 :         880 :     c1 = t;
      74                 :             : }
      75                 :             : 
      76                 :             : /** [c0,c1,c2] += a * b */
      77                 :     1986160 : inline void muladd3(limb_t& c0, limb_t& c1, limb_t& c2, const limb_t& a, const limb_t& b)
      78                 :             : {
      79                 :     1986160 :     double_limb_t t = (double_limb_t)a * b;
      80                 :     1986160 :     limb_t th = t >> LIMB_SIZE;
      81                 :     1986160 :     limb_t tl = t;
      82                 :             : 
      83                 :     1986160 :     c0 += tl;
      84         [ +  + ]:     1986160 :     th += (c0 < tl) ? 1 : 0;
      85                 :     1986160 :     c1 += th;
      86         [ +  + ]:     1986160 :     c2 += (c1 < th) ? 1 : 0;
      87                 :     1986160 : }
      88                 :             : 
      89                 :             : /**
      90                 :             :  * Add limb a to [c0,c1]: [c0,c1] += a. Then extract the lowest
      91                 :             :  * limb of [c0,c1] into n, and left shift the number by 1 limb.
      92                 :             :  * */
      93                 :       42768 : inline void addnextract2(limb_t& c0, limb_t& c1, const limb_t& a, limb_t& n)
      94                 :             : {
      95                 :       42768 :     limb_t c2 = 0;
      96                 :             : 
      97                 :             :     // add
      98                 :       42768 :     c0 += a;
      99         [ +  + ]:       42768 :     if (c0 < a) {
     100                 :         655 :         c1 += 1;
     101                 :             : 
     102                 :             :         // Handle case when c1 has overflown
     103         [ -  + ]:         655 :         if (c1 == 0) c2 = 1;
     104                 :             :     }
     105                 :             : 
     106                 :             :     // extract
     107                 :       42768 :     n = c0;
     108                 :       42768 :     c0 = c1;
     109                 :       42768 :     c1 = c2;
     110                 :       42768 : }
     111                 :             : 
     112                 :             : } // namespace
     113                 :             : 
     114                 :             : /** Indicates whether d is larger than the modulus. */
     115                 :        1690 : bool Num3072::IsOverflow() const
     116                 :             : {
     117         [ +  + ]:        1690 :     if (this->limbs[0] <= std::numeric_limits<limb_t>::max() - MAX_PRIME_DIFF) return false;
     118         [ +  + ]:         528 :     for (int i = 1; i < LIMBS; ++i) {
     119         [ +  - ]:         517 :         if (this->limbs[i] != std::numeric_limits<limb_t>::max()) return false;
     120                 :             :     }
     121                 :             :     return true;
     122                 :             : }
     123                 :             : 
     124                 :          11 : void Num3072::FullReduce()
     125                 :             : {
     126                 :          11 :     limb_t c0 = MAX_PRIME_DIFF;
     127                 :          11 :     limb_t c1 = 0;
     128         [ +  + ]:         539 :     for (int i = 0; i < LIMBS; ++i) {
     129                 :         528 :         addnextract2(c0, c1, this->limbs[i], this->limbs[i]);
     130                 :             :     }
     131                 :          11 : }
     132                 :             : 
     133                 :             : namespace {
     134                 :             : /** A type representing a number in signed limb representation. */
     135                 :             : struct Num3072Signed
     136                 :             : {
     137                 :             :     /** The represented value is sum(limbs[i]*2^(SIGNED_LIMB_SIZE*i), i=0..SIGNED_LIMBS-1).
     138                 :             :      *  Note that limbs may be negative, or exceed 2^SIGNED_LIMB_SIZE-1. */
     139                 :             :     signed_limb_t limbs[SIGNED_LIMBS];
     140                 :             : 
     141                 :             :     /** Construct a Num3072Signed with value 0. */
     142                 :        1080 :     Num3072Signed()
     143                 :        1080 :     {
     144                 :        1080 :         memset(limbs, 0, sizeof(limbs));
     145                 :        1080 :     }
     146                 :             : 
     147                 :             :     /** Convert a Num3072 to a Num3072Signed. Output will be normalized and in
     148                 :             :      *  range 0..2^3072-1. */
     149                 :         270 :     void FromNum3072(const Num3072& in)
     150                 :             :     {
     151                 :         270 :         double_limb_t c = 0;
     152                 :         270 :         int b = 0, outpos = 0;
     153         [ +  + ]:       13230 :         for (int i = 0; i < LIMBS; ++i) {
     154                 :       12960 :             c += double_limb_t{in.limbs[i]} << b;
     155                 :       12960 :             b += LIMB_SIZE;
     156         [ +  + ]:       26190 :             while (b >= SIGNED_LIMB_SIZE) {
     157                 :       13230 :                 limbs[outpos++] = limb_t(c) & MAX_SIGNED_LIMB;
     158                 :       13230 :                 c >>= SIGNED_LIMB_SIZE;
     159                 :       13230 :                 b -= SIGNED_LIMB_SIZE;
     160                 :             :             }
     161                 :             :         }
     162                 :         270 :         Assume(outpos == SIGNED_LIMBS - 1);
     163                 :         270 :         limbs[SIGNED_LIMBS - 1] = c;
     164                 :         270 :         c >>= SIGNED_LIMB_SIZE;
     165                 :         270 :         Assume(c == 0);
     166                 :         270 :     }
     167                 :             : 
     168                 :             :     /** Convert a Num3072Signed to a Num3072. Input must be in range 0..modulus-1. */
     169                 :         270 :     void ToNum3072(Num3072& out) const
     170                 :             :     {
     171                 :         270 :         double_limb_t c = 0;
     172                 :         270 :         int b = 0, outpos = 0;
     173         [ +  + ]:       13770 :         for (int i = 0; i < SIGNED_LIMBS; ++i) {
     174                 :       13500 :             c += double_limb_t(limbs[i]) << b;
     175                 :       13500 :             b += SIGNED_LIMB_SIZE;
     176         [ +  + ]:       13500 :             if (b >= LIMB_SIZE) {
     177                 :       12960 :                 out.limbs[outpos++] = c;
     178                 :       12960 :                 c >>= LIMB_SIZE;
     179                 :       12960 :                 b -= LIMB_SIZE;
     180                 :             :             }
     181                 :             :         }
     182                 :         270 :         Assume(outpos == LIMBS);
     183                 :         270 :         Assume(c == 0);
     184                 :         270 :     }
     185                 :             : 
     186                 :             :     /** Take a Num3072Signed in range 1-2*2^3072..2^3072-1, and:
     187                 :             :      *  - optionally negate it (if negate is true)
     188                 :             :      *  - reduce it modulo the modulus (2^3072 - MAX_PRIME_DIFF)
     189                 :             :      *  - produce output with all limbs in range 0..2^SIGNED_LIMB_SIZE-1
     190                 :             :      */
     191                 :         270 :     void Normalize(bool negate)
     192                 :             :     {
     193                 :             :         // Add modulus if this was negative. This brings the range of *this to 1-2^3072..2^3072-1.
     194                 :         270 :         signed_limb_t cond_add = limbs[SIGNED_LIMBS-1] >> (LIMB_SIZE-1); // -1 if this is negative; 0 otherwise
     195                 :         270 :         limbs[0] += signed_limb_t(-MAX_PRIME_DIFF) & cond_add;
     196                 :         270 :         limbs[FINAL_LIMB_POSITION] += (signed_limb_t(1) << FINAL_LIMB_MODULUS_BITS) & cond_add;
     197                 :             :         // Next negate all limbs if negate was set. This does not change the range of *this.
     198                 :         270 :         signed_limb_t cond_negate = -signed_limb_t(negate); // -1 if this negate is true; 0 otherwise
     199         [ +  + ]:       13770 :         for (int i = 0; i < SIGNED_LIMBS; ++i) {
     200                 :       13500 :             limbs[i] = (limbs[i] ^ cond_negate) - cond_negate;
     201                 :             :         }
     202                 :             :         // Perform carry (make all limbs except the top one be in range 0..2^SIGNED_LIMB_SIZE-1).
     203         [ +  + ]:       13500 :         for (int i = 0; i < SIGNED_LIMBS - 1; ++i) {
     204                 :       13230 :             limbs[i + 1] += limbs[i] >> SIGNED_LIMB_SIZE;
     205                 :       13230 :             limbs[i] &= MAX_SIGNED_LIMB;
     206                 :             :         }
     207                 :             :         // Again add modulus if *this was negative. This brings the range of *this to 0..2^3072-1.
     208                 :         270 :         cond_add = limbs[SIGNED_LIMBS-1] >> (LIMB_SIZE-1); // -1 if this is negative; 0 otherwise
     209                 :         270 :         limbs[0] += signed_limb_t(-MAX_PRIME_DIFF) & cond_add;
     210                 :         270 :         limbs[FINAL_LIMB_POSITION] += (signed_limb_t(1) << FINAL_LIMB_MODULUS_BITS) & cond_add;
     211                 :             :         // Perform another carry. Now all limbs are in range 0..2^SIGNED_LIMB_SIZE-1.
     212         [ +  + ]:       13500 :         for (int i = 0; i < SIGNED_LIMBS - 1; ++i) {
     213                 :       13230 :             limbs[i + 1] += limbs[i] >> SIGNED_LIMB_SIZE;
     214                 :       13230 :             limbs[i] &= MAX_SIGNED_LIMB;
     215                 :             :         }
     216                 :         270 :     }
     217                 :             : };
     218                 :             : 
     219                 :             : /** 2x2 transformation matrix with signed_limb_t elements. */
     220                 :             : struct SignedMatrix
     221                 :             : {
     222                 :             :     signed_limb_t u, v, q, r;
     223                 :             : };
     224                 :             : 
     225                 :             : /** Compute the transformation matrix for SIGNED_LIMB_SIZE divsteps.
     226                 :             :  *
     227                 :             :  * eta: initial eta value
     228                 :             :  * f:   bottom SIGNED_LIMB_SIZE bits of initial f value
     229                 :             :  * g:   bottom SIGNED_LIMB_SIZE bits of initial g value
     230                 :             :  * out: resulting transformation matrix, scaled by 2^SIGNED_LIMB_SIZE
     231                 :             :  * return: eta value after SIGNED_LIMB_SIZE divsteps
     232                 :             :  */
     233                 :       27169 : inline limb_t ComputeDivstepMatrix(signed_limb_t eta, limb_t f, limb_t g, SignedMatrix& out)
     234                 :             : {
     235                 :             :     /** inv256[i] = -1/(2*i+1) (mod 256) */
     236                 :       27169 :     static const uint8_t NEGINV256[128] = {
     237                 :             :         0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59,
     238                 :             :         0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31,
     239                 :             :         0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89,
     240                 :             :         0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61,
     241                 :             :         0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9,
     242                 :             :         0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91,
     243                 :             :         0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9,
     244                 :             :         0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1,
     245                 :             :         0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19,
     246                 :             :         0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1,
     247                 :             :         0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01
     248                 :             :     };
     249                 :             :     // Coefficients of returned SignedMatrix; starts off as identity matrix. */
     250                 :       27169 :     limb_t u = 1, v = 0, q = 0, r = 1;
     251                 :             :     // The number of divsteps still left.
     252                 :       27169 :     int i = SIGNED_LIMB_SIZE;
     253                 :      367699 :     while (true) {
     254                 :             :         /* Use a sentinel bit to count zeros only up to i. */
     255         [ +  - ]:      197434 :         int zeros = std::countr_zero(g | (MAX_LIMB << i));
     256                 :             :         /* Perform zeros divsteps at once; they all just divide g by two. */
     257                 :      197434 :         g >>= zeros;
     258                 :      197434 :         u <<= zeros;
     259                 :      197434 :         v <<= zeros;
     260                 :      197434 :         eta -= zeros;
     261                 :      197434 :         i -= zeros;
     262                 :             :          /* We're done once we've performed SIGNED_LIMB_SIZE divsteps. */
     263         [ +  + ]:      197434 :         if (i == 0) break;
     264                 :             :         /* If eta is negative, negate it and replace f,g with g,-f. */
     265         [ +  + ]:      170265 :         if (eta < 0) {
     266                 :       81536 :             limb_t tmp;
     267                 :       81536 :             eta = -eta;
     268                 :       81536 :             tmp = f; f = g; g = -tmp;
     269                 :       81536 :             tmp = u; u = q; q = -tmp;
     270                 :       81536 :             tmp = v; v = r; r = -tmp;
     271                 :             :         }
     272                 :             :         /* eta is now >= 0. In what follows we're going to cancel out the bottom bits of g. No more
     273                 :             :          * than i can be cancelled out (as we'd be done before that point), and no more than eta+1
     274                 :             :          * can be done as its sign will flip once that happens. */
     275         [ +  + ]:      170265 :         int limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
     276                 :             :         /* m is a mask for the bottom min(limit, 8) bits (our table only supports 8 bits). */
     277                 :      170265 :         limb_t m = (MAX_LIMB >> (LIMB_SIZE - limit)) & 255U;
     278                 :             :         /* Find what multiple of f must be added to g to cancel its bottom min(limit, 8) bits. */
     279                 :      170265 :         limb_t w = (g * NEGINV256[(f >> 1) & 127]) & m;
     280                 :             :         /* Do so. */
     281                 :      170265 :         g += f * w;
     282                 :      170265 :         q += u * w;
     283                 :      170265 :         r += v * w;
     284                 :      170265 :     }
     285                 :       27169 :     out.u = (signed_limb_t)u;
     286                 :       27169 :     out.v = (signed_limb_t)v;
     287                 :       27169 :     out.q = (signed_limb_t)q;
     288                 :       27169 :     out.r = (signed_limb_t)r;
     289                 :       27169 :     return eta;
     290                 :             : }
     291                 :             : 
     292                 :             : /** Apply matrix t/2^SIGNED_LIMB_SIZE to vector [d,e], modulo modulus.
     293                 :             :  *
     294                 :             :  * On input and output, d and e are in range 1-2*modulus..modulus-1.
     295                 :             :  */
     296                 :       27169 : inline void UpdateDE(Num3072Signed& d, Num3072Signed& e, const SignedMatrix& t)
     297                 :             : {
     298                 :       27169 :     const signed_limb_t u = t.u, v=t.v, q=t.q, r=t.r;
     299                 :             : 
     300                 :             :     /* [md,me] start as zero; plus [u,q] if d is negative; plus [v,r] if e is negative. */
     301                 :       27169 :     signed_limb_t sd = d.limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
     302                 :       27169 :     signed_limb_t se = e.limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
     303                 :       27169 :     signed_limb_t md = (u & sd) + (v & se);
     304                 :       27169 :     signed_limb_t me = (q & sd) + (r & se);
     305                 :             :     /* Begin computing t*[d,e]. */
     306                 :       27169 :     signed_limb_t di = d.limbs[0], ei = e.limbs[0];
     307                 :       27169 :     signed_double_limb_t cd = (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
     308                 :       27169 :     signed_double_limb_t ce = (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
     309                 :             :     /* Correct md,me so that t*[d,e]+modulus*[md,me] has SIGNED_LIMB_SIZE zero bottom bits. */
     310                 :       27169 :     md -= (MODULUS_INVERSE * limb_t(cd) + md) & MAX_SIGNED_LIMB;
     311                 :       27169 :     me -= (MODULUS_INVERSE * limb_t(ce) + me) & MAX_SIGNED_LIMB;
     312                 :             :     /* Update the beginning of computation for t*[d,e]+modulus*[md,me] now md,me are known. */
     313                 :       27169 :     cd -= (signed_double_limb_t)1103717 * md;
     314                 :       27169 :     ce -= (signed_double_limb_t)1103717 * me;
     315                 :             :     /* Verify that the low SIGNED_LIMB_SIZE bits of the computation are indeed zero, and then throw them away. */
     316                 :       27169 :     Assume((cd & MAX_SIGNED_LIMB) == 0);
     317                 :       27169 :     Assume((ce & MAX_SIGNED_LIMB) == 0);
     318                 :       27169 :     cd >>= SIGNED_LIMB_SIZE;
     319                 :       27169 :     ce >>= SIGNED_LIMB_SIZE;
     320                 :             :     /* Now iteratively compute limb i=1..SIGNED_LIMBS-2 of t*[d,e]+modulus*[md,me], and store them in output
     321                 :             :      * limb i-1 (shifting down by SIGNED_LIMB_SIZE bits). The corresponding limbs in modulus are all zero,
     322                 :             :      * so modulus/md/me are not actually involved here. */
     323         [ +  + ]:     1331281 :     for (int i = 1; i < SIGNED_LIMBS - 1; ++i) {
     324                 :     1304112 :         di = d.limbs[i];
     325                 :     1304112 :         ei = e.limbs[i];
     326                 :     1304112 :         cd += (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
     327                 :     1304112 :         ce += (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
     328                 :     1304112 :         d.limbs[i - 1] = (signed_limb_t)cd & MAX_SIGNED_LIMB; cd >>= SIGNED_LIMB_SIZE;
     329                 :     1304112 :         e.limbs[i - 1] = (signed_limb_t)ce & MAX_SIGNED_LIMB; ce >>= SIGNED_LIMB_SIZE;
     330                 :             :     }
     331                 :             :     /* Compute limb SIGNED_LIMBS-1 of t*[d,e]+modulus*[md,me], and store it in output limb SIGNED_LIMBS-2. */
     332                 :       27169 :     di = d.limbs[SIGNED_LIMBS - 1];
     333                 :       27169 :     ei = e.limbs[SIGNED_LIMBS - 1];
     334                 :       27169 :     cd += (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
     335                 :       27169 :     ce += (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
     336                 :       27169 :     cd += (signed_double_limb_t)md << FINAL_LIMB_MODULUS_BITS;
     337                 :       27169 :     ce += (signed_double_limb_t)me << FINAL_LIMB_MODULUS_BITS;
     338                 :       27169 :     d.limbs[SIGNED_LIMBS - 2] = (signed_limb_t)cd & MAX_SIGNED_LIMB; cd >>= SIGNED_LIMB_SIZE;
     339                 :       27169 :     e.limbs[SIGNED_LIMBS - 2] = (signed_limb_t)ce & MAX_SIGNED_LIMB; ce >>= SIGNED_LIMB_SIZE;
     340                 :             :     /* What remains goes into output limb SINGED_LIMBS-1 */
     341                 :       27169 :     d.limbs[SIGNED_LIMBS - 1] = (signed_limb_t)cd;
     342                 :       27169 :     e.limbs[SIGNED_LIMBS - 1] = (signed_limb_t)ce;
     343                 :       27169 : }
     344                 :             : 
     345                 :             : /** Apply matrix t/2^SIGNED_LIMB_SIZE to vector (f,g).
     346                 :             :  *
     347                 :             :  * The matrix t must be chosen such that t*(f,g) results in multiples of 2^SIGNED_LIMB_SIZE.
     348                 :             :  * This is the case for matrices computed by ComputeDivstepMatrix().
     349                 :             :  */
     350                 :       27169 : inline void UpdateFG(Num3072Signed& f, Num3072Signed& g, const SignedMatrix& t, int len)
     351                 :             : {
     352                 :       27169 :     const signed_limb_t u = t.u, v=t.v, q=t.q, r=t.r;
     353                 :             : 
     354                 :       27169 :     signed_limb_t fi, gi;
     355                 :       27169 :     signed_double_limb_t cf, cg;
     356                 :             :     /* Start computing t*[f,g]. */
     357                 :       27169 :     fi = f.limbs[0];
     358                 :       27169 :     gi = g.limbs[0];
     359                 :       27169 :     cf = (signed_double_limb_t)u * fi + (signed_double_limb_t)v * gi;
     360                 :       27169 :     cg = (signed_double_limb_t)q * fi + (signed_double_limb_t)r * gi;
     361                 :             :     /* Verify that the bottom SIGNED_LIMB_BITS bits of the result are zero, and then throw them away. */
     362                 :       27169 :     Assume((cf & MAX_SIGNED_LIMB) == 0);
     363                 :       27169 :     Assume((cg & MAX_SIGNED_LIMB) == 0);
     364                 :       27169 :     cf >>= SIGNED_LIMB_SIZE;
     365                 :       27169 :     cg >>= SIGNED_LIMB_SIZE;
     366                 :             :     /* Now iteratively compute limb i=1..SIGNED_LIMBS-1 of t*[f,g], and store them in output limb i-1 (shifting
     367                 :             :      * down by SIGNED_LIMB_BITS bits). */
     368         [ +  + ]:      953176 :     for (int i = 1; i < len; ++i) {
     369                 :      926007 :         fi = f.limbs[i];
     370                 :      926007 :         gi = g.limbs[i];
     371                 :      926007 :         cf += (signed_double_limb_t)u * fi + (signed_double_limb_t)v * gi;
     372                 :      926007 :         cg += (signed_double_limb_t)q * fi + (signed_double_limb_t)r * gi;
     373                 :      926007 :         f.limbs[i - 1] = (signed_limb_t)cf & MAX_SIGNED_LIMB; cf >>= SIGNED_LIMB_SIZE;
     374                 :      926007 :         g.limbs[i - 1] = (signed_limb_t)cg & MAX_SIGNED_LIMB; cg >>= SIGNED_LIMB_SIZE;
     375                 :             :     }
     376                 :             :     /* What remains is limb SIGNED_LIMBS of t*[f,g]; store it as output limb SIGNED_LIMBS-1. */
     377                 :       27169 :     f.limbs[len - 1] = (signed_limb_t)cf;
     378                 :       27169 :     g.limbs[len - 1] = (signed_limb_t)cg;
     379                 :             : 
     380                 :       27169 : }
     381                 :             : } // namespace
     382                 :             : 
     383                 :         270 : Num3072 Num3072::GetInverse() const
     384                 :             : {
     385                 :             :     // Compute a modular inverse based on a variant of the safegcd algorithm:
     386                 :             :     // - Paper: https://gcd.cr.yp.to/papers.html
     387                 :             :     // - Inspired by this code in libsecp256k1:
     388                 :             :     //   https://github.com/bitcoin-core/secp256k1/blob/master/src/modinv32_impl.h
     389                 :             :     // - Explanation of the algorithm:
     390                 :             :     //   https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md
     391                 :             : 
     392                 :             :     // Local variables d, e, f, g:
     393                 :             :     // - f and g are the variables whose gcd we compute (despite knowing the answer is 1):
     394                 :             :     //   - f is always odd, and initialized as modulus
     395                 :             :     //   - g is initialized as *this (called x in what follows)
     396                 :             :     // - d and e are the numbers for which at every step it is the case that:
     397                 :             :     //   - f = d * x mod modulus; d is initialized as 0
     398                 :             :     //   - g = e * x mod modulus; e is initialized as 1
     399                 :         270 :     Num3072Signed d, e, f, g;
     400                 :         270 :     e.limbs[0] = 1;
     401                 :             :     // F is initialized as modulus, which in signed limb representation can be expressed
     402                 :             :     // simply as 2^3072 + -MAX_PRIME_DIFF.
     403                 :         270 :     f.limbs[0] = -MAX_PRIME_DIFF;
     404                 :         270 :     f.limbs[FINAL_LIMB_POSITION] = ((limb_t)1) << FINAL_LIMB_MODULUS_BITS;
     405                 :         270 :     g.FromNum3072(*this);
     406                 :         270 :     int len = SIGNED_LIMBS; //!< The number of significant limbs in f and g
     407                 :         270 :     signed_limb_t eta = -1; //!< State to track knowledge about ratio of f and g
     408                 :             :     // Perform divsteps on [f,g] until g=0 is reached, keeping (d,e) synchronized with them.
     409                 :       27169 :     while (true) {
     410                 :             :         // Compute transformation matrix t that represents the next SIGNED_LIMB_SIZE divsteps
     411                 :             :         // to apply. This can be computed from just the bottom limb of f and g, and eta.
     412                 :       27169 :         SignedMatrix t;
     413                 :       27169 :         eta = ComputeDivstepMatrix(eta, f.limbs[0], g.limbs[0], t);
     414                 :             :         // Apply that transformation matrix to the full [f,g] vector.
     415                 :       27169 :         UpdateFG(f, g, t, len);
     416                 :             :         // Apply that transformation matrix to the full [d,e] vector (mod modulus).
     417                 :       27169 :         UpdateDE(d, e, t);
     418                 :             : 
     419                 :             :         // Check if g is zero.
     420         [ +  + ]:       27169 :         if (g.limbs[0] == 0) {
     421                 :             :             signed_limb_t cond = 0;
     422         [ +  + ]:      533070 :             for (int j = 1; j < len; ++j) {
     423                 :      522144 :                 cond |= g.limbs[j];
     424                 :             :             }
     425                 :             :             // If so, we're done.
     426         [ +  + ]:       10926 :             if (cond == 0) break;
     427                 :             :         }
     428                 :             : 
     429                 :             :         // Check if the top limbs of both f and g are both 0 or -1.
     430                 :       26899 :         signed_limb_t fn = f.limbs[len - 1], gn = g.limbs[len - 1];
     431                 :       26899 :         signed_limb_t cond = ((signed_limb_t)len - 2) >> (LIMB_SIZE - 1);
     432                 :       26899 :         cond |= fn ^ (fn >> (LIMB_SIZE - 1));
     433                 :       26899 :         cond |= gn ^ (gn >> (LIMB_SIZE - 1));
     434         [ +  + ]:       26899 :         if (cond == 0) {
     435                 :             :             // If so, drop the top limb, shrinking the size of f and g, by
     436                 :             :             // propagating the sign to the previous limb.
     437                 :       13230 :             f.limbs[len - 2] |= (limb_t)f.limbs[len - 1] << SIGNED_LIMB_SIZE;
     438                 :       13230 :             g.limbs[len - 2] |= (limb_t)g.limbs[len - 1] << SIGNED_LIMB_SIZE;
     439                 :       13230 :             --len;
     440                 :             :         }
     441                 :             :     }
     442                 :             :     // At some point, [f,g] will have been rewritten into [f',0], such that gcd(f,g) = gcd(f',0).
     443                 :             :     // This is proven in the paper. As f started out being modulus, a prime number, we know that
     444                 :             :     // gcd is 1, and thus f' is 1 or -1.
     445                 :         270 :     Assume((f.limbs[0] & MAX_SIGNED_LIMB) == 1 || (f.limbs[0] & MAX_SIGNED_LIMB) == MAX_SIGNED_LIMB);
     446                 :             :     // As we've maintained the invariant that f = d * x mod modulus, we get d/f mod modulus is the
     447                 :             :     // modular inverse of x we're looking for. As f is 1 or -1, it is also true that d/f = d*f.
     448                 :             :     // Normalize d to prepare it for output, while negating it if f is negative.
     449                 :         270 :     d.Normalize(f.limbs[len - 1] >> (LIMB_SIZE  - 1));
     450                 :         270 :     Num3072 ret;
     451                 :         270 :     d.ToNum3072(ret);
     452                 :         270 :     return ret;
     453                 :             : }
     454                 :             : 
     455                 :         880 : void Num3072::Multiply(const Num3072& a)
     456                 :             : {
     457                 :         880 :     limb_t c0 = 0, c1 = 0, c2 = 0;
     458                 :         880 :     Num3072 tmp;
     459                 :             : 
     460                 :             :     /* Compute limbs 0..N-2 of this*a into tmp, including one reduction. */
     461         [ +  + ]:       42240 :     for (int j = 0; j < LIMBS - 1; ++j) {
     462                 :       41360 :         limb_t d0 = 0, d1 = 0, d2 = 0;
     463                 :       41360 :         mul(d0, d1, this->limbs[1 + j], a.limbs[LIMBS + j - (1 + j)]);
     464         [ +  + ]:      992640 :         for (int i = 2 + j; i < LIMBS; ++i) muladd3(d0, d1, d2, this->limbs[i], a.limbs[LIMBS + j - i]);
     465                 :       41360 :         mulnadd3(c0, c1, c2, d0, d1, d2, MAX_PRIME_DIFF);
     466         [ +  + ]:     1034000 :         for (int i = 0; i < j + 1; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[j - i]);
     467                 :       41360 :         extract3(c0, c1, c2, tmp.limbs[j]);
     468                 :             :     }
     469                 :             : 
     470                 :             :     /* Compute limb N-1 of a*b into tmp. */
     471         [ +  - ]:         880 :     assert(c2 == 0);
     472         [ +  + ]:       43120 :     for (int i = 0; i < LIMBS; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[LIMBS - 1 - i]);
     473                 :         880 :     extract3(c0, c1, c2, tmp.limbs[LIMBS - 1]);
     474                 :             : 
     475                 :             :     /* Perform a second reduction. */
     476                 :         880 :     muln2(c0, c1, MAX_PRIME_DIFF);
     477         [ +  + ]:       43120 :     for (int j = 0; j < LIMBS; ++j) {
     478                 :       42240 :         addnextract2(c0, c1, tmp.limbs[j], this->limbs[j]);
     479                 :             :     }
     480                 :             : 
     481         [ -  + ]:         880 :     assert(c1 == 0);
     482         [ -  + ]:         880 :     assert(c0 == 0 || c0 == 1);
     483                 :             : 
     484                 :             :     /* Perform up to two more reductions if the internal state has already
     485                 :             :      * overflown the MAX of Num3072 or if it is larger than the modulus or
     486                 :             :      * if both are the case.
     487                 :             :      * */
     488         [ +  + ]:         880 :     if (this->IsOverflow()) this->FullReduce();
     489         [ -  + ]:         880 :     if (c0) this->FullReduce();
     490                 :         880 : }
     491                 :             : 
     492                 :        2192 : void Num3072::SetToOne()
     493                 :             : {
     494                 :        2192 :     this->limbs[0] = 1;
     495         [ +  + ]:      105216 :     for (int i = 1; i < LIMBS; ++i) this->limbs[i] = 0;
     496                 :        2192 : }
     497                 :             : 
     498                 :         270 : void Num3072::Divide(const Num3072& a)
     499                 :             : {
     500         [ +  + ]:         270 :     if (this->IsOverflow()) this->FullReduce();
     501                 :             : 
     502                 :         270 :     Num3072 inv{};
     503         [ -  + ]:         270 :     if (a.IsOverflow()) {
     504                 :           0 :         Num3072 b = a;
     505                 :           0 :         b.FullReduce();
     506                 :           0 :         inv = b.GetInverse();
     507                 :             :     } else {
     508                 :         270 :         inv = a.GetInverse();
     509                 :             :     }
     510                 :             : 
     511                 :         270 :     this->Multiply(inv);
     512         [ -  + ]:         270 :     if (this->IsOverflow()) this->FullReduce();
     513                 :         270 : }
     514                 :             : 
     515                 :         390 : Num3072::Num3072(const unsigned char (&data)[BYTE_SIZE]) {
     516         [ +  + ]:       19110 :     for (int i = 0; i < LIMBS; ++i) {
     517                 :       18720 :         if (sizeof(limb_t) == 4) {
     518                 :             :             this->limbs[i] = ReadLE32(data + 4 * i);
     519                 :       18720 :         } else if (sizeof(limb_t) == 8) {
     520                 :       18720 :             this->limbs[i] = ReadLE64(data + 8 * i);
     521                 :             :         }
     522                 :             :     }
     523                 :         390 : }
     524                 :             : 
     525                 :         270 : void Num3072::ToBytes(unsigned char (&out)[BYTE_SIZE]) {
     526         [ +  + ]:       13230 :     for (int i = 0; i < LIMBS; ++i) {
     527                 :       12960 :         if (sizeof(limb_t) == 4) {
     528                 :             :             WriteLE32(out + i * 4, this->limbs[i]);
     529                 :       12960 :         } else if (sizeof(limb_t) == 8) {
     530                 :       12960 :             WriteLE64(out + i * 8, this->limbs[i]);
     531                 :             :         }
     532                 :             :     }
     533                 :         270 : }
     534                 :             : 
     535                 :         390 : Num3072 MuHash3072::ToNum3072(Span<const unsigned char> in) {
     536                 :         390 :     unsigned char tmp[Num3072::BYTE_SIZE];
     537                 :             : 
     538                 :         390 :     uint256 hashed_in{(HashWriter{} << in).GetSHA256()};
     539                 :         390 :     static_assert(sizeof(tmp) % ChaCha20Aligned::BLOCKLEN == 0);
     540                 :         390 :     ChaCha20Aligned{MakeByteSpan(hashed_in)}.Keystream(MakeWritableByteSpan(tmp));
     541                 :         390 :     Num3072 out{tmp};
     542                 :             : 
     543                 :         390 :     return out;
     544                 :             : }
     545                 :             : 
     546                 :         186 : MuHash3072::MuHash3072(Span<const unsigned char> in) noexcept
     547                 :             : {
     548                 :         186 :     m_numerator = ToNum3072(in);
     549                 :         186 : }
     550                 :             : 
     551                 :         270 : void MuHash3072::Finalize(uint256& out) noexcept
     552                 :             : {
     553                 :         270 :     m_numerator.Divide(m_denominator);
     554                 :         270 :     m_denominator.SetToOne();  // Needed to keep the MuHash object valid
     555                 :             : 
     556                 :         270 :     unsigned char data[Num3072::BYTE_SIZE];
     557                 :         270 :     m_numerator.ToBytes(data);
     558                 :             : 
     559                 :         270 :     out = (HashWriter{} << data).GetSHA256();
     560                 :         270 : }
     561                 :             : 
     562                 :         112 : MuHash3072& MuHash3072::operator*=(const MuHash3072& mul) noexcept
     563                 :             : {
     564                 :         112 :     m_numerator.Multiply(mul.m_numerator);
     565                 :         112 :     m_denominator.Multiply(mul.m_denominator);
     566                 :         112 :     return *this;
     567                 :             : }
     568                 :             : 
     569                 :          91 : MuHash3072& MuHash3072::operator/=(const MuHash3072& div) noexcept
     570                 :             : {
     571                 :          91 :     m_numerator.Multiply(div.m_denominator);
     572                 :          91 :     m_denominator.Multiply(div.m_numerator);
     573                 :          91 :     return *this;
     574                 :             : }
     575                 :             : 
     576                 :         203 : MuHash3072& MuHash3072::Insert(Span<const unsigned char> in) noexcept {
     577                 :         203 :     m_numerator.Multiply(ToNum3072(in));
     578                 :         203 :     return *this;
     579                 :             : }
     580                 :             : 
     581                 :           1 : MuHash3072& MuHash3072::Remove(Span<const unsigned char> in) noexcept {
     582                 :           1 :     m_denominator.Multiply(ToNum3072(in));
     583                 :           1 :     return *this;
     584                 :             : }
        

Generated by: LCOV version 2.0-1