A Straightforward Implementation of Dixon's Factoring Algorithm
When studying RSA cryptanalysis you may be asked to implement, as a learning exercise, Dixon’s random squares algorithm for factoring composite numbers of two primes. Unfortunately, many implementations demonstrate more powerful versions of the algorithm. They take shortcuts, make optimizations, or implement a different algorithm entirely (like the number field sieve). This is not helpful when you’re trying to learn how Dixon’s works in a straightforward way. In this post I’ll lay out a simple implementation of Dixon’s algorithm in C.
Note: If you’re here for your 695.641 Cryptology homework, several components of the working code are left as an exercise to the reader. The point is to see how the algorithm can be implemented in actual code. A lot of the bookkeeping is omitted, copying and pasting the snippets verbatim won’t work.
What is Dixon’s Factoring Algorithm?
I’m not going to entirely restate Dixon’s factoring algorithm here. There are plenty of resources available online for that, including the original paper.
In short, Dixon’s algorithm sets out to find some $x$ and $y$ such that $x \not\equiv y \pmod{n}$ and $x^2 \equiv y^2 \pmod{n}$. When these conditions are satisfied, we know we have a good chance of finding a non-trivial factor of $n$ by taking the greatest common divisor of $n$ and $(x+y)$ and/or $(x-y)$, i.e., $\text{gcd}(x+y, n)$ and $\text{gcd}(x-y, n)$. Finding non-trivial factors of a product $n$ of two large primes $p$ and $q$ is the first step in breaking RSA.
The algorithm also relies on a factor base, $B$, which is the set of all prime numbers less than some number $b$. The first step of the algorithm involves finding some $z$ values such that their $z^2 \pmod{n}$ are “smooth” over $B$. In plain English, we should be able to reconstruct $z^2 \pmod{n}$ with the primes in $B$. For example, if $z = 100$, $n=200000$, and $b = 10$ we have the following:
Let $a$ be a set of powers of each prime in $B$. Then we can rewrite $z^2 \pmod{n}$ as $a = \lbrace 0, 4, 0, 4, 0 \rbrace$. We can reconstruct $z^2 \pmod{n}$ with $1^0 \times 2^4 \times 3^0 \times 5^4 \times 7^0$. Because we can rewrite $z^2$ with all the primes in $B$, we can say it is $B$-smooth. When we find a bunch of $z$ values that satisfy this condition we can turn them into our $x$ and $y$ values and, hopefully, factor $n$.
The Algorithm
That should be enough background for the pseudocode makes sense. If not, read one of linked resources. We take the following pseudocode from Wikipedia:
input: positive integer $N$
output: non-trivial factor of $N$
Choose bound $B$
Let $P:=\lbrace p_1, p_2, \dots, p_k \rbrace$ be all primes $\leq B$
repeat
for $i = 1$ to $k+1$ do
Choose $0 < z_i < N$ such that $z_i^2\mod N$ is $B$-smooth
Let $a_i := \lbrace a_{i1}, a_{i2}, \dots, a_{ik} \rbrace$ such that $z_i^2\mod N = \prod_{p_j \in P}p_j^{a_{ij}}$
end for
Find non-empty $T \subseteq \lbrace 1, 2, \dots, k+1 \rbrace$ such that $\sum_{i \in T}a_i \equiv \vec{0} \pmod{2}$
Let $x := \left(\prod_{i\in T}z_i\right) \mod N$
$y : \left(\prod_{p_j\in P}p_j^{\left(\sum_{i\in T}a_{ij}\right)/2}\right)\mod N$
while $x \equiv \pm y \pmod{N}$
return gcd$(x+y, N)$
Everything up to the end of the for loop is what I explained in the previous section. The while condition keeps the loop going until we find our $x \not\equiv y \pmod{n}$ such that $x^2 \equiv y^2 \pmod{n}$. Then, we return the greatest common divisor of $x+y$ and $n$. The lines between the end for and while are how we use $z$ to find candidate $x$ and $y$.
The Code
Here’s my problem: all the implementations I can find online are not great learning tools. It appears (as Wikipedia suggest) that most of these implementations are actually the block Lanczos algorithm. Or they’re just doing something else entirely. Seriously, check it out:
- GeeksforGeeks’s implementation does not match their described steps. It works, but I’m not sure what this algorithm actually is.
- This paper works, but is using Kraitchik’s variation which is very similar, but slightly different.
- This blog gets close, but is a more “industrial scale” implementation. It’s more efficient, but try putting it side by side with the pseudocode. It doesn’t line up.
Let’s implement the algorithm according to the above pseudocode instead. It won’t be the most efficient (it’s actually quite slow), but it will be instructive. We’ll reuse the numbers from the Wikipedia example for consistency, $N = 84923$ and bound $B = 7$.
Note: My code will use the following shortcut for long long int:
typedef long long int l_int;
Let’s start at the top.
input: positive integer $N$
output: non-trivial factor of $N$Choose bound $B$
Let $P:=\lbrace p_1, p_2, \dots, p_k \rbrace$ be all primes $\leq B$
l_int n = 84923
l_int p[] = {2, 3, 5, 7}
size_t size = sizeof(p) / sizeof(p[0]);
Here we set $n = 84923$ and build our set $P$ from the bound $B = 7$. Moving on to the next lines…
repeat
The whole algorithm is in a while loop. Let’s hold off on that for now.
for $i = 1$ to $k+1$ do
This loop will move through $z_i$ values searching for candidates. Choosing $z_i$ values is outside the scope of this post, but to follow the wikipedia example let’s say we start with $z_1 = 500$ and $z_{100} = 599$. Then we’d have:
l_int z = 500;
for (int i = 0; i < 100; i++, z++) {
/* test if z is smooth */
}
Choose $0 < z_i < N$ such that $z_i^2\mod N$ is $B$-smooth
Let $a_i := \lbrace a_{i1}, a_{i2}, \dots, a_{ik} \rbrace$ such that $z_i^2\mod N = \prod_{p_j \in P}p_j^{a_{ij}}$
end for
So how do we test our $z_i^2 \mod n$ values for $B$-smoothness? Trial division! We also want to save the powers of each prime factor we find that can reconstruct $z_i^2$ as shown in the second line above.
/* first find z^2 mod n */
l_int z_square = (z * z) % n;
/* this vector will hold the powers of each prime */
int a[size];
memset(a, 0, size * sizeof(a[0]));
/* Now we do trial division to whittle down z_square with our primes in p. */
for (int i = 0; i < size; i++) {
while ((z_square % p[i]) == 0) {
a[i]++;
z_square /= p[i];
}
}
if (z_square == 1)
/* we got all the way down to 1 with primes in p, z^2 is b-smooth */
return 1;
else
/* there is still some z^2 left, z^2 is not b-smooth */
return 0;
Note: If you’re here for the homework, negative values in p[i] will mess
this code up. C does not do modulo the way you’d expect with negative numbers,
think of ways you can convert negative numbers to positive ones $\pmod{n}$.
You may also be able to form a solution with both the original $z_i^2 \pmod{n}$
and $(-1)z_i^2\pmod{n}$ being used.
You’ll want to put this code into a function you can call for each candidate $z_i$ value, and you’ll want a way to keep track of the $z_i$ values and subsequent $a$ vectors where $z_i^2$ is $B$-smooth. I’ll leave keeping track of $B$-smooth numbers and $a$ vectors up to the reader. If you’re on Linux, I recommend putting them together in a struct and linking them in a list with BSD’s queue. If you’re not, you can use a language with dynamically resizable sets/lists like Python.
Here is the struct I used:
LIST_HEAD(listhead, entry) head;
struct e_listhead *headp;
typedef struct entry {
l_int z;
l_int a[size];
LIST_ENTRY(entry) entries;
} candidate_entry;
LIST_INIT(&head);
Then when I find a $B$-smooth candidate, I simply:
candidate_entry *candidate = malloc(sizeof(candidate_entry));
/* insert here: assign z and a to candidate */
LIST_INSERT_HEAD(&head, candidate, entries);
Moving on, now we should have a list of $z_i$ values and corresponding $a$ vectors. The next step is:
Find non-empty $T \subseteq \lbrace 1, 2, \dots, k+1 \rbrace$ such that $\sum_{i \in T}a_i \equiv \vec{0} \pmod{2}$
As you might notice, this line leave a lot up to the implementer. Just find some subsets! East, right? In English, what we’re trying to do here is find all possible combinations of our $a_i$ vectors where each entry sums to an even number. For example:
would qualify, because the resulting vector $\lbrace 0, 2, 0, 2 \rbrace$ has values that are all even (i.e., equal $0 \pmod{2}$).
In other words, we need to sum every possible combination of all the $a$
vectors we just found. We can do this by finding the
power set
of the set of all our $a$ vectors, minus the empty set. The easiest way
to do this is by generating a set of bit-masks where a $1$ in position $i$
indicates we will include $a_i$ in the sum, and a $0$ means we will ignore
it. Note that this loop will iterate $2^n$ times where $n$ is the number of
$B$-smooth numbers we found. smooth_counter is the number of $B$-smooth
numbers we found, and $2^n$ is given by (1 << smooth_counter). We start
the loop at 1 to skip the empty set.
for (int counter = 1; counter < (1 << smooth_counter); counter++) {
/* make a mask for a candidate */
int mask[smooth_counter];
memset(mask, 0, smooth_counter * sizeof(int));
for (int i = 0; i < smooth_counter; i++) {
if (counter & (1 << i))
mask[i] = 1;
}
/* compute vector from mask */
int result_vector[size];
memset(result_vector, 0, size * sizeof(int));
for (int i = 0; i < smooth_counter; i++) {
if (mask[i]) {
for (int j = 0; j < size; j++)
result_vector[j] += a[j];
}
}
/* test if vector sums to zero, if we have an odd value skip the rest */
int flag = 0;
for (int i = 0; i < size; i++) {
if ((result_vector[i] % 2) != 0)
flag = 1;
}
if (!flag) {
/* save the mask and vector */
}
}
Remember, size is the number of primes we have in $B$.
Again, the exact method for saving the resulting non-empty $T$ in the last
if block is up to the reader. And again, I recommend using a struct and
BSD’s queue as shown previously. We’re almost there, just two more significant
lines to deal with.
Let $x := \left(\prod_{i\in T}z_i\right) \mod N$
$y : \left(\prod_{p_j\in P}p_j^{\left(\sum_{i\in T}a_{ij}\right)/2}\right)\mod N$
We will assume we saved mask[] and result_vector[] from the previous section
and that we have our $z$ values in some array z[].
/* calculate the x product */
l_int x_product = 1;
for (int i = 0; i < smooth_counter; i++) {
if (mask[i])
x_product = (x_product * z[i]) % n;
}
/* calculate the y product */
l_int y_product = 1;
for (int i = 0; i < size; i++) {
/* y = p^a[ij]/2 */
for (int j = 0; j < result_vector[i] / 2; j++)
y_product = (y_product * p[i]) % n;
}
You might notice some ways these calculations can be integrated into the
previous code snippet that would make the code run faster. Remember, these
x_product and y_product values need to be saved too.
while $x \equiv \pm y \pmod{N}$
Let’s address this part now. The pseudo code as written has the whole thing loop for each candidate $z$ value. The way I’ve written this code assumes that we’re operating on a single set of $z$ candidates. If we can’t find a result, we simply fail instead of trying again with a new set of $z$s. Still, we need to test our candidate $x$ and $y$ pairs to find one that satisfies $x \not\equiv y \pmod{n}$ and $x^2 \equiv y^2 \pmod{n}$. The following code assumes we are iterating over all the $x$ and $y$ values previously saved. That is, we keep looping until we find the answer.
return gcd$(x+y, N)$
We will also include our return value, it follows naturally from the test.
/* shorter variable names */
l_int x = x_product;
l_int y = y_product;
/* check if x is not congruent y mod n when x^2 is congruent y^2 mod n */
if ((x % n) != ((y + n) % n) && ((x * x) % n) == ((y * y) % n)) {
l_int gcd = GCD_ALGORITHM(x + y, n);
if (gcd != 1 && gcd != n) /* exclude trivial answers */
return gcd;
}
GCD_ALGORITHM is some function that calculate the greatest common divisor
of $x+y$ and $n$. The Extended Euclidean algorithm
is a good option. If we loop through all our x_products and y_products
without finding a non trivial factor we can return -1 to indicate an
error.
return -1;
We might want to return both factors, which we can get with
l_int gcd_one = GCD_ALGORITHM(x + y, n);
l_int gcd_two = GCD_ALGORITHM(x - y, n);
And that’s it, you can now factor numbers with Dixon’s random squares algorithm. Here’s what a trial run might look like with $N = 84923$ and bound $B = 7$. I added some prints so you can see what’s going on and test some of the numbers by hand.
user@linux:~$ cc dixon.c
user@linux:~$ ./a.out
Found b-smooth: 505^2
Found b-smooth: 513^2
Found b-smooth: 537^2
Found b-smooth: 655^2
Found b-smooth: 668^2
x=505, y=16, gcd=521
The factors of 84923 are 521 and 163
Note that in this case we only needed to let the code run until it found both $513$ and $537$, as we know these to make up the answer. For a composite with an unknown solution you can tweak the number of $z$ values searched until you find a solution.