Optimized modular multiplicative inverse for Bitcoin (secp256k1)

Question

I have written an application to brute-force attack Bitcoin addresses for OpenCL. It implements a simple exhaustive search starting from public key G (private key 1) and point increment (addition of G) to get next key.

It works, but its performance is one order of magnitude slower than its CPU implementation using GMP (same point increment algorithm but with GMP routines instead of custom ones). On GMP implementation, hashing is the bottleneck. On OpenCL implementation, point increment is the bottleneck.

My point increment implementation uses modular multiplicative inverse via binary extended GCD (algorithm 14.61 of Handbook of Applied Cryptography). It is one order of magnitude slower than GMP implementation using mpz_invert().

What I need is an optimized point increment (point addition with G) optimized for curve secp256k1. This is an educational project, so my knowledge about the subject is limited. I don't know if there is a way to accomplish it without using modular multiplicative inverse or if there is a way to speed up this calculation for my particular application.

I saw and implementation (by Segher Boessenkool) using modular exponentiation with precomputed values, but there was no information on how to get those values for other curves. I failed to replicate GMP implementation. Its source code is huge and a mess of mixed optimized C and assembly code impossible to understand for me.

Requirements:

• A fast point addition by G implementation.
• Safety is irrelevant due to its use.
• If possible, optimized for curve secp256k1.
• If possible, precompute as much as possible.

Question:

How can I improve speed of elliptic curve point addition in my specific case?

---

Notes on comments:

• The search is sequential inside each OpenCL kernel, but parallel considering the whole program. My worksize is optimized experimentally, getting close to theoretical GPU FLOPS performing SHA-256 only, so the bottleneck is not in OpenCL usage, but algorithmic implementation of elliptic curve cryptography.
• The performance one order of magnitude smaller I am talking about is measured on CPU, not GPU. I get about 30 kKeys/s with my own modular multiplicative inverse and 200 kKeys/s (a laptop) with GMP. On GPU I get about 2 MKeys/s. Computing unoptimized SHA-256+RIPEMD-160 hash only, without point addition, I get about 1 GHashes/s (AMD Radeon RX 580).
• Main program computes an array of public keys using scalar multiplication. It can compute them pseudo-randomly or equally spaced. This array is feed to the OpenCL kernel, and each instance computes RIPEMD160(SHA256(P)) == TARGET, RIPEMD160(SHA256(P+G)) == TARGET, RIPEMD160(SHA256(P+2G)) == TARGET, ... For now, there is only a single target address.
• OpenCL is not the issue here. Profiling on CPU reveals that 90% of the time is expend computing modular multiplicative inverse. Same algorithm with GMP mpz_invert() uses only 25% of the time to do same step, using the remaining time on hash computation.

Answer 1

Half-extended binary GCD

The extended binary GCD of the HAC's algorithm 14.61 mentioned in the question, when given positive integers $x$ and $y$ , computes integers $a$ , $b$ , $v$ with $ax+by=v=\gcd(x,y)$ . When $v=1$, it follows that $a$ is the multiplicative inverse of $x$ .

In our context, $0<x<y$, $y$ odd, $\gcd(x,y)=1$ (the last two conditions are met since $y$ is $p=2^{256}-2^{32}-977$ of secp256k1$v=1$ , which is prime); thus $a$ is always the inverse. Importantly, we do not need $b$ , which allows simplifications. This leads to the half-extended binary GCD, which I hope can help to some degree:

• If $x$ is odd then $u\gets x$ ; else $u\gets x+y$ ;
• $v\gets y$ ; $a\gets0$ ; $d\gets 1$ ;   [or $d\gets y-1$ in unsigned variant, see note]
• While $v\ne1$
  [invariant here and at start of the next while loop: $u$ and $v$ are odd and distinct]
  • While $v<u$
    • $u\gets u-v$ ; $d\gets d-a$ ;   [or $d\gets d+a$ in unsigned variant, see note]
    • While $u$ is even (that's at least once)
      • If $d$ is odd then $d\gets d+y$ ;
      • $u\gets u/2$ ; $d\gets d/2$ ;
  • $v\gets v-u$ ; $a\gets a-d$ ;   [or $a\gets a+d$ in unsigned variant, see note]
  • While $v$ is even (that's at least once)
    • If $a$ is odd then $a\gets a+y$ ;
    • $v\gets v/2$ ; $a\gets a/2$ ;
• $a\gets a\bmod y$ ; that's the desired inverse.

Note: As shown, $a$ and $d$ are signed. In an unsigned variant we can keep them non-negative, by making the three changes noted.

As shown, $a$ and $d$ can grow to several times $y$ (in absolute value for the signed variant). We can devote some extra bits to this. I conjecture without proof that $\max(|a|,|d|)<4y\log_2(y)$ (or $\max(a,d)<8y\log_2(y)$ for the unsigned variant). Alternatively, we can change $a\gets a+y$ into $a\gets a-y$ when $a\ge0$ for the signed variant (for the unsigned one, that can be when $a\ge y'$ for some convenient $y'\ge y$, e.g. $y=2^{256}$ ); and same for $d\gets d-y$ . This is enough to ensure everything stays below $4y'$, or $2^{258}$ in the context.

Possible optimization: the number of times $i$ that the "While $u$ is even" loop is run is predictable from the low-order bits of $u$, and we can group the overall effect on $u$ and $d$ into a single loop. What needs to be added to $d$ depends on $i$, the $i$ low-order bits of $d$, and $y$ . It can be precomputed since $y$ is fixed. Even if we do this for $i$ up to 2 or 3 only (then iterate if necessary), that saves significantly (on a classical CPU at least). Same for $v$ and $a$.

The (half)-extend Euclidean algorithm (see this answer) also computes the multiplicative inverse, and has a simple variant avoiding multiplication (the quotient is made a power of two). It can be competitive. The two can be blended.

---

GPUs are best at SIMD workloads

GPUs shine at workloads where the same instruction flow applies to different independent jobs launched in parallel. Taking different code paths in simultaneously launched jobs at best dramatically reduces performance. Thus GPU-friendly algorithms are those with long sequences without tests, and tests which predictably mostly take the same path in different jobs. Algorithms often need to be tweaked towards that.

[hopefully more to come there]

---

It's a brute force attack

The attack require about $2^{b-1}$ point addition to find a private key with $b$ unknown bits, or about $256-b$ known bits (e.g. at zero). I see no reason why there should be such keys, other than deliberately. Thus I'm doubtful on the overall usefulness of the thing. We can still think about it for educational purposes.

---

Starting from addresses forces using Cartesian coordinates

Because a version of the attack is said to be hash-bound, we know that the attack tries to find a private key from the hash (known as address) of the corresponding public key, rather than from the public key. This has been confirmed in comments of that other answer, which proposes a mathematical technique assuming the public key is known. I'm not exploring that avenue.

Thus, the bulk of the sequential version of the algorithm is to compute public keys $kG$ for incremental private key $k$ using point addition of $G$ to move from one public key to the other; hash the resulting public key (essentially) in its Cartesian uncompressed $X\|Y$ form; and check the hash.

---

The attack should be multi-target

There is a tremendous potential benefit in targeting multiple addresses simultaneously; only the comparison of the final hash is made slightly more compute intensive (it becomes a search among the addresses targeted), and all the rest is shared. That's too good to ignore.

---

Outer loop for GPU (it's fine)

When optimizing an algorithm for GPU or other highly parallel architecture, it is essential to split the task into independent jobs. That part of the question made me think it was not done:

search starting from public key G (private key 1) and point increment (addition of G) to get next key

With this algorithm, the overall search speed would not scale up with more GPU power: point increment is going to be a bottleneck no matter how optimized it is.

But the question now also has:

Main program computes an array of public keys using scalar multiplication.

and that suggests "starting from public key G (private key 1)" does not apply to the OpenCL jobs/kernels. Instead, the private key search interval is subdivided, the corresponding public key computed, and the search performed; that's fine.

Answer 2

Actually, if the reason you are computing elliptic curve points is to compare them to a target one, you don't need to compute a multiplicative inverse. Instead, if you keep the point your iterating in projective coordinates $H = (H_x, H_y, H_z$) and the target value in standard coordinates $T = (T_x, T_y)$, you can determine whether $H = \pm T$ simply by checking whether $H_x = T_x H_z$ (actually, if you're keeping $H$ in Jacobean coordinates, this is $H_x = T_x H_z^2$ instead). As the additional overhead from doing computations with projective coordinates is less than the overhead from doing a modular inverse, this is a win.

Answer 3

GMP also uses extended GCD. It is fast precisely because it uses a "a mess of mixed optimized C and assembly code". You can't expect to achieve similar performance without the same optimizations.

Alternatively, your implementation could be particularly slow, it's impossible to tell without access to it.

Answer 4

An answer is a combination of answers 1 and 2. Compute an array of consecutive points in projective coordinates (only needs addition, multiplication, and squaring), multiply them all together, perform a single modular multiplicative inverse of the product, and manipulate the result (see references) so you can transform points back to affine coordinates to hash them.

Overall, the cost of this algorithm is 13 multiplications, 2 squares, and 6 additions per point plus a single modular multiplicative inverse per batch. My actual implementation runs 10 times faster on CPU. There are more room for improvements, but this is not related to cryptography anymore.

Answer 5

I have implemented inverse by modulo algorithm from this paper and modified it to work with overflowing and underflowing uint256 type (on the first iteration both numbers are considered positive, later top bit is use to compare their signs): https://www.researchgate.net/publication/304417579_Modular_Inverse_Algorithms_Without_Multiplications_for_Cryptographic_Applications

function inverseShift(uint256 a, uint256 m) public pure returns(uint256) {
    unchecked {
        uint256 u;
        uint256 v;
        uint256 r;
        uint256 s;
        if (a < m) { u = m; v = a; r = 0; s = 1; }
        else       { v = m; u = a; s = 0; r = 1; }

        uint256 llu = ll(u, true);
        uint256 llv = ll(v, true);
        for (uint i = 0; v + 1 > 2; i++) { // not 0 +/-1
            uint256 f = llu - llv;
            if (i == 0 || (u >> 255) == (v >> 255)) {
                u = u - (v << f);
                r = r - (s << f);
            }
            else {
                u = u + (v << f);
                r = r + (s << f);
            }

            llu = ll(u);
            if (llu < llv) {
                (u,v,r,s,llu,llv) = (v,u,s,r,llv,llu);
            }
        }
        if (v == 0) { return 0; }
        if (v >> 255 == 1) { s = 0-s; }
        if (s > m) { return s - m; }
        if (s >> 255 == 1) { return s + m; }
        return s;
    }
}

Moreover pretty efficient but not yet tested uint256 recursive implementation (without divisions) can be found here: https://gist.github.com/k06a/fd96f05d38bd25b232d13a027822de0c