Special-purpose factorization algorithms

In this post, we will take a brief look at a nice collection of special-purpose factorization algorithms. Most of them old and well-known. Some of them new.

By a "special-purpose factorization algorithm", we mean an algorithm that can find one or more factors of a large integer

$n$ that satisfies a specific special form, without strongly depending on the size of

$n$ .

When the special conditions are met, a factor is usually found very fast. On the other hand, when the number does not have the required special form, the algorithm may fail to produce a factor in a reasonable amount of time, or fail entirely.

Special-purpose factorization algorithms are commonly used in general-purpose factorization routines as a pre-factorization test, trying to reduce the size of

$n$ by finding easy factors, before sending the number to the heavy general-purpose modern machineries, such as the Quadratic Sieve (QS) or the Number Field Sieve (NFS) algorithms, which strongly depend on the size of

$n$ .

The most commonly used special-purpose factorization algorithms are trial division, Pollard's rho method, Pollard's p-1 method and Lenstra's elliptic-curve factorization method (ECM).

# Perfect power check

$n = r^k$ , for some positive integer

$r$ and some integer

$k >= 2$ , then, in order to find the factorization of

$n$ , we only need to factorize

$r$ and duplicate each prime factor accordingly.

This can be efficiently checked by iterating over each

$k$ , starting from

$k = floor(log_2(n))$ , down to

$k=2$ , and checking if

$n^(1/k)$ is an integer.

For checking if

$n^(1/k)$ is an integer, let

$r = floor(n^(1/k))$ , which can be computed efficiently using Newton's method. If

$n = r^k$ , then

$n$ is a perfect power.

This is not considered a factorization algorithm, but it's an import thing to check when factoring an unknown integer.

# Trial division

The simplest way to try to find a small factor of a given integer

$n$ , is to set a bound

$B$ , at most

$B = floor(sqrt(n))$ , and iterate over the prime numbers

$p$ in the range

$[2, B]$ and check if any such

$p$ divides

$n$ .

This algorithm is of practical use only for numbers that have at least a small factor (say, less than

$10^8$ ).

If multiple numbers need to be checked for small factors, a simple optimization would be computing the primorial of

$B$ only once, which is the product of all the primes in the range

$[2, B]$ , say

$P = B#$ , and using the pre-computed value of

$P$ to compute

$d = gcd(P, n)$ for multiple values of

$n$ , where

$d$ is a divisor of

$n$ (not necessarily prime), or

$1$ if

$n$ has no prime factors in the range

$[2,B]$ .

# Difference of squares

By representing a given integer

$n$ as a difference of two squares:

$n = a^2 - b^2$

it can then be factorized as:

$n = (a + b) (a - b)$

In Fermat's factorization method, we start by setting

$a = floor(sqrt(n))$ and checking if

$(a+k)^2 - n$ is a perfect square, for some

$k >= 1$ .

This method is effective when

$n$ has a factor near

$sqrt(n)$ .

# Congruence of squares

An improvement over Fermat's factorization method, noticed by Maurice Kraitchik, is to try to find a congruence of squares modulo

$n$ , of the form

$x^2 ≡ y^2 mod n$ , with the condition that

$x ≢ ±y mod n$ , which factorizes

$n$ as:

$n = gcd(x+y, n) * gcd(x-y, n)$

Same as in Fermat's factorization method, we can try to find a congruence of squares by setting

$a = floor(sqrt(n))$ and checking if

$(a + k)^2 mod n$ is a perfect square, for some

$k >= 1$ .

However, a slightly better approach would be to use the continued fraction convergents

$A_k/B_k$ of

$sqrt(n)$ and testing if

$A_k^2 mod n$ is a perfect square, for some

$k >= 1$ .

When the continued fraction convergents

$A_k/B_k$ of

$sqrt(n)$ are used, we call it the "Pell factorization method", because it uses the same technique as in solving the Pell equation

$x^2 - n y^2 = 1$ , which has solutions of the form

$A_k^2 - n B_k^2 = 1$ , for some

$k >= 1$ .

This method is effective when

$n$ has a factor near

$sqrt(n)$ .

The idea of creating a congruence of squares, is used in many modern general-purpose factorization algorithms, including the CFRAC method, which we described in an earlier post.

# Difference of powers

A generalization of Fermat's factorization method, would be to check if

$n$ has the form

$n = a^k ± b^k$ , for some positive integers

$a$ and

$b$ , with odd

$k >= 3$ . Then either

$gcd(a-b, n)$ or

$gcd(a+b, n)$ is a non-trivial factor of

$n$ .

This method is based on the observation that, for odd

$k >= 3$ :

$(a^x - b^y)$ divides

$(a^(k x) - b^(k y))$

$(a^x + b^y)$ divides

$(a^(k x) + b^(k y))$

Without loss of generality, let's assume that:

$n = r_1^(e_1) ± r_2^(e_2)$

We can try to find such a representation by fixing

$B = floor(log_2(n))$ and iterating over each

$k$ in the range

$[2, B]$ .

For each

$k$ , we check if any of the following values is a perfect power (let's call it

$r_2^(e_2)$ ):

${n - k^x, k^(x+1) - n, n - y^k, (1+y)^k - n}$

where

$x = floor(log_(k)(n))$ and

$y = floor(n^(1/k))$ .

For each such perfect power found, we can represent

$n$ as

$n = r_1^(e_1) ± r_2^(e_2)$ .

Next, we iterate over each divisor

$d_1$ of

$e_1$ and over each divisor

$d_2$ of

$e_2$ , checking if

$gcd(r_1^(d_1) ± r_2^(d_2), n)$ is a non-trivial factor of

$n$ .

This method is effective for numbers of the form

$a^x ± b^y$ , where

$x$ and

$y$ have many common divisors.

# Congruence of powers

In this method, we generalize Kraitchik's idea to find congruence of powers of the form

$a^k ≡ b^k mod n$ , for some positive integers

$a$ and

$b$ , with odd

$k >= 3$ . Then either

$gcd(a-b, n)$ or

$gcd(a+b, n)$ is a non-trivial factor of

$n$ .

Same as in the previous method, we try to find such a representation by fixing

$B = floor(log_2(n))$ and iterating over each

$k$ in the range

$[2, B]$ .

For each

$k$ , we check if either one of the following congruences is true:

$x^k ≡ b^k mod n$

$(x+1)^k ≡ b^k mod n$

where

$x = floor(n^(1/k))$ .

For each such congruence found, we have

$a^k ≡ b^k mod n$ , where either

$gcd(a + b, n)$ or

$gcd(a - b, n)$ may be a non-trivial factor of

$n$ .

# Cyclotomic polynomials

For some numbers, we can also use the cyclotomic polynomials

$Φ_n(x)$ for finding factors of special form.

Given a composite integer

$n$ and a fixed value of

$b >= 2$ , we iterate over each

$k$ in the range

$2 <= k <= floor(1+\log_b(n))$ , checking if

$gcd(Φ_k(b), n)$ is a non-trivial factor of

$n$ .

In practice, we can check each value of

$b$ in the range

$2 <= b <= floor(log_2(n))$ .

# HOLF method

In Hart's One-Line Factoring Algorithm, we try to find a congruence of squares of the form:

$a^2 ≡ b^2 mod n$

where

$a = 1 + floor(sqrt(n k))$ , for some positive integer

$k$ .

Then

$gcd(a - b, n)$ may be a non-trivial factor of

$n$ .

This method is effective for numbers of the form

$n = x * (x k ± z)$ , where

$x$ can be arbitrary large, but

$k$ and

$z$ are small.

# Carmichael method

In some sense, this is a generalization of the HOLF method and can factorize numbers of the form:

$n = x * \prod_{k=1}^r ((x±1) a_k ± 1)$

for a given value

$r$ and a small range

$[L,R]$ that includes all the values

$a_k$ , where the set

${a_k}$ is a combination of

$r$ elements in the given range

$[L,R]$ .

This method can be efficiently implemented using the binary search algorithm to find the value

$x$ in the range:

$2 <= x <= floor(n^(1/(r+1)))$

Numbers of this form include the Carmichael numbers and the Lucas pseudoprimes, and can be factorized relatively fast by this method when

$r$ is small and the range

$[L,R]$ is not very wide.

# Miller-Rabin $n-1$ method

The Miller-Rabin primality test can sometimes produce a factor of

$n$ . By writing

$n-1$ as:

$n-1 = 2^s d$

with

$s >= 1$ and

$d$ odd, we try to find

$x ≡ a^d mod n$ , for some base

$a >= 2$ , such that

$x ≢ ±1 mod n$ , and

$x^(2^k d) ≡ 1 mod n$ , for some

$k$ in the range

$1 <= k < s$ .

If such congruence is found, then

$y = x^(2^(k-1) d) mod n$ factorizes

$n$ as:

$n = gcd(y+1, n) * gcd(y-1, n)$

Alternatively, a simpler version would be to compute

$x = a^(2^k d) mod n$ and check if

$x ± 1$ has a common factor with

$n$ , for some

$k$ in the range

$0 <= k < s$ and some random value of

$a$ .

This method works best in factoring Carmichael numbers and some Fermat pseudoprimes.

# Lucas-Miller $n+1$ method

In the same spirit as the Miller-Rabin factorization method, we can take advantage of the Lucas sequence

$U_n(P,Q)$ , writing

$n+1$ as:

$n+1 = 2^s d$

with

$s >= 1$ and

$d$ odd, and checking if

$U_(d 2^k)(P,Q)$ (computed modulo

$n$ ) has a common factor with

$n$ , for some

$k$ in the range

$0 <= k < s$ and some random values of

$P$ and

$Q$ .

This method works best in factoring Lucas-Carmichael numbers and some Lucas pseudoprimes.

Another version of this method, is the

$n-1$ variation, in which we write

$n-1 = 2^s d$ , with

$d$ odd, and we check

$V_(d 2^k)(P,Q)$ and

$V_(d 2^k)(P,Q) - P$ , which is able to factorize Carmichael numbers and Fermat pseudoprimes.

# FLT factorization method

Inspired by Fermat's Little Theorem (FLT), we can try to find a non-trivial factor of

$n$ by checking

$gcd((2^n mod n) - (2^k mod n), n)$ for some small odd

$k >= 1$ .

With

$n = p*q$ , let:

$a = 2^p mod n$

$b = 2^q mod n$

then:

$2^n ≡ a+b-2 mod n$

where:

$gcd((2^n mod n) - (2^(q-p + 1) mod n), n)$

is a non-trivial factor of

$n$ .

The method can be efficiently implemented iteratively, by choosing a base

$b >= 2$ and setting:

$z = b^n mod n$

$w_0 = b$

$z != b$ , then for

$k >= 1$ , we iteratively do:

$w_k = w_{k-1} * b^2 mod n$

at each step checking if

$gcd(z - w_k, n)$ is a non-trivial factor of

$n$ .

This method is effective, in particular, when

$n$ has all of its prime factors close to each other.

More generally, we can consider any C-finite sequence

$f(n)$ , trying to find a factor of

$n$ by checking:

$gcd((f(n) mod n) - (f(k) mod n), n)$

$gcd((f(n) mod n) - (f(1+floor(n/k)) mod n), n)$

for several small

$k >= 1$ .

# Pollard rho method

Pollard's rho method tries to find a congruence of the form

$f(x_i) ≡ f(x_j) mod p$ , for some polynomial

$f(x)$ , where

$p$ is a prime factor of

$n$ .

For example, let:

$f(x) = x^2 + 1$ . If

$x_i ≡ x_j mod p$ , then

$f(x_i) ≡ f(x_j) mod p$ .

Since we do not know the value of

$p$ , we have to work modulo

$n$ and assume that there exists some

$p$ that divides

$n$ which satisfies the wanted congruence.

If

$f(x_i) ≡ f(x_j) mod p$ , this means that

$gcd(f(x_i) - f(x_j), n)$ is a divisor of

$n$ .

This method can be efficiently implemented by using a cycle detection algorithm, such as Floyd's cycle detection algorithm, starting by initializing the following two sequences modulo

$n$ , for some value of

$a$ :

$x_0 = f(a) mod n$

$y_0 = f(f(a)) mod n$

then, for

$i >= 0$ , we iteratively do:

$x_(i+1) = f(x_i) mod n$

$y_(i+1) = f(f(y_i)) mod n$

at each step checking if

$gcd(x_i - y_i, n)$ is a non-trivial factor of

$n$ .

By the birthday paradox, the Pollard rho method is expected to find a prime factor

$p$ of

$n$ in

$O(sqrt(p))$ steps.

The polynomial

$f(x)$ can be optimized based on what is known about the factors of

$n$ . In general, the most used polynomial is

$f(x) = x^2 + c$ , for some constant

$c$ .

A variation of this method, that we call the "rho-exp" method, uses the following interesting polynomial:

$f(x) = x^k + 2 k - 1$

where

$k = lcm(1..B)$ , for some small bound

$B$ , which is, in some sense, a fusion between Pollard's rho and Pollard's p-1 methods.

# Pollard $p-1$ method

Pollard's p-1 method is based on Fermat's little theorem which states that:

$a^(p-1) ≡ 1 mod p$

for all prime numbers

$p$ with

$gcd(a, p) = 1$ .

Based on the identity

$(a^b)^c = a^(b c)$ , we have:

$(a^(p-1))^k ≡ a^((p-1) k) ≡ 1 mod p$

Therefore, if we can compute a multiple of

$p-1$ , then

$gcd((a^((p-1) k) mod n) - 1, n)$ may be a non-trivial factor of

$n$ , where

$p$ is a prime factor of

$n$ .

This method is effective when

$p-1$ has only small factors (i.e.:

$p-1$ is

$B$ -smooth, for some small bound

$B$ ).

Then

$(a^(B!) mod n) - 1$ is very likely to be a non-trivial factor of

$n$ , where

$a$ is any integer coprime to

$n$ .

A more optimized implementation, which does not require computing

$B!$ explicitly, is by setting:

$x_0 = 2$ (where

$2$ can be any other integer that is coprime to

$n$ )

then, for

$k >= 1$ , performing:

$x_k = x_(k-1)^k mod n$

at each step checking if

$gcd(x_k - 1, n)$ is a non-trivial factor of

$n$ .

Technically, we only need to be able to compute a multiple of the smallest divisor

$d$ of

$p-1$ , satisfying

$a^d ≡ 1 mod p$ . This value of

$d$ is called the multiplicative order of

$a$ modulo

$p$ .

Then

$gcd((a^(d m) mod n) - 1, n)$ is a divisor of

$n$ for all

$m >= 1$ .

# Fibonacci $p±1$ factorization method

The Fibonacci factorization method is based on the well-known property satisfied by all the prime numbers

$p$ :

$F_(p-k) ≡ 0 mod p$

where

$k = (5/p)$ is the Legendre symbol and

$F_n$ is the

$n$ -th Fibonacci number.

Therefore, if we can compute a multiple of the smallest divisor

$d$ of

$p - (5/p)$ , satisfying

$F_d ≡ 0 mod p$ , where

$p$ is a prime factor of

$n$ , then

$gcd(F_(d m) mod n, n)$ is a divisor of

$n$ for all

$m >= 1$ .

This method is similar in flavor to Pollard's p-1 and Williams' p+1 algorithms, and can efficiently find a prime factor

$p$ of

$n$ , if the smallest divisor

$d$ of

$p-(5/p)$ , satisfying

$F_d ≡ 0 mod p$ , is

$B$ -smooth for some small bound

$B$ .

# Lucas $p±1$ factorization method

A generalization of the Fibonacci factorization method, would be to use the Lucas sequence

$U_n(P,Q)$ , for which the following congruence is satisfied for all the prime numbers

$p$ :

$U_(p - k)(P,Q) ≡ 0 mod p$

where

$k = (D/p)$ is the Legendre symbol and

$D = P^2 - 4Q$ , with the condition that

$0 < P < p$ and

$Q < p$ and

$D != 0$ .

We can take advantage of this property by computing

$L = lcm(1..B)$ for a given small bound

$B$ , and iterating over

$P$ in a given range

$[a,b]$ , then checking if

$gcd(U_L(P,Q) mod n, n)$ is a non-trivial factor of

$n$ , where

$Q$ is a random integer value satisfying

$-n < Q < n$ with the condition that

$P^2 - 4Q != 0$ .

This method is effective when the smallest divisor

$d$ of

$p - (D/p)$ , satisfying

$U_d(P,Q) ≡ 0 mod p$ , is

$B$ -smooth for some small bound

$B$ .

If we can compute a multiple of

$d$ , then

$gcd(U_(d m)(P,Q) mod n, n)$ is a divisor of

$n$ for all

$m >= 1$ .

# Chebyshev $p±1$ factorization method

The Chebyshev factorization method uses the following nesting property of the Chebyshev polynomials

$T_n(x)$ :

$T_m(T_n(x)) = T_{m n}(x)$

By writing the Chebyshev polynomials

$T_n(x)$ in terms of the Lucas sequence

$V_n(P, Q)$ , we have:

$T_n(x) = \frac{V_n(2x, 1)}{2}$

This method is effective in finding a prime factor,

$p$ of

$n$ , when

$p-1$ or

$p+1$ is

$B$ -smooth, for some small value of

$B$ .

In factoring

$n$ , we start by setting an initial value for

$x_0$ , then for each

$k >= 1$ , we perform the following iteration:

$x_k = T_k(x_{k-1}) mod n$

at each step checking if

$gcd(x_k - 1, n)$ is a non-trivial factor

$n$ .

# Quadratic $p±1$ factorization method

The quadratic

$p±1$ factorization method is based on quadratic integers of the form

$a + b sqrt(w)$ , for some

$w >= 2$ .

By representing a quadratic integer

$a + b sqrt(w) = (a,b)$ , for some fixed

$w >= 2$ , we have the following multiplication properties:

$(a,b) * (c,d) = ((a c + b d w), (a d + b c))$

$(a,b)^2 = ((a^2 + b^2 w), (2 a b))$

This allows us to compute

$(a, b)^L mod n$ using exponentiation by squaring, where

$L = lcm(1..B)$ for some small bound

$B$ .

More efficiently, we can start by setting

$(a,b)$ to some random integers coprime to

$n$ , then, for each prime

$p$ in the range

$[2,B]$ , we perform the following iteration:

$(a,b) := (a,b)^(p^(floor(log_p(B)))) mod n$

at each step checking if either

$gcd(a-1, n)$ or

$gcd(b, n)$ is a non-trivial factor of

$n$ .

This method is effective when

$n$ has a prime factor

$p$ such that

$p-1$ or

$p+1$ is

$B$ -smooth.

# Sophie Germain's factorization identity

Sophie Germain's factorization identity is given as:

$x^4 + 4y^4 = (x^2 + 2xy + 2y^2) (x^2 - 2xy + 2y^2)$

This identity can also be used in factoring integers of the form

$a^4 + 4^(2b+1)$ , which can be written as

$x^4 + 4y^4$ , for

$x = a$ and

$y = 2^b$ , as illustrated bellow:

$a^4 + 4^(2b+1) = a^4 + 4 * 4^(2b) = a^4 + 4 * (2^b)^4$

The OEIS sequence of integers that can be represented as

$x^4 + 4y^4$ is A227855.

# Other factorization identities

Another special factorization identity can be created for integers of the form

$2^k (2^k - k + 1) + 1$ , for special values of

$k$ , similar in flavor to Sophie Germain's factorization identity.

For example, let:

$a(n) = 2^f(n) (2^f(n) - f(n) + 1) + 1$ , where

$f(n) = 4 n (n+1)$ .

Then

$a(n)$ can be factorized as:

$a(n) = (x - y - z + 1) (x + y + z + 1)$

where:

$x = 2^f(n)$

$y = 2^(f(n)/2)$

$z = n * 2^(f(n)/2 + 1)$

Based on this identity, we can also conclude that

$a(n)$ is composite for all

$n >= 1$ .

Also of interest may be the following identity (see: A231735):

$(2n+1)^2 (2k)^(2k) - 1 = (2^(k+1) k^k n + 2^k k^k - 1) (2^(k+1) k^k n + 2^k k^k + 1)$

# Implementations

This section includes source-code implementations of the presented methods (+ some additional methods), implemented in the Perl and Sidef programming languages.

Mathematics and computer science