How to calculate a mod-inverse: Two different ways

[NotATether]

I found myself having to write cryptography code on C++ and I needed to invert a key. Well, libsecp256k1 does not support that natively, so I was forced to find a different way. In the process, I have tried and successfully implemented two different methods you can use to invert your private keys and perform division on them if necessary.

Definition of mod-inverse

Modular Multiplicative Inverse in cryptography is when you find a number inside a group of numbers, that multiplies with another number to give 1. Then those two numbers are modular inverses of each other.

I will give you an example with small numbers. Lets suppose that our number line only supports:

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10

(edit: The first number was supposed to be zero not nine).

And then we "bend" the line into a circle, that is, if you add 10+1 you you would go back to 0, and vice versa for 0 - 1.

Now suppose we want to multiply 5*9. What would that give us in this group? Well normally, we would get 45, but that number is not defined here. So what can we do?

Well, let's see what happens if we subtract "1" to it. Then we will have 44, which happens to be evenly divisible by 11 (representing the total number of elements in the group), and 44 mod 11 is zero - you can open any calculator to verify this.

So we have just proven that 5 and 9 are multiplicative inverses (modulo 11 elements). We have found a very powerful way of finding mod-inverses by expressing them as:

a*x + m*y = gcd(a, m) = 1

Where we are trying to find the inverse of a mod m, and x and y are unknown - in an elliptic curve, A and M are not supposed to have any common factors, hence why it's 1.

Y is not very important here - it just tells us how many times you need to add or subtract m from a*x to get 1. It is X we are really interested in, for this is the modular inverse of A.

Hence it gives rise to the first solution to find modular inverse:

1. Using EGCD

EGCD stands for "extended greatest common divisor" but most people just call it the "extended Euclidean algorithm" because it is just an extension of Euler's ancient "sieves" algorithm for finding greatest common factors of two numbers. As I have just demonstrated above, EGCD also finds the multiplicative inverse of the number and also the "shift" that needs to be applied to the group size number to get the GCD, which, in crypto's case, is just 1.

Below I provide an EGCD algorithm you can use with any integral type with high enough resolution (here I use my own fixed-point library available on github at: https://github.com/ZenulAbidin/xFD). This code is BSD 3-clause licensed and is easily modified to support generalized EGCD by returning all the variables in an array or as reference parameters.

Code:

// Returns x such that xa + by = gcd(a, b) (= 1 for prime b)
Decimal MMI(const Decimal& a, const Decimal& b) {
    auto old_r = a, r = b;
    auto old_s = 1_D, s = 0_D;
    auto old_t = 0_D, t = 1_D;
    Decimal tmp;

    while (r != 0_D) {
        auto quotient = xFD::Floor(old_r / r);
        tmp = r;
        r = old_r - quotient * r;
        old_r = tmp;
        tmp = s;
        s = old_s - quotient * s;
        old_s = tmp;
        tmp = t;
        t = old_t - quotient * t;
        old_t = tmp;
    }

    // s, t = a/gcd(a,b) and b/gcd(a,b) respectively
    // old_r is gcd
    // r is zero obviously
    // old_s, old_t are the Bézout coefficients x and y for a and b respectively.

    // The multiplicative inverse should not be negative.
    return (old_s < 0_D) ? old_s + b : old_s; 
}

(here, m is just swapped for a different name "b").

2. Using Fermat's little theorem

This was described to me by a fellow Bitcoiner on this forum, and arises from the fact that at the curve order P:

a^P === a(^1)

This identity does not look very interesting by itself, until you see what happens when we take some arithmetic on both sides:

a^(P-1) ===a^(1-1) === a^0 === 1

This is called Euler's theorem (and it only works when a and m have a GCD of 1). It's not very useful for us either, but if we apply this again:

a^(P-2) ===a^(1-2) === a^-1

Then we get the multiplicative inverse of a without EGCD.

Now naturally, this only works if you already have an elliptic curve library with support for exponentiation, which is really just repeated self-multiplication. It is also slower than the first method, but one big advantage of this algo is that there are no branches with different timings, making this algo very secure against side-channel attacks.


So pick your poison. Writing a wallet? Use exponentiation. EGCD is much faster, but is unsuitable for use in wallets.

Additional reading:

https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

[1714744831]

1714744831

[1714744831]

1714744831

[ymgve2]

With Fermat's little theorem, you have no guarantees about the immunity to side channel attacks, because it all relies on the underlying library's exponentiation function being safe against side channel attacks. Most math libraries are not safe against this.

[COBRAS]

I found myself having to write cryptography code on C++ and I needed to invert a key. Well, libsecp256k1 does not support that natively, so I was forced to find a different way. In the process, I have tried and successfully implemented two different methods you can use to invert your private keys and perform division on them if necessary.

Definition of mod-inverse

Modular Multiplicative Inverse in cryptography is when you find a number inside a group of numbers, that multiplies with another number to give 1. Then those two numbers are modular inverses of each other.

I will give you an example with small numbers. Lets suppose that our number line only supports:

9, 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10

And then we "bend" the line into a circle, that is, if you add 10+1 you you would go back to 0, and vice versa for 0 - 1.

Now suppose we want to multiply 5*9. What would that give us in this group? Well normally, we would get 45, but that number is not defined here. So what can we do?

Well, let's see what happens if we subtract "1" to it. Then we will have 44, which happens to be evenly divisible by 11 (representing the total number of elements in the group), and 44 mod 11 is zero - you can open any calculator to verify this.

So we have just proven that 5 and 9 are multiplicative inverses (modulo 11 elements). We have found a very powerful way of finding mod-inverses by expressing them as:

a*x + m*y = gcd(a, m) = 1

Where we are trying to find the inverse of a mod m, and x and y are unknown - in an elliptic curve, A and M are not supposed to have any common factors, hence why it's 1.

Y is not very important here - it just tells us how many times you need to add or subtract m from a*x to get 1. It is X we are really interested in, for this is the modular inverse of A.

Hence it gives rise to the first solution to find modular inverse:

1. Using EGCD

EGCD stands for "extended greatest common divisor" but most people just call it the "extended Euclidean algorithm" because it is just an extension of Euler's ancient "sieves" algorithm for finding greatest common factors of two numbers. As I have just demonstrated above, EGCD also finds the multiplicative inverse of the number and also the "shift" that needs to be applied to the group size number to get the GCD, which, in crypto's case, is just 1.

Below I provide an EGCD algorithm you can use with any integral type with high enough resolution (here I use my own fixed-point library available on github at: https://github.com/ZenulAbidin/xFD). This code is BSD 3-clause licensed and is easily modified to support generalized EGCD by returning all the variables in an array or as reference parameters.

Code:

// Returns x such that xa + by = gcd(a, b) (= 1 for prime b)
Decimal MMI(const Decimal& a, const Decimal& b) {
    auto old_r = a, r = b;
    auto old_s = 1_D, s = 0_D;
    auto old_t = 0_D, t = 1_D;
    Decimal tmp;

    while (r != 0_D) {
        auto quotient = xFD::Floor(old_r / r);
        tmp = r;
        r = old_r - quotient * r;
        old_r = tmp;
        tmp = s;
        s = old_s - quotient * s;
        old_s = tmp;
        tmp = t;
        t = old_t - quotient * t;
        old_t = tmp;
    }

    // s, t = a/gcd(a,b) and b/gcd(a,b) respectively
    // old_r is gcd
    // r is zero obviously
    // old_s, old_t are the Bézout coefficients x and y for a and b respectively.

    // The multiplicative inverse should not be negative.
    return (old_s < 0_D) ? old_s + b : old_s; 
}

(here, m is just swapped for a different name "b").

2. Using Fermat's little theorem

This was described to me by a fellow Bitcoiner on this forum, and arises from the fact that at the curve order P:

a^P === a(^1)

This identity does not look very interesting by itself, until you see what happens when we take some arithmetic on both sides:

a^(P-1) ===a^(1-1) === a^0 === 1

This is called Euler's theorem (and it only works when a and m have a GCD of 1). It's not very useful for us either, but if we apply this again:

a^(P-2) ===a^(1-2) === a^-1

Then we get the multiplicative inverse of a without EGCD.

Now naturally, this only works if you already have an elliptic curve library with support for exponentiation, which is really just repeated self-multiplication. It is also slower than the first method, but one big advantage of this algo is that there are no branches with different timings, making this algo very secure against side-channel attacks.


So pick your poison. Writing a wallet? Use exponentiation. EGCD is much faster, but is unsuitable for use in wallets.

Additional reading:

https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

this work for publick keys ?

I was try sub and div. Sub is more powerfull for daungeade range.


some examples of my codes, if you need pubs replace priv to pubs and  get ranges smaler then original pubkey.

Code:

N =    115792089237316195423570985008687907852837564279074904382605163141518161494337

def inv(v): return pow(v, N-2, N)
def divnum(a, b): return ( (a * inv(b) ) % N )
i= 0#2**10#2**10#*4275#219788
priv=0xa333d84649e1afa9c17a74c04f6c32  #- (0x634adc*2**56)#2***70
#-5-0x7962e67f4ee462dbb2d11009e66
#-7 x 0xe5a300a670ed5ab8dce16c1530e
#0x1ceba4b7b228839f71ed9af 2^00
#######0211220xa333d84649e1afa9c17a74c04f6c32=0x2fe0a51aefdd932cbce3d - 2^94
0x378fc97ca23fd4356c48969
0x246aab53748a772e927d3023efbb#Pr 0xa333d84649e1afa9c17a74c04f6c32
x =0

shitnim = 33554432
z=1
while i <  shitnim:
   
   s =(divnum(priv,2**94)-(shitnim)) %N
   while z < 1000:
        d= (divnum(z,2**94)) %N
        x = divnum(s,d)
        z = z +1
        
        if x <2**256and x >= 0:
            print("New Priv x",hex(x),z,i)
            
            
   i = i +1

Code:

from random import randint

N =    115792089237316195423570985008687907852837564279074904382605163141518161494337

def inv(v): return pow(v, N-2, N)
def divnum(a, b): return ( (a * inv(b) ) % N )

i=0
input = 0x9fd24b3abe244d6c443df56fa494dc

#input =  0x5f87127# +1

delta = 6#2+2#+2+2+2+2+2+2+2+2+2+2+2+2
print(delta)

gamma = 2
z =0
d1= 80

while i <= input:
   z=1
   i = i +1
   while z <= 1000:
      
    
     d= (divnum(input-z,delta)) 
     s = divnum(i,gamma) %N
     result = divnum(d,s)
     z = z + 1
     #print(i,z)  
     
     
     
    
    
     if result <=input and result >=0:
        print("result",hex(result),"i",hex(i),"input",hex(input),z)

Enother interesting , imho, is using property then div a / x = 1 only if a=x, in this way posible found a/x = in less number of operation then a.

ex:

Code:

from random import randint

N =    115792089237316195423570985008687907852837564279074904382605163141518161494337

def inv(v): return pow(v, N-2, N)
def divnum(a, b): return ( (a * inv(b) ) % N )

i=0
#input2^^120 = 0x9fd24b3abe244d6c443df56fa494dc

input = 0x5f87 +1

delta = 12

gamma = 2

d1= 80

while i < 2**61:
    d= (divnum(input,delta)) 
    s = divnum(i,gamma) %N
    result = divnum(d,s)
    
    if result <input and result >=0:
        print("result",hex(result),"i",hex(i),"input",hex(input))
        
    i = i +1

result:

Code:

result 0x0 i 0x0 input 0x5f88
result 0xfec i 0x1 input 0x5f88
result 0x7f6 i 0x2 input 0x5f88
result 0x3fb i 0x4 input 0x5f88
result 0x4 i 0x3fb input 0x5f88
result 0x2 i 0x7f6 input 0x5f88
result 0x1 i 0xfec input 0x5f88

as you see result = 1 geted in 0xfec operation instaed input 0x5f88 operation.

basepoint for bsgs or kangaroo vill = inversion(gamma)




example with modified base point:

Code:

priv 0x953f49ae70ed676   2**60

           953ec9ac70ed676
           753f49ae70ed676

/ 0x20000000 = dcf2e84ffffffffffffffffffffffffee73997ebdf8a334b693c56a1461840e7


- 1 - 0x4000 =
dcf2e84ffffffffffffffffffffffffee73997ebdf8a334b693c56a1461800e6

 pub =  0217e6a9925b659b97394039c7f355150a7be9f3adcf797bc5681f8584af0a69a6







- -_-------_---
2**29:
    
0x20000000

inv (2**29)
=

ac4589f7ffffffffffffffffffffffff25151e38f58f7d8301173ec3803d864e
    
    

G:
    
 02:   
    590670df6004cc1f1d1dda2f08a1045b77349c40409a63187336705ec5f292e4
    
04:
    
590670df6004cc1f1d1dda2f08a1045b77349c40409a63187336705ec5f292e4

b3f7df737c4ee700fd1ab386939343d70227b48df0e7ba438549695e30ea7578

a*x + m*y = 1

how to use this with pubkeys ? 1 is pubkey from priv 1, this is no problem .
a and m is pubkeys too(a is unknown for ex and m i known. x and y are need to find, but they are not unknown....

. what about way to find

?

[moderator's note: consecutive posts merged]

[NotATether]

I will give you an example with small numbers. Lets suppose that our number line only supports:

9, 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10

And then we "bend" the line into a circle, that is, if you add 10+1 you you would go back to 0, and vice versa for 0 - 1.

I'm a bit confused. Did you mean to write "0, 1, ..., 9 and 10"?

Yes, that was a typo. "0" and "9" are right next to each other on the keyboard. Thanks for pointing that out.

With Fermat's little theorem, you have no guarantees about the immunity to side channel attacks, because it all relies on the underlying library's exponentiation function being safe against side channel attacks. Most math libraries are not safe against this.

That is correct.

I am now wondering whether it is possible to eliminate side-channels in EGCD at all, considering that the loop uses expensive division and multiplication ops followed by cheap subtraction ops that could leak some info about keys via power frequency.

[COBRAS]

then divide

Code:

./md 11 / 4
Result: 3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a3340d9053

xman@localhost:~/ecctools$ ./md 1 / 4
Result: bfffffffffffffffffffffffffffffff0c0325ad0376782ccfddc6e99c28b0f1

xman@localhost:~/ecctools$ ./md 0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a3340d9053 + 0xbfffffffffffffffffffffffffffffff0c0325ad0376782ccfddc6e99c28b0f1
Result: 3

xman@localhost:~/ecctools$ ./md 99 / 50
Result: 75c28f5c28f5c28f5c28f5c28f5c28f52cea09745aded9113ea35473f960a325

xman@localhost:~/ecctools$ ./md 1 / 50
Result: 8a3d70a3d70a3d70a3d70a3d70a3d7098dc4d3725469c72a812f0a18d6d59e1e

xman@localhost:~/ecctools$ ./md 0x75c28f5c28f5c28f5c28f5c28f5c28f52cea09745aded9113ea35473f960a325 + 0x8a3d70a3d70a3d70a3d70a3d70a3d7098dc4d3725469c72a812f0a18d6d59e1e
Result: 2

then div 99 / 50 result = 50 + 49/50.so fo find 99 need, if dont known divisision (49), so calc from (1/50 to 50/50 ),+ 50 and privkey 99 will bi finded in 49 operations... mach then 99.

[_Counselor]

So pick your poison. Writing a wallet? Use exponentiation. EGCD is much faster, but is unsuitable for use in wallets.

Additional reading:

https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

FYI: https://gcd.cr.yp.to/papers.html
The fastest modular inversion with constant time that I know. About 4000 CPU clock cycles per inversion. Also check "software" tab for avx-optimized example sources.
It can be optimized further but at the cost of constant time implementation.

[COBRAS]

So pick your poison. Writing a wallet? Use exponentiation. EGCD is much faster, but is unsuitable for use in wallets.

Additional reading:

https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

FYI: https://gcd.cr.yp.to/papers.html
The fastest modular inversion with constant time that I know. About 4000 CPU clock cycles per inversion. Also check "software" tab for avx-optimized example sources.
It can be optimized further but at the cost of constant time implementation.

no way to use inv for downgrade priv  simple mult point  to mod inv(a) will give redult point * inv(a) bits !!!!

[NotATether]

So pick your poison. Writing a wallet? Use exponentiation. EGCD is much faster, but is unsuitable for use in wallets.

Additional reading:

https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm

FYI: https://gcd.cr.yp.to/papers.html
The fastest modular inversion with constant time that I know. About 4000 CPU clock cycles per inversion. Also check "software" tab for avx-optimized example sources.
It can be optimized further but at the cost of constant time implementation.

That's great, I'll check it out when I get some time.

I tried to implement this fast modular exponentation into C++, but it works only for small numbers. Obviously, it won't work for long numbers, because of overflow problems. This makes a dilemma. How can we even represent these numbers properly in C++?

I was thinking using libsecp256k1 for that kind of thing, but I'm at an utter loss of how to do a modinverse with a byte array as the private key, because that particular function exists but it's not exported.