Bernstein’s Factorization Method Helped Factor RSA-240 in 2020

Here are free GPU credits. First come first serve :)

A newsletter about theoretical computer science and applied mathematics. Subscribe.

We code one of the integer factorization methods used to factor a 795-bit RSA number in 2020 for our CADO-NFS clone. Jump to the bottom for Python code :)

This is part of our Practical Number Field Sieve for Programmers series:

Part 1: Discrete Logarithms and the Index Calculus Solution.

Part 2: Solving Index Calculus Equations over Integers and Finite Fields.

Part 3: Computation of Discrete Logarithms in Algebraic Number Fields.

Part 4: 2 Dimensional Lattice Basis Reduction.

Part 5: Individual Reduction Phase and Logarithm Collection.

Part 6: Lattice Sieving and Special Q Descent.

Part 7: Continued Fractions and Fast Lattice Sieving.

Part 8: Descending Multiple Integers using Bernstein’s Algorithm.

Code is available on Colab.

1.0 Paper Introduction

How To Find Small Factors of Integers (Bernstein, 2002) introduces a fast algorithm to find all small factors less than a bound B, for a given list of integers.

Bernstein’s algorithm is also called batch smoothness detection and contributed to 25% less time spent in factoring RSA-240, a 795-bit number (Boudot et al., 2020).

In our case, we desire to find multiple individual discrete logarithms for our number field sieve and this is called batch DLP or delayed-target DLP (Guillevic, 2015).

This paper uses product trees like those encountered in the Pollard P-1 factorization algorithm and implementations focus on finding smooth parts of integers (Bernstein, 2004).

2.0 Bernstein’s Algorithm

Bernstein’s factorization algorithm takes as input a list of integers of length n, each integer of bit-length y, and outputs for each integer a list of small prime factors less than a bound B (Bernstein, 2002).

The complete algorithm is presented below:

We code the Python skeleton below then the next sections present the necessary auxilliary algorithms.

#Algorithm 7.1 : factor many integers by small primes
def BernsteinFactorization(N, P):
    """
    Algorithm 7.1.
    Input: N, list of targets, P list of small primes
    Output: list of (n, primeList) pairs, 
    """
    #Step 1: base case
    if not N:
        return []
    if not P:
        return [(n, []) for n in N]
    #Step 2: multiply all n
    x = 1
    for n in N:
        x *= n
    #Step 3: product tree of P
    T = ConstructProductTree(P)
    #Step 4: find dividing primes
    P_prime = FindDividingPrimes(x, T)
    return BernsteinFactorization_Recursive(N, P_prime)


def BernsteinFactorization_Recursive(N, P_prime):
    """Recursion in Algorithm 7.1."""
    #Step 5: Base case
    if len(N) == 1:
        n = N[0]
        return [(n, [p for p in P_prime if n % p == 0])]
    mid = len(N) // 2
    M = N[:mid]
    rest = N[mid:]
    return BernsteinFactorization_Recursive(M, P_prime) + BernsteinFactorization_Recursive(rest, P_prime)

2.1 Algorithm 6.1: Finding Product Trees

Step 2 and 3 of Bernstein factorization demand us find a product tree: a binary tree of positive integers where each non-leaf vertex is a product of its two children (Bernstein, 2002).

Product trees are constructed recursively as given below:

We are provided an example of a product tree as well:

Python code resembles:

#Algorithm 6.1 : product tree
def ConstructProductTree(P):
    """
    Algorithm 6.1.
    Input: P, a list of integers
    Output: Tree given as list of tuples
    """
    if len(P) == 1:
        return P[0]
    mid = len(P) // 2
    left = ConstructProductTree(P[:mid])
    right = ConstructProductTree(P[mid:])
    leftRoot = left if isinstance(left, int) else left[2]
    rightRoot = right if isinstance(right, int) else right[2]
    return (left, right, leftRoot * rightRoot)

2.2 Algorithm 6.3: Find Cases x mod p = 0

Step 4 incorporates the intermediate products obtained in Step 2 and Step 3 to find small divisors:

Python code resembles:

#Algorithm 6.3 : find primes dividing x
def FindDividingPrimes(x, T):
    """
    Algorithm 6.3.
    Input: Integer x >= 0, Product tree T of a nonempty list of odd primes,
    Output: The list of primes p in T such that x mod p == 0.
    """
    result = []
    FindDividingPrimes_Recursive(x, T, result)
    return result

def FindDividingPrimes_Recursive(x, T, result):
    def FindTreeRoot(T):
        return T if isinstance(T, int) else T[2]
    u = FindTreeRoot(T)
    c = u.bit_length() #ceil(lg(u+1))
    d = x.bit_length() if x > 0 else 0
    if d > c + 1:
        r = TwoAdicDivision(d - c, c, u, x)
    else:
        r = x
    #Step 5: Now 0 <= r < 4u and 2^k * r ≡ x (mod u) for some k.
    if isinstance(T, int):
        if r % u == 0:
            result.append(u)
        return
    left, right, _ = T
    FindDividingPrimes_Recursive(r, left, result)
    FindDividingPrimes_Recursive(r, right, result)

2.3 Algorithm 5.3: 2-adic Division

Step 3 of Algorithm 6.3 of the previous section introduces a fast alternative to division for find a good r. It’s given below:

In Python:

#Algorithm 5.3 : 2-adic division
def TwoAdicDivision(b, c, u, x):
    """
    Algorithm 5.3.
    Input: b,c > 0, odd u < 2^c, and 0 <= x < 2^(c+b),
    Output: r < 2^(c+1) with 2^b * r ≡ x (mod u).
    """
    v = TwoAdicInverse(b, u)
    mask_b = (1 << b) - 1
    x0 = x & mask_b
    x1 = x >> b
    q = (v * x0) & mask_b
    r = x1 + ((x0 + u * q) >> b)
    return r

2.4 Algorithm 5.1: 2-adic Inverse

Step 1 of Algorithm 5.3 demands us find the 2-adic inverse. This is given in algorithm 5.1 as shown below:

In Python,

#Algorithm 5.1 : 2-adic inverse
def TwoAdicInverse(b, u):
    """
    Algorithm 5.1.
    Input: b > 0 and odd u > 0
    Output: v < 2^b with 1 + u*v ≡ 0 (mod 2^b).
    """
    if b == 1:
        return 1
    c = (b + 1) // 2                      # ceil(b/2)
    v0 = TwoAdicInverse(c, u)
    mask_c = (1 << c) - 1
    u0 = u & mask_c
    u1 = (u >> c) & mask_c
    z = (((1 + u0 * v0) >> c) + u1 * v0) & mask_c
    v = (v0 + (z * v0 << c)) & ((1 << b) - 1)
    return v

3.0 Numerical Example

The algorithm is proven to terminate correctly and (Bernstein, 2002) provides an example as shown below:

The algorithm in its entirety is given below:

"""
Implementation of Algorithm 7.1 from:
    Daniel J. Bernstein, "How to Find Small Factors of Integers"

Finds all small prime divisors of many integers simultaneously.
"""
#Algorithm 5.1 : 2-adic inverse
def TwoAdicInverse(b, u):
    """
    Algorithm 5.1.
    Input: b > 0 and odd u > 0
    Output: v < 2^b with 1 + u*v ≡ 0 (mod 2^b).
    """
    if b == 1:
        return 1
    c = (b + 1) // 2                      # ceil(b/2)
    v0 = TwoAdicInverse(c, u)
    mask_c = (1 << c) - 1
    u0 = u & mask_c
    u1 = (u >> c) & mask_c
    z = (((1 + u0 * v0) >> c) + u1 * v0) & mask_c
    v = (v0 + (z * v0 << c)) & ((1 << b) - 1)
    return v
    
#Algorithm 5.3 : 2-adic division
def TwoAdicDivision(b, c, u, x):
    """
    Algorithm 5.3.
    Input: b,c > 0, odd u < 2^c, and 0 <= x < 2^(c+b),
    Output: r < 2^(c+1) with 2^b * r ≡ x (mod u).
    """
    v = TwoAdicInverse(b, u)
    mask_b = (1 << b) - 1
    x0 = x & mask_b
    x1 = x >> b
    q = (v * x0) & mask_b
    r = x1 + ((x0 + u * q) >> b)
    return r

#Algorithm 6.1 : product tree
def ConstructProductTree(P):
    """
    Algorithm 6.1.
    Input: P, a list of integers
    Output: Tree given as list of tuples
    """
    if len(P) == 1:
        return P[0]
    mid = len(P) // 2
    left = ConstructProductTree(P[:mid])
    right = ConstructProductTree(P[mid:])
    leftRoot = left if isinstance(left, int) else left[2]
    rightRoot = right if isinstance(right, int) else right[2]
    return (left, right, leftRoot * rightRoot)
    
#Algorithm 6.3 : find primes dividing x
def FindDividingPrimes(x, T):
    """
    Algorithm 6.3.
    Input: Integer x >= 0, Product tree T of a nonempty list of odd primes,
    Output: The list of primes p in T such that x mod p == 0.
    """
    result = []
    FindDividingPrimes_Recursive(x, T, result)
    return result

def FindDividingPrimes_Recursive(x, T, result):
    def FindTreeRoot(T):
        return T if isinstance(T, int) else T[2]
    u = FindTreeRoot(T)
    c = u.bit_length() #ceil(lg(u+1))
    d = x.bit_length() if x > 0 else 0
    if d > c + 1:
        r = TwoAdicDivision(d - c, c, u, x)
    else:
        r = x
    #Step 5: Now 0 <= r < 4u and 2^k * r ≡ x (mod u) for some k.
    if isinstance(T, int):
        if r % u == 0:
            result.append(u)
        return
    left, right, _ = T
    FindDividingPrimes_Recursive(r, left, result)
    FindDividingPrimes_Recursive(r, right, result)

#Algorithm 7.1 : factor many integers by small primes
def BernsteinFactorization(N, P):
    """
    Algorithm 7.1.
    Input: N, list of targets, P list of small primes
    Output: list of (n, primeList) pairs, 
    """
    #Step 1: base case
    if not N:
        return []
    if not P:
        return [(n, []) for n in N]
    #Step 2: multiply all n
    x = 1
    for n in N:
        x *= n
    #Step 3: product tree of P
    T = ConstructProductTree(P)
    #Step 4: find dividing primes
    P_prime = FindDividingPrimes(x, T)
    return BernsteinFactorization_Recursive(N, P_prime)


def BernsteinFactorization_Recursive(N, P_prime):
    """Recursion in Algorithm 7.1."""
    #Step 5: Base case
    if len(N) == 1:
        n = N[0]
        return [(n, [p for p in P_prime if n % p == 0])]
    mid = len(N) // 2
    M = N[:mid]
    rest = N[mid:]
    return BernsteinFactorization_Recursive(M, P_prime) + BernsteinFactorization_Recursive(rest, P_prime)

    
if __name__ == "__main__":
    P = [3, 5, 7, 11, 13, 17, 19]
    N = [492, 2567, 3135, 5889]
    results = BernsteinFactorization(N, P)
    for n, primes in results:
        print(f"{n}: {primes}")

Running this code yields, as expected:

Subscribe if you made it this far!

References

Bernstein, D. (2002). How To Find Small Factors of Integers. Mathematics of Computation. PDF.

Boudot, F., Gaudry, P., Guillevic, A., Heninger, N., Thomé, E., & Zimmermann, P. (2020). Comparing The Difficulty Of Factorization And Discrete Logarithm: A 240-Digit Experiment. IACR. PDF.

Guillevic, A. (2015). Computing Individual Discrete Logarithms Faster in with the NFS-DL Algorithm. In: Iwata, T., Cheon, J. (eds)., Advances in Cryptology -- ASIACRYPT 2015. Lecture Notes in Computer Science(), vol 9452. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-48797-6_7. PDF.

Bernstein, D. (2002). How To Find Smooth Parts of Integers. CRYPTO. PDF.