QS method added

.gitignore

@@ -1,2 +1,3 @@
 __pycache__
-var61.txt
+var61.txt
+quadratic_sieve.log

QS.py

@@ -0,0 +1,180 @@
+# Метод квадратичного решета
+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 SQ(n):
+    def gauss(M):
+        marks = [False] * len(M)
+        for j in range(len(M[0])):
+            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)
+        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()
+
+    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:
+                    print(f"[STEP_1] {len(smooth_vals)}/{needed}")
+                    smooth_vals.append(((a * x) + b, y))
+                    M.append(exp)
+                    seen.add(tuple(exp))
+
+    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)):
+        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:
+                    logging.info(f"Found non-trivial divisor: {d}")
+                    return d
+    return None
+
+if __name__ == "__main__":
+    N1 = 13611197472111783959 # takes 2 seconds
+    N2 = 1191515026104746183243378937330489098579 # does not compute
+    N3 = 74048093444435937986114388960912781233885985702403356033834092312625704192350369 # does not compute
+
+    number = N1
+
+    start_time = time.time()
+    print('d=', SQ(number))
+    elapsed_time = time.time() - start_time
+    logging.info(f"Total program execution time: {elapsed_time} seconds.")

dixon.py

@@ -1,3 +1,4 @@
+# Алгоритм Диксона
 from random import randint
 from math import sqrt, exp, log, gcd
 import numpy as np

utils.py

@@ -1,5 +1,3 @@
-from collections import defaultdict
-
 def is_prime(N):
     """
     Тест Миллера-Рабина проверки числа на простоту.
@@ -69,4 +67,39 @@ def is_smooth(n, base):
             fact[i] += 1
     if n != 1:
         return None
-    return fact
+    return fact
+
+def legendre(a, p):
+    """Вычисление символа Лежандра"""
+    if a % p == 0:
+        return 0
+    return pow(a, (p - 1) // 2, p)
+
+def tonelli(n, p):
+    """Реализует алгоритм Тонелли-Шенкса для нахождения квадратного корня числа"""
+    q = p - 1
+    s = 0
+    while q % 2 == 0:
+        q //= 2
+        s += 1
+    if s == 1:
+        return pow(n, (p + 1) // 4, p)
+    z = 2
+    while legendre(z, p) != p - 1:
+        z += 1
+    c = pow(z, q, p)
+    r = pow(n, (q + 1) // 2, p)
+    t = pow(n, q, p)
+    m = s
+    while t != 1:
+        t2 = t
+        i = 0
+        while t2 != 1 and i < m:
+            t2 = pow(t2, 2, p)
+            i += 1
+        b = pow(c, 2 ** (m - i - 1), p)
+        r = (r * b) % p
+        c = (b * b) % p
+        t = (t * c) % p
+        m = i
+    return r