# Метод квадратичного решета
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.")