Factorization

12. Factorization#

12.1. Naive Factorization#

How can we find the prime factorization of a number $N$? If one knows all the primes less than or equal to $N$, this is not so bad: we test each prime. If it divides $N$ we see how many times it does, and save the both the prime and exponent. We then keep dividing until we get $1$. Here is an implementation of this process:

def naive_factor(n):
    """
    Given n returns the factiorization of n as a list of tuples (p, k),
    where p is prime, and k is the power of p dividing n.
    """
    factorization = []  # result
    p = 1  # initialized (not used)
    while n > 1:
        p = next_prime(p)  # prime
        k = 0  # exponent
        while n % p == 0:
            k += 1
            n /= p
        if k != 0:
            factorization.append((p, k))
    return factorization

Note that we’ve used Sage’s next_prime. If not using Sage, we would need to use some primality test, which can be time consuming, to find the next prime. But note that this is a case in which a probabilistic test is fine: it might slow things down in the unlikely event that it gives a false positive, but it will not produce a fake prime factor.

Indeed, we would be testing element after element, so if some number $q$ is composite and the test has told us that it is prime, then $q$ will not divide the number, since it must be made of smaller prime factors that do not divide the number anymore at this point.

Thus, we would have an unnecessary division, but no error is introduced in the result.

Let’s test it in a couple of simple examples.

naive_factor(6)

[(2, 1), (3, 1)]

naive_factor(17)

[(17, 1)]

Now, let’s run a more extensive test:

lower_bound = 2
upper_bound = 10^6

number_of_tests = 10^2
for _ in range(number_of_tests):
    n = randint(lower_bound, upper_bound + 1)
    factorization = naive_factor(n)
    m = prod(p^k for (p, k) in factorization)
    if n != m:
        print(f"Failed for {n = }.  Factorization given: {factorization}.")
else:
    print("It works!")

It works!

It seems to work, but if you run it you will likely see that it is slow. Therefore, we need more efficient ways to factor large integers, specially if they only have large prime factors.

12.2. Pollard’s $p-1$ Factorization Algorithm#

The following idea and resulting algorithm were introduced by J. Pollard in 1974.

Suppose you want to factor $N = pq$, where $p$ and $q$ are very large distinct primes, as used with the RSA encryption. Suppose that somehow we find an integer $L$ such that:

$L = (p - 1)i$, for some integer $i$, and so $p - 1$ divides $L$;
$L = (q - 1)j + k$, for integers $j$ and $k$, with $k \in \{ 1, 2, \ldots, (q-2) \}$, and so $q-1$ does not divide $L$.

Then, if $\gcd(a, N) = 1$, i.e., $a$ not divisible by either $p$ or $q$, then, by Fermat's Little Theorem (or Euler's Theorem):

\[a^L = a^{(p-1)i} = \left( a^{p-1} \right)^i \equiv 1^i = 1 \pmod{p}\]

and

\[a^L = a^{(q-1)j + k} = \left( a^{q-1} \right)^j \cdot a^k \equiv 1^j \cdot a^k = a^k \pmod{q}.\]

For a randomly chosen $a$ relatively prime to $pq$, we have that the probability that $a^k \not\equiv 1 \pmod{q}$ is high, and in this case we have that $p \mid a^L - 1$ and $q \nmid a^L - 1$, so $\gcd(a^L - 1, N) = \gcd(a^L - 1, pq) = p$. Using the Euclidean algorithm we can quickly compute this GCD and find $p$. Dividing $N$ by $p$ we find $q$ and have factored $N$.

So, the question is really how do we find this $L$ (with the properties above)? Here is Pollard’s observation: if $p - 1$ is a product of small primes only, then $p-1$ divides $n!$ for some $n$ “not too large”. Then, if $\gcd(a, p) = 1$, we have that

\[a^{n!} = \left( a^{p-1} \right)^{n!/(p-1)} \equiv 1^{n!/(p-1)} = 1 \pmod{p}\]

(noting that in this case $n!/(p-1)$ is an integer by assumption).

Hoping that $a^{n!} \not\equiv 1 \pmod{q}$, which can happen with a reasonable likelihood, as observed above, we can reduce $a^{n!}$ modulo $N$, obtaining, say, $r$, and then compute $\gcd(r, N)$. If it is neither $1$ nor $N$, then this GCD is $p$, one the prime factors, and we can factor $N$.

Here is an example: consider $p = 76801$ and $q = 55837$.

p, q = 76801, 55837
is_prime(p), is_prime(q)

(True, True)

Let’s factor $p-1$ and $q-1$:

factor(p - 1)

2^10 * 3 * 5^2

factor(q - 1)

2^2 * 3^3 * 11 * 47

We need $n$ such that $n!$ is divisible by $p - 1$, so we need $n \geq 12$:

factorial(11) % (p - 1)

factorial(12) % (p - 1)

(So, indeed $n$ is relatively small!)

On the other hand, $q - 1$ does not divide $12!$:

factorial(12) % (q - 1)

So, we see that:

Mod(2, p)^factorial(12)

Mod(2, q)^factorial(12)

Here is a possible factorization procedure that uses this idea:

Algorithm 12.1 (Pollard’s $p-1$ Factorization)

Given a composite integer $N$, we want to find a proper factor (i.e., a factor different from $1$ and $N$) of $N$.

Choose a bound $B$ for $n$ (as in the $n!$ above).
Initialize $a \leftarrow 2$ (in $\mathbb{Z}/N\mathbb{Z}$) and $b \leftarrow a$.
For $j$ in $\{ 2, 3, \ldots, B \}$:
1. Set $b$ as the reduction modulo $N$ of $b^j$.
2. Compute $d = \gcd(b-1, N)$.
3. If $d \neq 1, N$, then it is a factor, and we return it (and stop). If not, go to the next iteration of the loop.
If the loop reaches $B$ with not success, increase $a$, avoiding perfect powers, e.g., going to the next prime, and repeat (from 2).

First, observe that in step 3, we are computing $a^{j!}$, as the values of $b$ are:

\[b \leftarrow a = a^{1!}, \quad b \leftarrow b^2 = a^2 = a^{2!}, \quad b \leftarrow b^3 = (a^2)^3 = a^{3!}, \quad b \leftarrow b^4 = (a^{3!})^4 = a^{4!}, \quad \text{etc.}\]

One can also try to reduce the number of operations: note that if $(p-1)$ divides $n!$, then it divides $m!$ for any $m > n$ (as $n!$ is then a factor of $m!$). So, instead of computing $\gcd(b-1, N)$ at every step, we can do it at every $k$ steps, reducing the number of operations. So, we would only compute $\gcd(a^{k!} - 1, N)$, $\gcd(a^{(2k)!} - 1, N)$, $\gcd(a^{(3k)!} - 1, N)$, etc.

The downside of this approach is that increases the chances that both $p-1$ and $q-1$ divide $a^{(rk)!}$, which we need to avoid. But, if $k$ is not too large, this is not very likely and it might be worth the risk.

Homework

You will implement this algorithm in your homework.

12.2.1. Number of Operations#

If we perform the powers in the algorithm using Fast Powering, the products in the loop will require, at worst,

\[\begin{split}\begin{align*} (2 \log_2(2) + 2) &+ (2 \log_2(3) + 2) + \cdots + (2 \log_2(B) + 2) \\ &= 2 (\log_2(2) + \log_2(3) + \cdots + \log_2(B)) + 2B - 2 \\ &= 2 \log_2(B!) + 2B - 2 \end{align*}\end{split}\]

products. So, indeed, this $B$ cannot be large!

We also need $B$ computations of GCDs, and this will require, in the original algorithm, at most

\[(B-1) \cdot (2 \log_2(N) + 2)\]

long divisions. If we skip in step-sizes of $k$, as explained above, this number would be divided by $k$. (But note that $k$ should not be large!)

So, although this is a great improvement in factorization of large numbers (with a prime factor $p$ such that $p-1$ has only small prime factors), it is still a very time consuming algorithm.

12.2.2. Example#

Let’s try to use this method to factor $N = 4{,}288{,}337{,}437$:

N = 4288337437

We initialize $a$ and $b$ as $2$ modulo $N$:

a = b = Mod(2, N)

Now, comes the loop. Here, we will not explicitly choose $B$, but just see how many iterations we need until we find a factor.

We start with $j=2$, and set $b$ to the residue of $b^2$ modulo $N$.

j = 2
b ^= j  # or b = b^j

Important

Note that we use Mod(b, N)^j to efficiently compute the power $b^j$ and reduce it modulo $N$. Do not use Mod(b^j, N), as it is a lot less efficient.

Now, we compute the GCD of $b-1$ and $N$, first converting $b$ back to $\mathbb{Z}$, and check if we get something different from $1$ and $N$. (Note we need to convert $b$ to an integer.)

d = gcd(ZZ(b) - 1, N)
d

We didn’t, so we repeat:

j = 3
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 4
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 5
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 6
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 7
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 8
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 9
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 10
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 11
b ^= j
d = gcd(ZZ(b) - 1, N)
d

j = 12
b ^= j
d = gcd(ZZ(b) - 1, N)
d

Since we’ve got $768011$, this means that $76801$ is a factor of $N$:

divmod(N, 76801)

(55837, 0)

Indeed, we have that $N = 76801 \cdot 55837$, and we have factored $N$.

N == 76801 * 55837

True

Note that in principle we do not know if this is the prime factorization of $N$, but if we know that $N$ is the modulus for some RSA cryptosystem, then it must be.

Let’s now illustrate the method if we compute the GCD only every $k$ steps. Let’s use $k=3$ here.

We initialize it the same way:

a = b = Mod(2, N)

But now we compute $k$ consecutive powers, without computing the GCD:

j = 2
b ^= j

j = 3
b ^= j

j = 4
b ^= j

Then, we compute the GCD.

d = gcd(ZZ(b) - 1, N)
d

Since we got $1$, we repeat: three consecutive powers, and then the GCD:

j = 5
b ^= j

j = 6
b ^= j

j = 7
b ^= j

d = gcd(ZZ(b) - 1, N)
d

Repeat it:

j = 8
b ^= j

j = 9
b ^= j

j = 10
b ^= j

d = gcd(ZZ(b) - 1, N)
d

j = 11
b ^= j

j = 12
b ^= j

j = 13
b ^= j

d = gcd(ZZ(b) - 1, N)
d

Again, we found the factor $76801$.

Note that $p=76801$ is such that $p-1$ has only small factors:

factor(76800)

2^10 * 3 * 5^2

On the other hand, if $q = 55837$, then $q-1$ does have some “large-ish” factors:

factor(55836)

2^2 * 3^3 * 11 * 47

12.3. Choosing $p$ and $q$ for the RSA Cryptosystem#

Firstly, we must note:

Important

Since Algorithm 12.1 is a lot more efficient than the naive factorization, when producing $N = pq$ for the RSA cryptosystem, we need to choose $p$ and $q$ such $p-1$ and $q-1$ have large prime factors!

Let’s introduce some useful terminology:

Definition 12.1 ($B$-Smooth Numbers)

For some $B >0$, a number $n$ is called $B$-smooth if all its prime factors are less than or equal to $B$.

So, with the RSA, we don’t want $p-1$ and $q-1$ to $B$-smooth for any relatively small value of $B$. But, if we choose a random large prime $p$, how likely it is that $p-1$, or equivalently, $(p-1)/2$, will be $B$-smooth for a given $B$?

The number $(p-1)/2$, for some random prime $p$, is “more ore less” like a random number of size about $p/2$. (We use $(p-1)/2$ instead of $p-1$, since the latter is always even, so not really a random integer.) Thus, the question becomes:

Question

Given $B$, how likely is that $n$, within a certain bound, is $B$-smooth?

Or, alternatively:

Question

Given a certain range, how large must $B$ be so that a random element $n$ in that range has a “good chance” of being $B$-smooth?

These questions are in fact relevant to all modern factorization methods and we will come back to them in the Smooth Numbers section below.

12.4. Factorization via Difference of Squares#

We will now introduce another more efficient factorization algorithm. We start with the following definition:

Definition 12.2 (Perfect Square)

An integer $x$ is a perfect square if there is another integer $y$ such that $x = y^2$.

Remember from Algebra the factorization

\[x^2 - y^2 = (x-y)(x+y).\]

We use this idea to try to factor a large integer $N$: we can try to find some integer $b$ such that $N + b^2$ is a perfect square, i.e., there is another integer $a$ such that $N + b^2 = a^2$. (Note that this implies that $a > b$.) If we do, then

\[N = a^2 - b^2 = (a-b)(a+b).\]

If $a > b + 1$, then $a - b$ is a proper factor!

We already know that $a > b$, so $a \geq b+1$, but in fact we cannot have that $a = b + 1$, as in that case we must have that

\[\begin{split}\begin{align*} a - b &= 1, \\ a + b &= N. \end{align*}\end{split}\]

Solving this system we get $a = (N+1)/2$ and $b=(N-1)/2$. But this then implies that

\[N = (a-b)(a+b) = \frac{N^2 - 1}{4} \quad \Longrightarrow N^2 - 4N - 1 = 0.\]

Using the quadratic formula, this means that

\[N = \frac{4 \pm \sqrt{20}}{2} = 2 \pm \sqrt{5},\]

which is not an integer. Therefore, we must have that $a > b+1$, and $a-b$ would indeed be a proper factor.

Thus, we can try to take some random values for $b$ and then check if $N + b^2$ is a perfect square. (But how does one check if an integer is a perfect square efficiently?) Unfortunately this method takes too long for large $N$.

But here is an idea that can help: choose random $b$ and $k$ and see if $kN + b^2$ is a perfect square (again, one might wonder how to do that), i.e., $kN + b^2 = a^2$, for some integer $a$. This means that $kN = a^2 - b^2 = (a-b)(a+b)$. Often, either $a-b$ or $a+b$ contains some factor of $N$, and so we compute $\gcd(N, a - b)$ and $\gcd(N, a + b)$, hoping to find a proper factor.

Note that $kN = a^2 - b^2$ means that $a^2 \equiv b^2 \pmod{N}$, or $a^2 = b^2$ in $\mathbb{Z}/N\mathbb{Z}$. Moreover, since $a^2 = kN + b^2 > N$, we must have that $a > \lfloor \sqrt{N} \rfloor$. So, we want

(12.1)#\[a^2 = b^2 \quad \text{in $\mathbb{Z}/N\mathbb{Z}$, with $a > \lfloor \sqrt{N} \rfloor$.}\]

But still, finding $a$ and $b$ from (12.1) will likely take too long, so we need a method to make it more efficient.

Here is the general procedure to find $a$ and $b$ as in (12.1):

Relation Building: Find “many” $a_1, a_2, \ldots, a_r > \lfloor \sqrt{N} \rfloor$ such that the reduction modulo $N$ of $a_i^2$, which we shall denote by $c_i$, is B-smooth for some relatively small $B$.
Elimination: Take a product of $c_{i_1} c_{i_2} \cdots c_{i_s}$ of some of the $c_i$’s such that every prime in the product appears an even number of times, i.e., with an even power in the prime factorization, so that $c_{i_1} c_{i_2} \cdots c_{i_s} = b^2$ for some integer $b$.
GCD Computation: Let $a = a_{i_1} a_{i_2} \cdots a_{i_s}$ (with the same indices, i.e., we take the products of the $a_i$’s corresponding to the $c_i$’s we’ve used in the Elimination). Then

\[a^2 = a_{i_1}^2 a_{i_2}^2 \cdots a_{i_s}^2 \equiv c_{i_1} c_{i_2} \cdots c_{i_s} = b^2 \pmod{N},\]

so we compute $\gcd(N, a - b)$. If this GCD is neither $1$ nor $N$, we’ve found a factor of $N$.

Note

The fact that $a_i > \lfloor \sqrt{N} \rfloor$ implies that $a_i^2 > N$, so when we reduce $a_i^2$ modulo $N$ (to get a number between $0$ and $N-1$), we are getting a different integer than $a_i^2$ itself. This is important for this method! If not, in the Elimination process, each $c_i$ would automatically always have even powers for each prime divisor, making the process futile.

Important

In the Elimination process, we will need to factor each $c_i$, and that is why it is important that they are all $B$-smooth for some small $B$, so that these can be quickly factored.

Now, if we do find $a$ and $b$ such that $a^2 \equiv b^2 \pmod{N}$, how likely is it that $\gcd(N, a-b)$ gives a proper factor of $N$?

The worst case is when $N=pq$, with $p$ and $q$ distinct prime factors (so, $N$ has only two proper factors). Then,

\[(a-b)(a+b) = a^2 - b^2 = kN = kpq,\]

so $p \mid (a-b)$ or $p \mid (a+b)$ (and possibly dividing both), with (about) the same probability and similarly for $q$. So, there is about a $50\%$ chance that, in this case, $p \mid (a-b)$ and $q \nmid (a-b)$. And note that if that is the case, then $q \mid (a+b)$, and quite likely $p \nmid (a+b)$, since for the latter to happen, we would need that $p \mid k$, and $k$ should not be that large. So, we do not need to compute both $\gcd(N, a+b)$ and $\gcd(N, a-b)$: quite likely one gives us a factor if and only if the other does as well. Since $a-b$ is smaller, we use it instead.

Note that this $50\%$ change of finding a factor when (12.1) is satisfied is quite good, which means we would likely not need to check many pairs $(a, b)$.

12.4.1. Example#

We will discuss the Relation Building process in a future section. The GCD Computation is straight forward. So, here we will illustrate how the Elimination process works, through an example.

Suppose that the Relation Building process gave us the following $11$-smooth $c_i$’s:

\[\begin{split}\begin{align*} c_1 &= 3 \cdot 5^2 \cdot 7^3 \cdot 11, \\ c_2 &= 3 \cdot 5^5 \cdot \phantom{7^3} \cdot 11, \\ c_3 &= \phantom {3 \cdot} 5\phantom{^5} \cdot 7\phantom{^3} \cdot 11^3. \end{align*}\end{split}\]

So, we want to make products of these $c_i$’s that produce an integer that is a perfect square. That is the same as to say that all primes in their prime factorizations appear with an even power: as

\[x = p_1^{2r_1} p_2^{2r_2} \cdots p_k^{2 r_k} \quad \Longleftrightarrow \quad x = \left( p_1^{r_1} p_2^{r_2} \cdots p_k^{r_k} \right)^2.\]

Now, a product of the $c_i$’s can be written as

\[c_1^{v_1} c_2^{v_2} c_3^{v_3}, \qquad \text{where $v_i$ is either $0$ or $1$.}\]

If $v_i = 0$, it means that $c_i$ is not a part of the product (it is left out), and if $v_i$ is $1$, then $c_i$ is a factor in the product. So, collecting the powers, we have that

\[c_1^{v_1} c_2^{v_2} c_3^{v_3} = 3^{v_1 + v_2} \cdot 5^{2v_1 + 5 v_2 + v_3} \cdot 7^{3v_1 + v_3} \cdot 11^{v_1 + v_2 + 3v_3}\]

and we want all the powers of the prime factors to be even.

Note that the powers $v_i$ (which can only be $0$ or $1$) are our unknowns! We want to find these $v_i$’s that make the product a perfect square. So, we want to find $v_i$’s, each either $0$ or $1$, such that

(12.2)#\[\begin{split}\begin{align*} \phantom{2}v_1 + \phantom{5}v_2 \phantom{+ 3v_3} & \equiv 0 \pmod{2}, \\ 2v_1 + 5v_2 + \phantom{3}v_3 & \equiv 0 \pmod{2}, \\ 3v_1 + \phantom{5v_2 + } \phantom{3}v_3 & \equiv 0 \pmod{2}, \\ \phantom{2}v_1 + \phantom{5}v_2 + 3v_3 & \equiv 0 \pmod{2}. \end{align*}\end{split}\]

But, since $v_i$’s are either $0$ or $1$, we can think of this as a system in $\mathbb{F}_2 = \mathbb{Z}/2\mathbb{Z}$, and can then also reduce the coefficients modulo $2$, obtaining the system (in $\mathbb{F}_2$):

\[\begin{split}\begin{align*} v_1 + v_2 \phantom{+ v_3} & =0, \\ \phantom{v_1 +} v_2 + v_3 & =0, \\ v_1 + \phantom{v_2 +} v_3 & =0, \\ v_1 + v_2 + v_3 & =0. \end{align*}\end{split}\]

Systems in $\mathbb{F}_2$ are very easy to solve. We can use the same techniques we use with regular systems (with real numbers as coefficients), but remembering that in $\mathbb{F}_2$ we have that $x = -x$ (as $2$ divides $x - (-x) = 2x$).

So, the first equation gives us that $v_1 = -v_2 = v_2$. The second gives us that $v_2 = -v_3 = v_3$, so $v_1 = v_2 = v_3$. Now the fourth equation gives us that

\[0 = v_1 + v_2 + v_3 = v_1 + v_1 + v_1 = 3 v_1 = v_1,\]

and so $v_1 = v_2 = v_3 = 0$. So, the only way to get even exponents is to have all $v_i$’s equal to $0$, which means we use no $c_i$. But this does not help us! We need at least one $v_i$ to be $1$.

This means we need more $c_i$’s, so we go back to Relation Building to get more $c_i$’s. Suppose we get now:

\[\begin{split}\begin{align*} c_4 &= 3\phantom{^3} \cdot \phantom{5 \cdot 7^2 \cdot} 11^2, \\ c_5 &= 3^3 \cdot 5 \cdot 7^2. \end{align*}\end{split}\]

So, now we have

(12.3)#\[c_1^{v_1} c_2^{v_2} c_3^{v_3} c_4^{v_4} c_5^{v_5}= 3^{v_1 + v_2 + v_4 + 3v_5} \cdot 5^{2v_1 + 5 v_2 + v_3 + v_5} \cdot 7^{3v_1 + v_3 + 2v_5} \cdot 11^{v_1 + v_2 + 3v_3 + 2v_4}.\]

Again, we need the exponents to be even, so we set the system in $\mathbb{F}_2$ (already reducing the coefficients odulo $2$) having these exponents as $0$:

\[\begin{split}\begin{align*} v_1 + v_2 +\phantom{v_3 +} v_4 + v_5& =0, \\ \phantom{v_1 +} v_2 + v_3 + \phantom{v_4 +} v_5 & =0, \\ v_1 + \phantom{v_2 +} v_3 \phantom{+ v_4 + v_5}& =0, \\ v_1 + v_2 + v_3 \phantom{ + v_4 + v_5}& =0. \end{align*}\end{split}\]

The third equation gives us that $v_1 = v_3$. Then, the fourth gives us:

\[0 = v_1 + v_2 + v_3 = 2v_1 + v_2 = v_2,\]

so $v_2 = 0$. Now, substituting in the second, we get

\[0 = 0 + v_1 + v_5 \qquad \Longrightarrow \qquad v_5 = v_1,\]

and hence $v_1 = v_3 = v_5$ (and $v_2 = 0$). Finally, substituting it in the first, we have

\[0 = v_1 + 0 + v_4 + v_1 = 2v_1 + v_4 = v_4,\]

and hence $v_4 = 0$. So, all solutions of the system are of the form $(v_1, 0, v_1, 0, v_1)$, where $v_1$ is either $0$ or $1$. If we take $v_1 = 0$, then all exponents become zero again, which is pointless, so we choose $v_1 = 1$. This means that

\[c_1^1 \cdot c_2^0 \cdot c_3^1 \cdot c_4^0 \cdot c_5^1 = c_1 c_3 c_5\]

is a perfect square. In fact, from (12.3), we have

\[c_1 c_3 c_5 = 3^4 \cdot 5^4 \cdot 7^6 \cdot 11^4 = \left( 3^2 \cdot 5^2 \cdot 7^3 \cdot 11^2 \right)^2.\]

12.4.1.1. Solving the System with Sage#

How could we solve the system in Sage? Let’s first enter the $c_i$’s as a list in Sage. We will start with the system that did not work:

list_c = [3 * 5^2 * 7^3 * 11, 3 * 5^5 * 11, 5 * 7 * 11^3]

Now we need to create a matrix for our system. A matrix is simply a “table” containing the coefficients of the variable for the system. (More details about matrices can be found in the later section Matrices.) For example, here is our original system

\[\begin{split}\begin{align*} v_1 + v_2 \phantom{+ v_3} & =0, \\ \phantom{v_1 +} v_2 + v_3 & =0, \\ v_1 + \phantom{v_2 +} v_3 & =0, \\ v_1 + v_2 + v_3 & =0. \end{align*}\end{split}\]

The corresponding matrix is:

\[\begin{split}\begin{bmatrix} 1 & 1 & 0 \\ 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{bmatrix}.\end{split}\]

The first column of this matrix contains the coefficients (in $\mathbb{F}_2$) of $v_1$, the second contains the coefficients of $v_2$, and the third the coefficients of $v_3$.

We can create this matrix in Sage with:

M = matrix(GF(2), [[1, 1, 0], [0, 1, 1], [1, 0, 1], [1, 1, 1]])
M

[1 1 0]
[0 1 1]
[1 0 1]
[1 1 1]

So, a matrix is created using a list containing each row (as another list) as its elements.

Remember that GF(2), used as the first argument, is the same as FiniteField(2) and (basically the same as) Zmod(2), and gives the field $\mathbb{F}_2$. Using it as the first argument of matrix tells Sage that all the entries of the matrix are in the $\mathbb{F}_2$. In particular, if the entries were not already zeros and ones, they would be automatically reduced modulo $2$ when creating the matrix. So, we could have entered the original coefficients of the system before manually reducing them, i.e., (12.2):

M = matrix(GF(2), [[1, 1, 0], [2, 5, 1], [3, 0, 1], [1, 1, 3]])
M

[1 1 0]
[0 1 1]
[1 0 1]
[1 1 1]

To get the solutions, we use the .right_kernel() method:

M.right_kernel()

Vector space of degree 3 and dimension 0 over Finite Field of size 2
Basis matrix:
[]

The [] means that only all zeros is a solution, so this means that we need more $c_i$’s. In practice, we can see when this happens by saving the result in a variable, say K and then using the .dimension method. If we get $0$, then it means we only have the solution with all zeros.

K = M.right_kernel()
K.dimension()

So, we added $c_4$ and $c_5$:

list_c.extend([3 * 11^2, 3^3 * 5 * 7^2])

Then, we had the following matrix:

M = matrix(GF(2), [[1, 1, 0, 1, 3], [2, 5, 1, 0, 1], [3, 0, 1, 0, 2], [1, 1, 3, 2, 0]])
M

[1 1 0 1 1]
[0 1 1 0 1]
[1 0 1 0 0]
[1 1 1 0 0]

Then, we find the solutions again:

K = M.right_kernel()
K

Vector space of degree 5 and dimension 1 over Finite Field of size 2
Basis matrix:
[1 0 1 0 1]

The [1 0 1 0 1] is our solution non-zero solution, the same we found before! We can also see that we have found a solution from the dimension, which is not zero now:

K.dimension()

We can print all solutions by looping over the kernel:

for v in K:
    print(v)

(0, 0, 0, 0, 0)
(1, 0, 1, 0, 1)

To loop over non-zero solutions, we can use a slice:

for v in K[1:]:
    print(v)

(1, 0, 1, 0, 1)

But now we need to automate the creation of the matrix M from the $c_i$’s. We need the factorization of each one, but again, we would know a priori, from the Relation Building step, that all of them are $B$-smooth for some $B$. In this case, we know that they are all $11$-smooth. So, we test all primes less than $B$ and use valuation to get the exponent for each of these primes. These exponents will give us the coefficients of the system.

Let’s start by getting a list of all primes less than or equal to $11$:

list_p = list(primes(12))  # all primes less than or equal to 11

Now, let’s find the first row. It contains the powers of the first prime (i.e., $2$) for all $c_i$’s. The second row contains the powers of the second prime (i.e., $3$) for all $c_i$’s, and so on. So, each row is given by

[valuation(c, p) for c in list_c]

for a given p in list_p. Hence:

coef_matrix = matrix(ZZ, [[valuation(c, p) for c in list_c] for p in list_p])
coef_matrix

[0 0 0 0 0]
[1 1 0 1 3]
[2 5 1 0 1]
[3 0 1 0 2]
[1 1 3 2 0]

(The ZZ specifies that the entries are integers.)

Note that we are starting with coef_matrix that contains the exponents as integers, i.e., before being reduced modulo $2$, as it will be useful later. But now, we can obtain M (with entries in $\mathbb{F}_2$):

M = matrix(GF(2), coef_matrix)
M

[0 0 0 0 0]
[1 1 0 1 1]
[0 1 1 0 1]
[1 0 1 0 0]
[1 1 1 0 0]

Note that this matrix is not quite the same as before, since we introduced the prime $2$ that does not appear in our $c_i$’s. But that does not affect the result:

K = M.right_kernel()
K

Vector space of degree 5 and dimension 1 over Finite Field of size 2
Basis matrix:
[1 0 1 0 1]

And now, to produce our $b$’s. Since the solution is $(1, 0, 1, 0)$, we have that

\[b^2 = c_1^1 \cdot c_2^0 \cdot c_3^1 \cdot c_4^0 \cdot c_5^1 = c_1 \cdot c_3 \cdot c_5.\]

So, to get $b$ itself, we need to take a square root:

v = K[1]  # first non-zero solution
b = sqrt(prod(c^x for (c, x) in zip(list_c, v)))
b

Note that if we had list_a, with the corresponding $a_i$’s, we would obtain $a$ in a similar manner (without the square root):

a = prod(ai^x for (ai, x) in zip(list_a, v))

Note that when computing $b$ we need to take a square root. But we can speed up this calculation since we can easily now find the powers of the primes that appear in $b^2$ (as in (12.3)) and hence we can simply divide them by $2$.

The trick to find these exponents of the prime factors of $b^2$ is to use matrix multiplication. We just have to use each solution found and multiply coef_matrix (and not M) by it:

coef_matrix * vector(ZZ, list(v))  # note the *list* there!

(0, 4, 4, 6, 4)

(Note that for this multiplication to work, we need to convert the solution (out list of zeros and ones) to a vector with integer entries.)

So, this is saying, in terms of matrix multiplications, that

\[\begin{split}\underbrace{\begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 3 \\ 2 & 5 & 1 & 0 & 1 \\ 3 & 0 & 1 & 0 & 2 \\ 1 & 1 & 3 & 2 & 0 \end{bmatrix}}_{\text{coefficient matrix}} \cdot \underbrace{\begin{bmatrix} 1 \\ 0 \\ 1 \\ 0 \\ 1 \end{bmatrix}}_{\text{solution}} = \underbrace{\begin{bmatrix} 0 \\ 4 \\ 4 \\ 6 \\ 4 \end{bmatrix}}_{\text{exponents of the primes}}\end{split}\]

This means that

\[b^2 = 2^0 \cdot 3^4 \cdot 5^4 \cdot 7^6 \cdot 11^4,\]

and hence

\[b = 2^0 \cdot 3^2 \cdot 5^2 \cdot 7^3 \cdot 11^2.\]

So, we can divide these exponents by $2$ to obtain the correct exponents for the prime factors of $b$.

exponents = 1/2 * coef_matrix * vector(ZZ, list(v))
b = prod(p^x for (p, x) in zip(list_p, exponents))
b

We get the same result (but faster)!

Note

This method only works for $b$ and not for $a$, which we would still compute as above:

a = prod(ai^x for (ai, x) in zip(list_a, v))

Important

If we are trying to factor $N$, as we would in practice, we should reduce each $a$ and $b$ modulo $N$, as in

a = ZZ(prod(Mod(ai, N)^x for (ai, x) in zip(list_a, v)))
b = ZZ(prod(Mod(p, N)^x for (p, x) in zip(list_p, exponents)))

to make speed up the computation of the GCD. Note that we compute the products in $\mathbb{Z}/N\mathbb{Z}$ and then convert them back to $\mathbb{Z}$ to make this reduction modulo $N$ more efficient.

In this case, we only have one $b$, as we’ve seen above, since we have only one non-zero solution, but we might have more if there is more than a single solution.

What we do is to loop over the solutions in K[1:]. At each iteration we find $a$ and $b$ (as shown above), and compute $\gcd(N, a - b)$. It is neither $1$ not $N$, we found our divisor. If not, we proceed the next iteration (i.e., the next solution). If we run through all solutions, never finding a divisor, the method failed.

Homework

You will implement this Elimination process in your homework.

12.5. Smooth Numbers#

As previously mentioned, smooth numbers appear in all modern factorization methods (as well as other computational problems). In this section we will give some mathematical properties of $B$-smooth numbers. The proofs are complex and beyond the scope of these notes, so we will simply state the properties. But let’s start with some examples.

Example 12.1 ($5$-Smooth Numbers)

We can break down the integers into $5$-smooth numbers and not $5$-smooth:

\[\begin{split}\begin{align*} \text{$5$-smooth:} &\quad 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, \ldots \\ \text{not $5$-smooth:} &\quad 7, 11, 13, 14, 17, 19, 21, 22, 23, 26, 28, 29, \ldots \end{align*}\end{split}\]

We now introduce a function to count the number of $B$-smooth numbers up to a given bound:

Definition 12.3 (Number of Smooth Numbers)

We define $\psi(x, B)$ to be the number of smooth numbers less than or equal to $x$.

So, we can see from the example above that $\psi(30, 5) = 17$.

We also need the following function:

Definition 12.4 (Function $L$)

We define function $L(x)$, for $x >0$, as

\[L(x) = e^{\sqrt{\log(x) \cdot \log(\log(x))}}.\]

(Remember that $\log(x)$ is the natural log, i.e., the log base $e \approx 2.7183$).

We then have the following technical result:

Theorem 12.1

For any real number $c$ with $0 < c < 1$ and for “large” values of $x$, we have

\[\psi(x, L(x)^c) \approx x \cdot L(c)^{- \frac{1}{2c} \cdot (1 + s_x)},\]

where $s_x$ approaches $0$ as $x$ gets large.

When using the difference of squares factorization method, in order to find solution we need that the system must have more equations than variables. The number of equations is given by the number of primes less than or equal to $B$. When we require the $c_i$’s to be $B$-smooth this number is $\pi(B)$ (the number of primes less than or equal to $B$). The number of variables is given by the number of $a_i$’s (or $c_i$’s). So, we need the number of $B$-smooth $c_i$’s to be greater than $\pi(B)$.

The theorem above gives us the an idea of how difficult finding that many $B$-smooth numbers is:

Theorem 12.2

Let $N$ be a large integer and set $B = L(N)^{\frac{1}{\sqrt{2}}}$. We expect to check about $B^2 = L(N)^{\sqrt{2}}$ random numbers modulo $N$ to find $\pi(B)$ numbers that are $B$-smooth. Thus, we expect to check $B^2 = L(N)^{\sqrt{2}}$ random numbers of the form $a^2$ reduced modulo $N$ to be able to factor $N$.

Here is a table to give you an idea of the size of $L(N)$ and $L(N)^{\sqrt{2}}$:

Table 12.1 Sizes of $L(N)$ and $L(N)^{\sqrt{2}}$ compared to $N$.#
$N$	$L(N)$	$L(N)^{\sqrt{2}}$
$2^{100}$	$2^{24.73}$	$2^{34.97}$
$2^{250}$	$2^{43.12}$	$2^{60.98}$
$2^{500}$	$2^{64.95}$	$2^{91.85}$
$2^{1000}$	$2^{97.14}$	$2^{137.38}$
$2^{5000}$	$2^{242.48}$	$2^{342.91}$

Or, in terms of the number of digits

Table 12.2 Numer of digits of $L(N)$ and $L(N)^{\sqrt{2}}$ compared to $N$.#
$\log_{10}(N)$	$\log_{10}(L(N))$	$\log_{10}(L(N)^{\sqrt{2}})$
$31$	$8$	$11$
$76$	$13$	$19$
$151$	$20$	$28$
$302$	$30$	$42$
$1506$	$73$	$104$

One can improve this process by taking the $a_i$’s, which we square and reduce modulo $N$, only slightly larger than $\sqrt{N}$ (instead of purely randomly), which reduces the number of $a_i$’s to try to about $L(N)$, which is a good improvement.

We will come back to this process (the Relation Building) in a later section.

12.5.1. Checking for $B$-Smoothness in Sage#

If you want to check that some number $a$ is $B$-smooth in Sage, one should not use the factor function. If the number is $B$-smooth, for a relatively small $B$, that could work.

a = 2^6 * 7^3 * 11^4 * 13^5 * 23^2
a

63127307789850304

B = 23
a = 63127307789850304

# test
prime_divisors = [x[0] for x in factor(a)]

if max(prime_divisors) <= B:
    print(f"{a} is {B}-smooth.")
else:
    print(f"{a} is not {B}-smooth.")

63127307789850304 is 23-smooth.

Note that Sage actually has a prime_factors function that could make this even easier:

B = 23
a = 63127307789850304

# test
if max(prime_factors(a)) <= B:
    print(f"{a} is {B}-smooth.")
else:
    print(f"{a} is not {B}-smooth.")

63127307789850304 is 23-smooth.

But, if not, it will take a lot longer than necessary!

%%time

a = 505315045110283216972911558435118337
B = 23

# test
if max(prime_factors(a)) <= B:
    print(f"{a} is {B}-smooth.")
else:
    print(f"{a} is not {B}-smooth.")

505315045110283216972911558435118337 is not 23-smooth.
CPU times: user 8.32 ms, sys: 30 µs, total: 8.35 ms
Wall time: 8.33 ms

The idea is to simply see how many times each prime less than or equal to $B$ divides the number, and see if only those primes raised to their corresponding powers and multiplied give the number.

Remember that we can use valuation(a, p) in Sage to see how many times the prime p divides a. Hence, to check for $B$-smoothness, we can do

a == prod(p^(valuation(a, p)) for p in primes(B + 1))

For example:

%%time

a = 505315045110283216972911558435118337
B = 23

# test
if a == prod(p^(valuation(a, p)) for p in primes(B + 1)):
    print(f"{a} is {B}-smooth.")
else:
    print(f"{a} is not {B}-smooth.")

505315045110283216972911558435118337 is not 23-smooth.
CPU times: user 57 µs, sys: 0 ns, total: 57 µs
Wall time: 58.4 µs

It’s a lot faster, as it only tests divisibility by small numbers!

Here is some information about the computer used for the computations above:

!inxi --system --cpu --memory

System:
  Host: dell7010 Kernel: 6.17.8-1-siduction-amd64 arch: x86_64 bits: 64
  Desktop: KDE Plasma v: 6.5.3 Distro: siduction 22.1.2 Masters_of_War -
    kde - (202303220911)
Memory:
  System RAM: total: 128 GiB available: 125.51 GiB used: 29.51 GiB (23.5%)
  Array-1: capacity: 128 GiB slots: 4 modules: 4 EC: None
  Device-1: DIMM1 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
    actual: 3600 MT/s
  Device-2: DIMM2 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
    actual: 3600 MT/s
  Device-3: DIMM3 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
    actual: 3600 MT/s
  Device-4: DIMM4 type: DDR5 size: 32 GiB speed: spec: 4800 MT/s
    actual: 3600 MT/s

CPU:
  Info: 24-core model: 13th Gen Intel Core i9-13900 bits: 64 type: MCP cache:
    L2: 32 MiB
  Speed (MHz): avg: 5000 min/max: 800/5300:5600:4200 cores: 1: 5000 2: 5000
    3: 5000 4: 5000 5: 5000 6: 5000 7: 5000 8: 5000 9: 5000 10: 5000 11: 5000
    12: 5000 13: 5000 14: 5000 15: 5000 16: 5000 17: 5000 18: 5000 19: 5000
    20: 5000 21: 5000 22: 5000 23: 5000 24: 5000

\(N\)	\(L(N)\)	\(L(N)^{\sqrt{2}}\)
\(2^{100}\)	\(2^{24.73}\)	\(2^{34.97}\)
\(2^{250}\)	\(2^{43.12}\)	\(2^{60.98}\)
\(2^{500}\)	\(2^{64.95}\)	\(2^{91.85}\)
\(2^{1000}\)	\(2^{97.14}\)	\(2^{137.38}\)
\(2^{5000}\)	\(2^{242.48}\)	\(2^{342.91}\)

\(\log_{10}(N)\)	\(\log_{10}(L(N))\)	\(\log_{10}(L(N)^{\sqrt{2}})\)
\(31\)	\(8\)	\(11\)
\(76\)	\(13\)	\(19\)
\(151\)	\(20\)	\(28\)
\(302\)	\(30\)	\(42\)
\(1506\)	\(73\)	\(104\)

Factorization

Contents

12. Factorization#

12.1. Naive Factorization#

12.2. Pollard’s \(p-1\) Factorization Algorithm#

12.2.1. Number of Operations#

12.2.2. Example#

12.3. Choosing \(p\) and \(q\) for the RSA Cryptosystem#

12.4. Factorization via Difference of Squares#

12.4.1. Example#

12.4.1.1. Solving the System with Sage#

12.5. Smooth Numbers#

12.5.1. Checking for \(B\)-Smoothness in Sage#