RSA Optimizations: Think Twice

I will start here straight ahead with the summary:

Resist your temptations and do not implement your customised and "optimized" RSA library.

We are overwhelmed today with technical breakthrough announcements, revolutionary ideas that will change our lives, tens of alert announcements on our social feeds ... LLMs are at our disposal, available for explaining everything in a couple of sentences. Who said that we must read and understand thoroughly? We have our AI personal assistants that will provide us concise and short summary! What a power we have, let's implement and optimise something, we don't want being left behind our competition, they've just announced that they've cut down their processing time by 65% due to a smart algorithm optimisation.

Historically, the RSA key generation process has been "optimised" many times. Vendors were seeking how to generate keys more quickly, how to encrypt faster, how to decrypt faster. You know, time is money... They have tried everything: key generation based on patterns, large public keys/ small private keys.... And all these attempts have failed miserably, introducing vulnerabilities that allow you to factor RSA and derive the private key based on the public key information.

Example: Estonia is considered as a leader in the digital transformation processes. In 2017 they've faced national wide cyber crisis and were forced to default and withdraw about 750,000 identity cards because of the vulnerable smart card chips due to the ROCA vulnerability. The same Infineon RSA chip was used in many smart cards and TPMs worldwide.

It is my personal opinion that these days we are rushing and pushing the boundaries up to the limit (wait a minute, we don't like pesimists, there are no boundaries, sky is the limit!) without thorough understanding of the engineering and mathematical principles. At the end of the day, who needs a math, we have our LLM copilot, a creativity booster full of ideas, empathy and politness (she never makes me feel stupid and uneducated).

Nevertheless, I will put myself in the shoes of an old fashioned engineer for a moment, maybe somewhere, sometime it will pay off.

Mathematics has been always available for those wanting to understand how something works under the hood.

Here is the scenario:

A vendor wants to optimise the RSA public encryption/decryption process. The focus is having considerable fast decryption. Use case: large number of participants will encrypt a message for the same recipient, but that recipient must decrypt ALL of those messages, it must be fast for sure!. The vendor's solution is creating long RSA public keys used for encryption and relatively short private keys used for decryption.

Let's us have a short summary how the RSA crypto system works:

Key generation
- - - - - - - - - - -
Pick two large prime numbers: p and q

The modulus N is: 

N = p * q

Calculate Euler's phi function:

phi = (p - 1) * (q - 1)

Choose a public key e such that it is co-prime with phi:

gcd(e, phi) = 1

Calculate the private key d as:

e * d = 1 mod(phi)

Encryption
- - - - - - - - - - -

c = m ^ e mod(N)

Decryption
- - - - - - - - - - -

m = c ^ d mod(N)        

Statement: If the private key is small enough, i.e. d < 1/3 * N^(1/4), then we can efficiently factor the modulus N and derive the private key d from the public key information (e, N).

The above statement is known as the Wiener's attack. The attack is simple. Here is how it goes.

e * d = 1 mod(phi)

e * d = k * phi + 1

e * d = k * (p - 1) * (q - 1) + 1 

(e * d)/(d * p * q)= (k * (p - 1) * (q - 1) + 1) / (d * p * q)

e / (p * q) = (k / d) * ((p - 1) * (q - 1) + 1 / k)/(p * q)

e / (p * q) = (k / d) * (p * q - (p + q) + 1 + 1 / k) / (p * q)

e / N = (k / d) * (1 - delta)

where:

delta = ((p + q) - 1 - 1 / k) / (p * q)
        

The value of delta is very small. That means e/N is slightly smaller than k/d. We don't know what k/d is, but we can calculate all e/N approximations and assume each approximation to be k/d. For that purpose we can use the method of Continued fraction.

The example from the Wiener's attack Wikipedia article gives us:

e/N = 17993/90581 = 
= 1/(5 + 1/(29 + 1/(4 + .... )))        

Every next approximation is smaller than the previous because we introduce additional term in the denominator. These approximations or bounds are:

1 /5
29 / 146
117 / 589
...        

So, e/N is slightly smaller than 1/5 or e/N is slightly smaller than 29/146 and so on. In other words, k/d = 0 or k/d = 1/5 or k/d = 29/146 ...

We can calculate the fractions and bounds in Python as:

def calculate_bound(fractions):

    bound = (1, fractions[-1])
    
    for f in fractions[-2:0:-1]:
    
        bound = (bound[1], bound[1] * f + bound[0])
            
    return bound

def continued_fractions(a,b):

    assert a < b
    
    fractions = [0]
    
    bounds = []
    
    while a != 0:
    
        q = b // a
        
        fractions.append(q)
        
        bounds.append(calculate_bound(fractions))
        
        r = b % a
        
        a, b = r, a
    
    return (fractions, bounds)        

Once we have the bounds we will iterate each of them until RSA conditions are satisfied:

k / d = bound

# we've assumed values for k and d at this point
# for 1/5, k = 1 and d = 5

e * d = 1 + k * phi

phi = (e * d - 1) // k        

At this point, we known the actual values for e, N and assumed value for phi calculated from the assumed values for k and d.

A quick remainder about the quadratic equation solutions:

a * x^2 + b * x + c = 0

D = sqrt(b^2 - 4 * a * c)

x1 = (-b + D) / (2 * a)

x2 = (-b - D) / (2 * a)        

We want to find p and q with the available and assumed information at our disposal. If we know p and q, we can calculate N' and compare if it matches the real, known N. If they match we are finished, we can calculate the private key d. If not, we shall try the next bound (k/d).

What are the solutions for quadratic equation like the following one?

x^2 - (p + q) * x + p * q = 0

D = sqrt ((p + q)^2 - 4 * p * q)

= sqrt(p^2 + 2 * p * q + q^2 - 4 * p * q)

= sqrt((p^2 - 2 * p * q + q^2)) = 

= sqrt((p - q)^2) 

= p - q

x1 = (p + q + p - q) / 2 = p

x2 = (p + q - p + q) / 2 = q        

Note that:

p + q = p * q - (p - 1) * (q - 1) + 1

p + q = p * q - p * q + p + q -1 + 1 = p + q

=>

N = p * q

phi = (p - 1) * (q - 1)

p + q = N - phi + 1        

We need to solve this equation that yields p and q:

x^2 - (N - phi + 1) * x + N = 0        

For that purpose we will use one more beautiful numerical math method, the Newton's method for finding polynomial roots.

In Python:

def derivative(x, a2, a1, a0):

    return a2 * 2 * x + a1

def value(x, a2, a1, a0):

    return a2 * x ** 2 + a1 * x + a0
    
def newton(x, a2, a1, a0):

    # derivative = (y1 - y0) / (x1 - x0)
    # y1 - y0 = derivative * (x1-x0)
    # x1 - x0 = (y1 - y0) / derivative 
    # x1 = x0 - y0/derivative
    
    delta = 1
    
    with localcontext() as ctx:
    
        ctx.prec = 2048
        
        x = Decimal(x)
        
        a2 = Decimal(a2)
        
        a1 = Decimal(a1)
        
        a0 = Decimal(a0)
    
        while delta > 0.01:
            
            x_1 = x - value(x, a2, a1, a0)/derivative(x, a2, a1, a0)
            
            delta = abs(x_1 -x)
            
            x = x_1
            
        p = round(x)
        
    return p        

The bounds iteration logic, the phi approximation, the N calculation and comparison are:

for bound in bounds:

    d = bound[1]
    
    k = bound[0]
    
    phi = (e * d - 1)//k
    
    a2 = 1
    
    a1 = phi - n - 1
    
    a0 = n
    
    p = newton(x, a2, a1, a0)
    
    q = n // p
    
    if p * q == n:
    
        break        

Let us see how it works in practice for large primes and modulus. The complete code for the Wiener's attack is:

import math

from decimal import *

def calculate_bound(fractions):

    bound = (1, fractions[-1])
    
    for f in fractions[-2:0:-1]:
    
        bound = (bound[1], bound[1] * f + bound[0])
            
    return bound

def continued_fractions(a,b):

    assert a < b
    
    fractions = [0]
    
    bounds = []
    
    while a!=0:
    
        q = b // a
        
        fractions.append(q)
        
        bounds.append(calculate_bound(fractions))
        
        r = b % a
        
        a, b = r, a
    
    return (fractions,bounds)
        

def derivative(x, a2, a1, a0):

    return a2 * 2 * x + a1

def value(x, a2, a1, a0):

    return a2 * x ** 2 + a1 * x + a0
    
def newton(x, a2, a1, a0):

    # derivative = (y1 - y0) / (x1 - x0)
    # y1 - y0 = derivative * (x1-x0)
    # x1 - x0 = (y1 - y0) / derivative 
    # x1 = x0 - y0/derivative
    
    delta = 1
    
    with localcontext() as ctx:
    
        ctx.prec = 2048
        
        x = Decimal(x)
        
        a2 = Decimal(a2)
        
        a1 = Decimal(a1)
        
        a0 = Decimal(a0)
    
        while delta > 0.01:
            
            x_1 = x - value(x, a2, a1, a0)/derivative(x, a2, a1, a0)
            
            delta = abs(x_1 -x)
            
            x = x_1
            
        p = round(x)
        
    return p

n =int('0xb8af3d3afb893a602de4afe2a29d7615075d1e570f8bad8ebbe9b5b9076594cf06b6e7b30905b6420e950043380ea746f0a14dae34469aa723e946e484a58bcd92d1039105871ffd63ffe64534b7d7f8d84b4a569723f7a833e6daf5e182d658655f739a4e37bd9f4a44aff6ca0255cda5313c3048f56eed5b21dc8d88bf5a8f8379eac83d8523e484fa6ae8dbcb239e65d3777829a6903d779cd2498b255fcf275e5f49471f35992435ee7cade98c8e82a8beb5ce1749349caa16759afc4e799edb12d299374d748a9e3c82e1cc983cdf9daec0a2739dadcc0982c1e7e492139cbff18c5d44529407edfd8e75743d2f51ce2b58573fea6fbd4fe25154b9964d', 16)
e = int('0x9ab58dbc8049b574c361573955f08ea69f97ecf37400f9626d8f5ac55ca087165ce5e1f459ef6fa5f158cc8e75cb400a7473e89dd38922ead221b33bc33d6d716fb0e4e127b0fc18a197daf856a7062b49fba7a86e3a138956af04f481b7a7d481994aeebc2672e500f3f6d8c581268c2cfad4845158f79c2ef28f242f4fa8f6e573b8723a752d96169c9d885ada59cdeb6dbe932de86a019a7e8fc8aeb07748cfb272bd36d94fe83351252187c2e0bc58bb7a0a0af154b63397e6c68af4314601e29b07caed301b6831cf34caa579eb42a8c8bf69898d04b495174b5d7de0f20cf2b8fc55ed35c6ad157d3e7009f16d6b61786ee40583850e67af13e9d25be3', 16)

fractions,bounds = continued_fractions(e, n)

print(fractions)

print(bounds)

x = 1

match = False

for bound in bounds:

    d = bound[1]
    
    k = bound[0]
    
    phi = (e * d - 1)//k
    
    a2 = 1
    
    a1 = phi - n - 1
    
    a0 = n
    
    p = newton(x, a2, a1, a0)
    
    q = n // p
    
    if p * q == n:

        match = True

        break
        
if match:

    phi = (p - 1) * (q - 1)

    d = pow(e, -1 , phi)

    print('p:', p)

    print('q:', q)

    print('d:', d)        

The RSA factorisation result reveals:

p: 134669927709128070756424801419548805501808076912262801434800605920827612464368906595661348393409080650056282638489243851059781971132159889896843018381187994628859917822755789986092763460463295405651440816348815008245093856412009397970192375577360567873185141159375280522236909060526334123001733587717969177157
q: 173121513995818161102245832683211062320154182361182024148671930683926870711405647497524667705258311163551127156955342410807842326402024536548989691926348678692995530897791794818646241971728281417961389048493180792474943296919266335463768525410560161755731138916915767148609523790355387780727531897114371948649
d: 79434351637397000170240219617391501050474168352481334243649813782018808904459        

Wow, RSA number factorisation is really easy and fast under these circumstances. The math used is not that complex too. All of this reminds me how important are the well established engineering and cyber security principles that rely on thorough analysis, critical observation and evaluation.

But that is just me :-)

Marjan