Elliptic curves and how to use them (part four)

Go to part 3

One of the basic parameters of processors is the word length of a single arithmetic operation. This value can range from 8 bits for simple micro-controllers through 32 bits for 80386 derivatives to the most popular 64-bit length today. In addition, processors may have extended instruction sets (such as MMX, SSE, etc.) that can extend the processor’s capabilities to 128-bit (or larger) words.
The word lengths provided are therefore too short for the implementation of cryptographic algorithms. For RSA, the suggested (in terms of security) key length is today (year 2019) 4096 bits, which is many times more than standard processors provide. For elliptic curves, the length of the prime number in the suggested curve is 521 bits (no, it’s not an error: $2^{521}-1$ is one of Mersenne prime numbers).
So, if we want to create an implementation of elliptic curves, we must have implementations of large integer numbers. We can use an existing one (mathematical libraries usually implement not only ‘large numbers’ but also the elliptic curves themselves) or write it yourself (it is enough to implement addition, subtraction, multiplication, modulo and division). The presented examples are written in Python, which has built-in ‘large numbers’ arithmetic, but when writing an implementation in a typical embedded we do not have this comfort – we have to implement everything in C / C ++.
Let’s check the performance of the implementation presented in the previous part, but let’s use the more demanding curve:

#3.3.  521-Bit Random ECP Group
print "ECPoint 521-Bit"
p = (2**521)-1
tv = FieldElement(p, 1)
neutral = ECPoint(p, -3, 0x0051953EB9618E1C9A1F929A21A0B68540EEA2DA725B99B315F3B8B489918EF109E156193951EC7E937B1652C0BD3BB1BF073573DF883D2C34F1EF451FD46B503F00)
gx = FieldElement.factoryMethod(tv, 0x00C6858E06B70404E9CD9E3ECB662395B4429C648139053FB521F828AF606B4D3DBAA14B5E77EFE75928FE1DC127A2FFA8DE3348B3C1856A429BF97E7E31C2E5BD66)
gy = FieldElement.factoryMethod(tv, 0x011839296A789A3BC0045C8A5FB42C7D1BD998F54449579B446817AFBD17273E662C97EE72995EF42640C550B9013FAD0761353C7086A272C24088BE94769FD16650)
g = ECPoint.factoryMethod(neutral, gx, gy)

#8.3.  521-Bit Random ECP Group
i = 0x0037ADE9319A89F4DABDB3EF411AACCCA5123C61ACAB57B5393DCE47608172A095AA85A30FE1C2952C6771D937BA9777F5957B2639BAB072462F68C27A57382D4A52
r = 0x0145BA99A847AF43793FDD0E872E7CDFA16BE30FDC780F97BCCC3F078380201E9C677D600B343757A3BDBF2A3163E4C2F869CCA7458AA4A4EFFC311F5CB151685EB9
start = monotonic.time.time()
gi = g*i
stop = monotonic.time.time()
print "Time: %s"%(stop-start)
start = monotonic.time.time()
gr = g*r
stop = monotonic.time.time()

ECPoint 521-Bit
Time: 0.146764039993
Time: 0.149250984192

So, one point multiplication on one Intel G4560 3.50GHz core amounted to over 400 ms, and at the time of establishing the connection you need to perform at least 4 multiplications (key negotiation and authentication), i.e. mathematical calculations alone would take almost 2 seconds. For processors whose speeds are measured not in GHz, but in MHz, the results would be measured in minutes, which is completely unacceptable. So we can try to optimize several implementation elements:

large numbers arithmetic,
operations in the field of $F_p$ ,
operations in the EC group ( $F_p$ ),
EC point multiplication ( $F_p$ ).

We can optimize the arithmetic of large numbers by standard in terms of computational complexity, conditional jumps, data copy elimination, etc. In our examples we assume for the purposes of our examples that in Python this implementation is well optimized.
Operations in the field of $F_p$ use the arithmetic of large numbers with additional requirements, which results in the need to use operations: modulo and element inversion algorithm. Division and modulo are the most demanding operations for the processor (division is at least 30 times slower than adding). In the implementation of part 3, each operation required the use of modulo. In addition, one could replace modulo operation by comparing to the number p and subtraction, but the conditional instruction would break the instruction piping in the processor, so that despite the prediction of jumps, the gain could not be large. For multiplication, however, there is a way to get rid of the expensive modulo operation – Montgomery’s reduction algorithm.
To be able to use the Montgomery algorithm, we must first calculate the residual $R=2^k$ which is the smallest approximation greater than p ( $2^{k-1} \; \leq \; p \; \lt \; 2^k$ ).
The operation of converting the number to a residuum form looks like this:

$a' \equiv a\ast R \;(mod\; p)$

The operation of the return conversion from the residuum form looks like this:

$a \equiv a' \ast R^{-1}\;(mod\;p)$

The operation of adding numbers in the form of residuum looks identical to the usual addition in $F_p$ :

$c' \equiv a' + b'\;(mod\;p)$

because:

$a'+b' \equiv a \ast R + b \ast R \equiv (a+b) \ast R \equiv c \ast R \equiv c' \;(mod\;p)$

The operation of multiplying numbers in the form of residuum looks as follows:

$c' \equiv a \ast b \ast R^{-1} \;(mod\;p)$

because:

$a'\ast b' \equiv a \ast R \ast b \ast R \equiv (a \ast b) \ast R^2 \equiv c \ast R^2 \equiv c' \ast R \;(mod\;p)$

However, the above multiplication formula is much less optimal than multiplying in $F_p$ . Therefore, the Montgomery algorithm is used to optimize the multiplication of numbers in the form of residuum, for which it is necessary to calculate $inv_p = \frac{R \ast R^{-1} - 1}{ p }$

The algorithm looks like this:

$T=a' \ast b'$
$m= ((T \; mod \;R) \ast inv_p \; )mod \; R$
$t=\frac{T+m \ast p}{R}$
If $t\geq p$
then $c'=t-p$
else $c'=t$

In the above algorithm we have 3 demanding operations: twice modulo and once dividing. However, the basis for modulo is the number $R=2^k$ , so we can simplify these operations:

$c \; mod \; R=c\wedge(R^{-1})$ [bit mask]
$\frac{c}{R}=c\gg k$ [shifting bits to the right]

So our implementation of $F_p$ arithmetic using the Montgomery algorithm might look like this:

#!/usr/bin/python
#...
class FieldElement:
    def __init__(self, p, rvalue, Rmask, k, invR, invp):
        self.p = p
        self.rvalue = rvalue
        self.Rmask = Rmask
        self.k = k
        self.invR = invR
        self.invp = invp
        
    def __str__(self):
        if self.rvalue < 0:
            return "-0x%x"%(self.convertFromResiduum()*(-1))
        return "0x%x"%self.convertFromResiduum()
    
    def __add__(self, other):
        if type(other) == int or type(other) == long:
            return self.add((other << self.k)%p)
        else:
            return self.add(other.rvalue)
    
    def add(self, rvalue):
        rvalue += self.rvalue
        if rvalue >= self.p:
            rvalue -= self.p
        return FieldElement(self.p, rvalue, self.Rmask, self.k, self.invR, self.invp)
    
    def __sub__(self, other):
        if type(other) == int or type(other) == long:
            return self.sub((other << self.k)%p)
        else:
            return self.sub(other.rvalue)
    
    def sub(self, rvalue):
        res = self.rvalue
        if rvalue > res:
            res += self.p
        res -= rvalue
        return FieldElement(self.p, res, self.Rmask, self.k, self.invR, self.invp)
    
    def __mul__(self, other):
        if type(other) == int or type(other) == long:
            return self.mul((other << self.k)%p)
        else:
            return self.mul(other.rvalue)
    
    def mul(self, rvalue):
        T = rvalue * self.rvalue
        #m = ( (T & self.Rmask) * self.invp) & self.Rmask
        m = ( T * self.invp) & self.Rmask
        t = (T + m * self.p) >> self.k
        if t >= self.p:
            rvalue = t - self.p
        else:
            rvalue = t
        return FieldElement(self.p, rvalue, self.Rmask, self.k, self.invR, self.invp)
    
    def __eq__(self, other):
        return (self.rvalue == other.rvalue)
    
    def __ne__(self, other):
        return (self.rvalue != other.rvalue)

    def convertFromResiduum(self):
        return ((self.rvalue*self.invR)%self.p)
    
    @staticmethod
    def createMethod(value, p):
        R = 1
        k = 0
        while R <= p:
            R <<= 1
            k += 1
        invR = invers(R, p)
        invp = ((R * invR) - 1) / p
        rvalue = (value * R)%p
        return FieldElement(p, rvalue, R - 1, k, invR, invp)
    
    @staticmethod
    def factoryMethod(otherFieldElement, value):
        rvalue = (value << otherFieldElement.k)%otherFieldElement.p
        return FieldElement(otherFieldElement.p, rvalue, otherFieldElement.Rmask, otherFieldElement.k, otherFieldElement.invR, otherFieldElement.invp)
#...

ECPoint 521-Bit
Time: 0.161954164505
Time: 0.163055896759

And … it’s worse 🙂
My first version of this code was even slower. Therefore, the above optimization makes sense only in implementations of large numbers in which operations on bits (shifting, masking) are efficient – it seems that in Python these operations are not. So you have to take my word for it that in C / C ++ implementation the above algorithm noticeably speeds up operations on the $F_p$ field.
In the next part there will be presented optimization, which will significantly accelerate the operations on the curves: Jacobian coordinates.

ec_part4 (full code without optimization)
ec_part4m (full code with Montgomery multiplication)

Post Views: 66

Leave a Reply Cancel reply