Naive Quadratic Sieve
The Quadratic Sieve is a factorization algorithm surpassed only by the General Number Field Sieve.
Here is a naive implementation in SageMath based on Eric Landquist’s lecture notes. It is the most basic version possible I could realize, it is actually painfully slow but easy to read and understand. Below are some improvements that should be considered as exercices by a student :
- Avoid exact arithmetic.
- Stop collecting relations when enough are found.
- Split the sieving interval.
- Implement Gaussian elminiation.
- Implement Tonelli-Shanks algorithm.
from tqdm import tqdm
def get_factor_base(n: int, bound: int) -> list[int]:
"""
Compute a factor base for `n` bounded by `bound`
Args:
n : Integer
bound : Integer
Returns:
The list of primes less than or equal to bound subject to n is a square mod p
"""
primes = Primes()
return [2] + [p for i in range(1, bound) if legendre_symbol(n, p := primes.unrank(i)) == 1]
def L(n: int)-> float:
"""
Compute the complexity function L[1/2,1] for `n` in L-notation
"""
return exp(sqrt(ln(n) * ln(ln(n))))
def smoothness_bound(n: int, scale: float = 0.5) -> int:
"""
Compute a smoothness bound for `n`, `scale` should me fine tuned to fit your memory constraints.
"""
return int(pow(L(n), scale))
def interval_bound(n: int, scale: float = 3 * sqrt(2) / 4) -> int:
"""
Compute an interval size for `n`, `scale` should be fine tuned to fit your memory constraints.
"""
return int(pow(L(n), scale))
def residues(x: int, p: int) -> list[int]:
"""
Compute quadratic residues of `x` mod `p`
"""
s = mod(x, p).sqrt()
return list(map(int,[s, p-s]))
def matrix_relations(smooth_numbers: list[int], factor_base: list[int]) -> sage.matrix.matrix_mod2_dense.Matrix_mod2_dense:
"""
Constructs a binary matrix representing factorization relations between a list of smooth numbers and a factor base.
This function generates a matrix where each row corresponds to a number in the `smooth_numbers` list,
and each column corresponds to a factor in the `factor_base` list. The entry M[i, j] is set to 1 if the
factor `factor_base[j]` divides `smooth_numbers[i]`, modulo 2.
Args:
smooth_numbers: List of integers whose factorizations are to be analyzed.
factor_base: List of prime factors serving as the reference factor base.
Returns:
The matrix where each row represents exponents mod 2 in the factorization of a number in `smooth_numbers` in terms of the `factor_base`.
"""
M = matrix(ZZ, len(smooth_numbers), len(factor_base))
for i, x in enumerate(smooth_numbers):
factorization = list(factor(x))
for p, np in factorization:
M[i, factor_base.index(p)] = np
return M.change_ring(GF(2))
def candidates(v, smooth_numbers: list[int], xis: list[int], s: int) -> list[int]:
"""
Computes candidate values (Px, Py) based on a binary vector and given lists of smooth numbers and xis.
Args:
v: A binary vector indicating which elements to include in the product.
smooth_numbers: List of integers representing the smooth numbers.
xis: List of integers representing the xis values.
s: The ceiling of square root of n.
Returns:
A list containing two integers: [isqrt(Px), Py], where:
- Px is the product of selected smooth numbers.
- Py is the product of (xis[i] + s) for selected indices.
- isqrt(Px) is the integer square root of Px.
"""
Px , Py = 1, 1
for i, bit in enumerate(v):
if bit:
Px *= smooth_numbers[i]
Py *= xis[i] + s
return [isqrt(Px), Py]
def check_candidates(test_x: int, test_y: int, n: int) -> list[int]:
"""
Checks if the candidate values (test_x, test_y) yield a non-trivial factorization of n.
This function verifies whether the difference between `test_x` and `test_y` shares a common factor with `n`.
If a non-trivial factorization of `n` is found, it returns the factors; otherwise, it returns an empty list.
Args:
test_x: An integer candidate value.
test_y: An integer candidate value.
n: The integer to be factored.
Returns:
A list containing two non-trivial factors of `n` if found, otherwise an empty list.
The factors satisfy: n = f1 * f2, where neither f1 nor f2 is 1.
"""
if test_x % n not in [test_y % n, -test_y % n]:
f1 = gcd(test_x-test_y, n)
f2 = n // f1
if n == f1 * f2 and (1 not in [f1, f2]):
return [f1, f2]
return []
def sieve(factor_base, M, n, s) -> list[int]:
"""
Performs a quadratic sieve to find smooth values of the form |(x + s)^2 - n| for x in [-M, M].
This function computes the values of |(x + s)^2 - n| for each x in the range [-M, M] and then
removes all prime factors in the `factor_base` from these values. The result is a list of partially
factored (smooth) values, which are candidates for further factorization of `n`.
Args:
factor_base: A list of prime numbers used as the factor base for sieving.
M: The range parameter, defining the interval [-M, M] for x values.
n: The integer to be factored.
s: The offset used in the quadratic expression (x + s)^2 - n.
Returns:
A list of integers representing the partially factored (smooth) values of |(x + s)^2 - n|
for x in [-M, M], after removing all factors in the `factor_base`.
"""
Qxis = [abs(pow(x + s, 2) - n) for x in range(-M, M + 1)]
for p in factor_base:
for r in set(residues(n, p)):
a = (r - s) % p
while a < M + 1:
Qxis[M + a] //= p
a += p
a = (r - s) % p - p
while a > - M :
Qxis[M + a] //= p
a -= p
return Qxis
def quadratic_sieve(n: int) -> list[int]:
B = smoothness_bound(n)
factor_base = get_factor_base(n, B)
M = interval_bound(n)
s = isqrt(n) + 1
Qxis = sieve(factor_base, M, n, s)
xis = [- M + i for i, Qx in enumerate(Qxis) if abs(Qx) == 1]
smooth_numbers = [abs(pow(x + s, 2) - n) for x in xis]
relations = matrix_relations(smooth_numbers, factor_base)
ker = relations.kernel()
for v in ker.basis():
test_x, test_y = candidates(v, smooth_numbers, xis, s)
if f := check_candidates(test_x, test_y, n):
return f
if __name__ == "__main__":
count = 0
BENCH_SIZE = 10**1
BIT_SIZE = 25
results = f""
for _ in tqdm(range(BENCH_SIZE)):
p = random_prime(2**BIT_SIZE, False, 2**(BIT_SIZE-1))
q = random_prime(2**BIT_SIZE, False, 2**(BIT_SIZE-1))
n = p * q
f = quadratic_sieve(n)
if f:
results += f"{n} = {f[0]} * {f[1]}\n"
count += 1
print(f"Rate of success : {float(count / BENCH_SIZE) * 100:.2f} %")
print("Factorisations found :\n", results, sep = "")