182 lines
6.0 KiB
Python
182 lines
6.0 KiB
Python
# Метод квадратичного решета
|
|
import logging
|
|
import utils
|
|
from math import isqrt, sqrt, exp as expo, log, gcd
|
|
import time
|
|
|
|
# Configure logging
|
|
logging.basicConfig(level=logging.INFO,
|
|
format='%(asctime)s - %(message)s',
|
|
filename='quadratic_sieve.log',
|
|
filemode='w')
|
|
|
|
def iammain(): return __name__ == "__main__"
|
|
|
|
def SQ(n):
|
|
def gauss(M):
|
|
marks = [False] * len(M)
|
|
for j in range(len(M[0])):
|
|
if iammain(): print(f"[STEP_2] {j + 1}/{len(M[0])}")
|
|
for i in range(len(M)):
|
|
if M[i][j] == 1:
|
|
marks[i] = True
|
|
for k in range(j):
|
|
if M[i][k] == 1:
|
|
for row in range(len(M)):
|
|
M[row][k] = (M[row][k] + M[row][j]) % 2
|
|
for k in range(j + 1, len(M[0])):
|
|
if M[i][k] == 1:
|
|
for row in range(len(M)):
|
|
M[row][k] = (M[row][k] + M[row][j]) % 2
|
|
break
|
|
return marks, M
|
|
|
|
def get_dep_cols(row):
|
|
return [i for i, val in enumerate(row) if val == 1]
|
|
|
|
def row_add(new_row, current):
|
|
return [current[i] ^ M[new_row][i] for i in range(len(M[new_row]))]
|
|
|
|
def is_dependent(cols, row):
|
|
return any(row[i] == 1 for i in cols)
|
|
|
|
def find_linear_deps(row):
|
|
ret = []
|
|
dep_cols = get_dep_cols(M[row])
|
|
current_rows = [row]
|
|
current_sum = M[row][:]
|
|
for i in range(len(M)):
|
|
if i == row:
|
|
continue
|
|
if is_dependent(dep_cols, M[i]):
|
|
current_rows.append(i)
|
|
current_sum = row_add(i, current_sum)
|
|
if sum(current_sum) == 0:
|
|
ret.append(current_rows[:])
|
|
return ret
|
|
|
|
def testdep(dep):
|
|
x = y = 1
|
|
for row in dep:
|
|
x *= smooth_vals[row][0]
|
|
y *= smooth_vals[row][1]
|
|
s = x
|
|
t = isqrt(y)
|
|
if iammain(): logging.info(f"Found s and t such that s^2 = t^2 mod n: s = {s}, t = {t}")
|
|
return gcd(s - t, n)
|
|
|
|
def create_base(n, B):
|
|
base = []
|
|
i = 2
|
|
while len(base) < B:
|
|
if utils.legendre(n, i) == 1:
|
|
base.append(i)
|
|
i += 1
|
|
while not utils.is_prime(i):
|
|
i += 1
|
|
return base
|
|
|
|
def poly(x, a, b, n):
|
|
return ((a * x + b) ** 2) - n
|
|
|
|
def solve(a, b, n):
|
|
start_vals = []
|
|
for p in base:
|
|
ainv = 1
|
|
if a != 1:
|
|
g, ainv, _ = gcd(a, p)
|
|
assert g == 1
|
|
r1 = utils.tonelli(n, p)
|
|
r2 = (-1 * r1) % p
|
|
start1 = (ainv * (r1 - b)) % p
|
|
start2 = (ainv * (r2 - b)) % p
|
|
start_vals.append([start1, start2])
|
|
return start_vals
|
|
|
|
def trial(n, base):
|
|
ret = [0] * len(base)
|
|
if n > 0:
|
|
for i in range(len(base)):
|
|
while n % base[i] == 0:
|
|
n //= base[i]
|
|
ret[i] = (ret[i] + 1) % 2
|
|
return ret
|
|
|
|
N = n
|
|
a = 1
|
|
b = isqrt(N) + 1
|
|
bound = int(sqrt(expo(sqrt(log(n)*log(log(n))))))
|
|
base = create_base(N, bound)
|
|
needed = len(base) + 1
|
|
|
|
sieve_start = 0
|
|
sieve_stop = 0
|
|
sieve_interval = bound
|
|
|
|
M = []
|
|
smooth_vals = []
|
|
start_vals = solve(a, b, N)
|
|
seen = set()
|
|
|
|
if iammain():
|
|
logging.info(f"Number of elements in the base: {len(base)}")
|
|
logging.info(f"Last element in the base: {base[-1]}")
|
|
|
|
while len(smooth_vals) < needed:
|
|
sieve_start = sieve_stop
|
|
sieve_stop += sieve_interval
|
|
interval = [poly(x, a, b, N) for x in range(sieve_start, sieve_stop)]
|
|
|
|
for p in range(len(base)):
|
|
t = start_vals[p][0]
|
|
while start_vals[p][0] < sieve_start + sieve_interval:
|
|
while interval[start_vals[p][0] - sieve_start] % base[p] == 0:
|
|
interval[start_vals[p][0] - sieve_start] //= base[p]
|
|
start_vals[p][0] += base[p]
|
|
if start_vals[p][1] != t:
|
|
while start_vals[p][1] < sieve_start + sieve_interval:
|
|
while interval[start_vals[p][1] - sieve_start] % base[p] == 0:
|
|
interval[start_vals[p][1] - sieve_start] //= base[p]
|
|
start_vals[p][1] += base[p]
|
|
|
|
for i in range(sieve_interval):
|
|
if interval[i] == 1:
|
|
x = sieve_start + i
|
|
y = poly(x, a, b, N)
|
|
exp = trial(y, base)
|
|
if tuple(exp) not in seen:
|
|
if iammain(): print(f"[STEP_1] {len(smooth_vals)}/{needed}")
|
|
smooth_vals.append(((a * x) + b, y))
|
|
M.append(exp)
|
|
seen.add(tuple(exp))
|
|
|
|
if iammain():
|
|
logging.info(f"Used range of x values: from {sieve_start} to {sieve_stop}")
|
|
logging.info(f"Result of sieving: found {len(smooth_vals)} smooth numbers")
|
|
logging.info(f"Example of x and f(x) values: {smooth_vals[:5]}")
|
|
logging.info(f"Example of exponent vectors: {[v[:20] + (['...'] if len(v) > 20 else []) for v in M[:5]]}")
|
|
|
|
marks, M = gauss(M)
|
|
|
|
for i in range(len(marks)):
|
|
if iammain(): print(f"[STEP_3] {i + 1}/{len(marks)}")
|
|
if not marks[i]:
|
|
deps = find_linear_deps(i)
|
|
for dep in deps:
|
|
d = testdep(dep)
|
|
if d != 1 and d != N:
|
|
if iammain(): logging.info(f"Found non-trivial divisor: {d}")
|
|
return d
|
|
return None
|
|
|
|
if iammain():
|
|
N1 = 13611197472111783959 # takes 2 seconds
|
|
N2 = 1191515026104746183243378937330489098579 # does not compute
|
|
N3 = 74048093444435937986114388960912781233885985702403356033834092312625704192350369 # does not compute
|
|
|
|
number = N1
|
|
|
|
start_time = time.time()
|
|
SQ(number)
|
|
elapsed_time = time.time() - start_time
|
|
logging.info(f"Total program execution time: {elapsed_time} seconds.") |