From 32e2e4d57ae768042a643cd2b028e49c9bf3c3b9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Mon, 7 Jun 2021 23:59:33 -0300
Subject: [PATCH 01/56] rename functions name to more suitable ones, fix
 minifactor3

---
 MPQS.sage | 34 +++++++++++++++++++---------------
 1 file changed, 19 insertions(+), 15 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 01d9355..117b248 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -139,8 +139,9 @@ def minifactor2(x, P, D = 1):
 
 def minifactor3(x, P, smooth_base):
     """ minifactor """
+    x = abs(x)
     smooth = gcd(x, smooth_base)
-    if x > smooth > 1:
+    if smooth > 1:
         not_smooth = x // smooth
         a1, r1, n1 = trial_division(smooth, P)
         a2, r2, n2 = trial_division(not_smooth, P)
@@ -158,7 +159,8 @@ def minifactor(x, P):
 
 
 class Poly:
-    """ Quadratic polynomial helper class """
+    """ Quadratic polynomial helper class: 
+    type Ax + 2Bx + C """
     def __init__(self, n, P, x_max, search = 0, verbose = None):
         self.n = int(n)
         self.P = P
@@ -242,10 +244,10 @@ class Poly:
         
 
 
-def rels_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = None):
+def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = None):
     """ relations search funcion """
     #print(N,start,stop)
-    sys.stderr.write("rels_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
+    sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
     if (stop-start) < 0:
         return [] 
 
@@ -262,6 +264,7 @@ def rels_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = No
     ltd = time.time()
     ld = len(Diffs)
     m = 1000
+    msg = ""
     for i in range(ld):
         if len(Rels) > required_relations:
             break
@@ -284,21 +287,22 @@ def rels_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = No
             eta = td * (ld/m)
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
-            sys.stderr.write("rels_find: range(%d, %d), inverval: %d of %d, found: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),tds,etas))
-    sys.stderr.write("rels_find: range(%d, %d, %d), interval: %d of %d, found: %d with prime_base: %d\n" % (start, stop,i, (stop-start),len(Rels),len(P)))
+            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),required_relations,tds,etas)
+            sys.stderr.write(msg)
+    sys.stderr.write(msg)
      
     #Rels += Found_Rels
     
     return Found_Rels
 
 
-def lin_alg(Rels, P):
+def linear_algebra(Rels, P):
     """ linear algebra, it generates a matrix in GF(2) then computes it's null-space"""
     M = matrix(GF(2), len(Rels), len(P), lambda i, j:P[j] in Rels[i][0][0])
     return M.left_kernel().basis()
 
 
-def proc_basis(N, basis, Rels):
+def process_basis_vectors(N, basis, Rels):
     """ process each basis vector """
     for K in basis:
         lhs = rhs = Ahs = 1
@@ -330,7 +334,7 @@ def find_primebase(n, bound):
     return primes, log_p
 
 
-def recalc_min_prime_thresh(thresh, Prime_base, log_p):
+def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
     """ recalculates the min efective prime base bound """
     min_prime = int(thresh * 3)
     fudge = sum(log_p[i] for i, p in enumerate(Prime_base) if p < min_prime)
@@ -340,7 +344,7 @@ def recalc_min_prime_thresh(thresh, Prime_base, log_p):
     return min_prime, thresh, fudge
 
 
-def gen_polys(N, Prime_base, x_max, needed):
+def generate_polys(N, Prime_base, x_max, needed):
     """ It searchs for distinct needed polys congruent to n."""
     n=1
     cpolys = 0
@@ -398,7 +402,7 @@ def _MPQS(N, verbose=True, M = 1):
     start = 0
     stop = B1
 
-    polys, early_factor = gen_polys(N, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
+    polys, early_factor = generate_polys(N, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
     if polys == None and early_factor != None:
         sys.stderr.write("Found small factor: %d\n" % early_factor)
         return [early_factor] + _MPQS(N // early_factor)
@@ -410,7 +414,7 @@ def _MPQS(N, verbose=True, M = 1):
         # trim primes, recalc min
         Prime_base = [p for p in Prime_base if p > min_prime]
         smooth_base = prod(Prime_base)
-        min_prime, thresh, fudge = recalc_min_prime_thresh(thresh, Prime_base, log_p)
+        min_prime, thresh, fudge = recalculate_min_prime_thresh(thresh, Prime_base, log_p)
 
         t1 = time.time()
         sys.stderr.write("Data collection with %d threads...\n" % T)
@@ -427,7 +431,7 @@ def _MPQS(N, verbose=True, M = 1):
         pols = []
         workpool = Pool(T)
         with workpool:
-            R = workpool.starmap(rels_find, inputs)  
+            R = workpool.starmap(relations_find, inputs)  
         workpool.close()
         workpool.join() 
         
@@ -440,12 +444,12 @@ def _MPQS(N, verbose=True, M = 1):
             sys.stderr.write("Found %d enough relations of %d needed relations...\n" % (len(Rels),required_relations))
             sys.stderr.write("Matrix creation...")
             
-            basis = lin_alg(Rels, Prime_base)
+            basis = linear_algebra(Rels, Prime_base)
 
             t3 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
             sys.stderr.write("Matrix reduction...\n")
-            result = proc_basis(N,basis, Rels)
+            result = process_basis_vectors(N,basis, Rels)
             t4 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t4-t3))
             

From cfe7095e4c17a0b7c445ae29e95e604029960077 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Tue, 8 Jun 2021 00:15:55 -0300
Subject: [PATCH 02/56] Add better comments

---
 MPQS.sage | 62 ++++++++++++++++++++++++++++++-------------------------
 1 file changed, 34 insertions(+), 28 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 117b248..85c8874 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -11,8 +11,9 @@ from multiprocessing import cpu_count, Pool, Manager
 from itertools import repeat
 import humanfriendly
 
+
 def mod_sqrt(n, p):
-    """ tonelli shanks algorithm """
+    """ Tonelli shanks algorithm """
     a = n % p
     if p % 4 == 3:
         return pow(a, (p+1) >> 2, p)
@@ -46,6 +47,7 @@ def mod_sqrt(n, p):
 
 
 def trial_division(n, P):
+    """ A simple trial division, factors are given by P-list. """
     a = []
     r = n
     l = len(P)
@@ -63,15 +65,14 @@ def trial_division(n, P):
     return a,r,n
 
 
-def filter_p(ppws):
-    """ filter out even powers """
+def filter_out_even_powers(ppws):
+    """ Filter out even powers. """
     d = {}
     for p,pw in ppws:
         if p not in d:
             d[p] = pw
         else:
             d[p] += pw
-    #print(d)
     tmp2 = []
     for p in d:
         if d[p] & 1 != 0: # only keep odd powers
@@ -81,7 +82,7 @@ def filter_p(ppws):
 
 def trial_division_minus_even_powers(n, P):
     a, r, n = trial_division(n, P)
-    a = filter_p(a)
+    a = filter_out_even_powers(a)
     return [a,r,n]
 
 
@@ -133,12 +134,12 @@ def minifactor2(x, P, D = 1):
         D//=x # reduce D
     #print(tmp)
     if R == 1: # if remainder is 1 then we found a B-smooth number then relation.
-        tmp = filter_p(a)  
+        tmp = filter_out_even_powers(a)  
         return (tmp, R, x), D                 
 
 
 def minifactor3(x, P, smooth_base):
-    """ minifactor """
+    """ Minifactor """
     x = abs(x)
     smooth = gcd(x, smooth_base)
     if smooth > 1:
@@ -147,12 +148,12 @@ def minifactor3(x, P, smooth_base):
         a2, r2, n2 = trial_division(not_smooth, P)
         if r2 == 1:
             a = a1 + a2
-            a3 = filter_p(a)
+            a3 = filter_out_even_powers(a)
             return a3, r2, x
         
                          
 def minifactor(x, P):
-    """ minifactor algo, finds odd-power primes in composites """
+    """ Minifactor algo, finds odd-power primes in composites """
     p = trial_division_minus_even_powers(x,P)
     if p[1] == 1: # if 1 x is B-smooth
         return p
@@ -171,6 +172,7 @@ class Poly:
         self.solve_for_x()
  
     def create(self):
+        """ Construct our polynomial. """
         n = self.n
         x_max = self.x_max
         root_2n = isqrt(2*n)
@@ -211,6 +213,7 @@ class Poly:
 
 
     def solve_for_x(self):
+        """ Get the mininmum x value of the polynomial where y=f(x) is a guaranteed relation. """
         A = self.A
         B = self.B
         C = self.C
@@ -234,7 +237,9 @@ class Poly:
                 self.start_vals = start_vals
         return start_vals
 
+
     def eval(self, x):
+        """ Eval the poly: y=f(x), return y and radical (Ax+B)"""
         A = self.A
         B = self.B
         C = self.C
@@ -243,16 +248,18 @@ class Poly:
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
         
 
-
 def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = None):
     """ relations search funcion """
-    #print(N,start,stop)
     sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
-    if (stop-start) < 0:
-        return [] 
 
-    I = [x for x in range(start, stop) if not is_prime(x) and not is_square(x)]
-   
+    #if (stop-start) < 0:
+    #    return [] 
+
+    if (stop-start) > 0:
+        I = [abs(x) for x in range(start, stop) if not is_prime(x) and not is_square(x)]
+    else: # this may not happen
+        I = [abs(x) for x in range(start, stop, -1) if not is_prime(x) and not is_square(x)]
+
     #D = reduce(lambda x, y: x * y, I)
     Diffs = [(pol.eval(x),x) for x in I]
     Found_Rels = []
@@ -273,13 +280,13 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
         #r = minifactor2(y, P, D//y)
         #f = minifactor(y, P)
         f = minifactor3(y, P, smooth_base)
-
         #print(r,trial_division(y,P),y)
         if 1:
         #if r != None:
             #f, D = r
             if f != None and f[1] == 1:
                 Rels.append((f,(y,Rad,A,x)))
+
         if i % m == 0:
             lt = time.time()
             td = lt - ltd
@@ -292,18 +299,17 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
     sys.stderr.write(msg)
      
     #Rels += Found_Rels
-    
     return Found_Rels
 
 
 def linear_algebra(Rels, P):
-    """ linear algebra, it generates a matrix in GF(2) then computes it's null-space"""
+    """ Linear algebra, it generates a matrix in GF(2) then computes it's left null-space. """
     M = matrix(GF(2), len(Rels), len(P), lambda i, j:P[j] in Rels[i][0][0])
     return M.left_kernel().basis()
 
 
 def process_basis_vectors(N, basis, Rels):
-    """ process each basis vector """
+    """ Process each basis vector, construct (a^2)-(b^2) mod n congruence. """
     for K in basis:
         lhs = rhs = Ahs = 1
         I = [f for f, k in zip(Rels, K) if k == 1]
@@ -319,7 +325,7 @@ def process_basis_vectors(N, basis, Rels):
 
 
 def find_primebase(n, bound):
-    """ finds the base prime for given n and bound """
+    """ Finds the base prime for given n and bound. """
     primes, mod_root, log_p, num_prime, p = [], [], [], 0, 3
     while p < bound or num_prime < 3:
         leg = legendre(n % p, p)
@@ -335,7 +341,7 @@ def find_primebase(n, bound):
 
 
 def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
-    """ recalculates the min efective prime base bound """
+    """ Recalculates the min efective prime base bound. """
     min_prime = int(thresh * 3)
     fudge = sum(log_p[i] for i, p in enumerate(Prime_base) if p < min_prime)
     fudge = fudge // 4
@@ -345,7 +351,7 @@ def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
 
 
 def generate_polys(N, Prime_base, x_max, needed):
-    """ It searchs for distinct needed polys congruent to n."""
+    """ It searchs for distinct needed polys congruent to n. """
     n=1
     cpolys = 0
     polys = []
@@ -369,7 +375,7 @@ def generate_polys(N, Prime_base, x_max, needed):
 
 
 def _MPQS(N, verbose=True, M = 1):
-    """ main MPQS function """
+    """ Main MPQS function. """
     bN, lN = int(log2(N)), len(str(N))
     i2N = isqrt(N) 
     i2Np1 = i2N + 1 
@@ -388,7 +394,6 @@ def _MPQS(N, verbose=True, M = 1):
         return Prime_base + _MPQS(N // Prime_base[0])
 
     x_max = B2 *60  # size of the sieve
-
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
     min_prime = int(thresh * 3)
@@ -400,7 +405,7 @@ def _MPQS(N, verbose=True, M = 1):
     sys.stderr.write("Need %d relations\n" % (required_relations))
 
     start = 0
-    stop = B1
+    stop = B1 # range to sieve
 
     polys, early_factor = generate_polys(N, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
     if polys == None and early_factor != None:
@@ -408,7 +413,7 @@ def _MPQS(N, verbose=True, M = 1):
         return [early_factor] + _MPQS(N // early_factor)
 
     manager = Manager()
-    Rels = manager.list()
+    Rels = manager.list() # placeholder list for relations shareable between child processes.
 
     while True:
         # trim primes, recalc min
@@ -425,7 +430,7 @@ def _MPQS(N, verbose=True, M = 1):
         for poly in polys:
             s1 = min(poly.start_vals[0]) 
             s2 = max(poly.start_vals[0]) 
-            inputs += [(N, start + s1, stop + s1 , Prime_base, smooth_base, Rels, required_relations,  poly)]
+            inputs += [(N, start + s1, stop + s2 , Prime_base, smooth_base, Rels, required_relations,  poly)]
 
         # deploy tasks to every cpu core.
         pols = []
@@ -444,7 +449,7 @@ def _MPQS(N, verbose=True, M = 1):
             sys.stderr.write("Found %d enough relations of %d needed relations...\n" % (len(Rels),required_relations))
             sys.stderr.write("Matrix creation...")
             
-            basis = linear_algebra(Rels, Prime_base)
+            basis = linear_algebra(Rels, Prime_base) # calculate left nullspace of a GF(2) matrix.
 
             t3 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
@@ -463,6 +468,7 @@ def _MPQS(N, verbose=True, M = 1):
       
 
 def MPQS(N):
+    """ Iterative version of MPQS. """
     result = _MPQS(N)
     if result != None:
         R= []

From f981729b89d629cf36c198034f27f182e2f22600 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 13:56:10 -0300
Subject: [PATCH 03/56] borrowed choose_multiplier from msieve

---
 MPQS.sage | 131 ++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 123 insertions(+), 8 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 85c8874..d3ea029 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -12,6 +12,106 @@ from itertools import repeat
 import humanfriendly
 
 
+def choose_multiplier(n, prime_list):
+    """
+    Code borrowed from msieve/mpqs.c
+    """
+    mult_list = [1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 26, 29, 30, 31, 33, 34, 35, 37, 38,39, 41, 42, 43, 46, 47, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73]
+    #mult_list = [1] + prime_list
+    #fb_size = len(P)
+    MAX_MP_WORDS = 30
+    NUM_MULTIPLIERS=len(mult_list)
+    NUM_TEST_PRIMES = 300
+
+    i =  j = 0
+    #num_primes = min(2 * fb_size, NUM_TEST_PRIMES);
+    best_score = 0.0
+    best_mult = 0
+    scores = [0.0] * NUM_MULTIPLIERS
+    num_multipliers = 0
+    #log2n = log(n)
+    M_LN2 = log(2)
+
+    #num_primes = min(num_primes, len(prime_list))
+    num_primes = len(prime_list)
+
+    """ measure the contribution of 2 as a factor of sieve
+       values. The multiplier itself must also be taken into
+       account in the score. scores[i] is the correction that
+       is implicitly applied to the size of sieve values for
+       multiplier i; a negative score makes sieve values 
+       smaller, and so is better """
+
+    for i in range(0, NUM_MULTIPLIERS): 
+        curr_mult = mult_list[i]
+        knmod8 = (curr_mult * n) % 8
+        logmult = log(curr_mult)
+        # only consider multipliers k such than
+        #   k*n will not overflow an mp_t */
+
+        #if (log2n + logmult > (32 * MAX_MP_WORDS - 2) * M_LN2)
+        #    break;
+
+        scores[i] = 0.5 * logmult;
+        if knmod8 == 1:
+            scores[i] -= 2 * M_LN2;
+            break
+        if knmod8 == 5:
+            scores[i] -= M_LN2;
+            break;
+        if knmod8 == 3 or knmod8 == 7:
+            scores[i] -= 0.5 * M_LN2;
+            break;
+        # even multipliers start with a handicap
+    num_multipliers = i;
+
+    # for the rest of the small factor base primes 
+
+    for i in range (1, num_primes):
+        prime = prime_list[i]
+        contrib = log(prime) / (prime - 1);
+        modp = n % prime
+
+        for j in range(0, num_multipliers):
+            curr_mult = mult_list[j];
+            knmodp = (curr_mult * modp) % prime
+            # if prime i is actually in the factor base
+            #   for k * n ... */
+            if (knmodp == 0 or legendre(knmodp, prime) == 1):
+
+                """ ...add its contribution. A prime p con-
+                   tributes log(p) to 1 in p sieve values, plus
+                   log(p) to 1 in p^2 sieve values, etc. The
+                   average contribution of all multiples of p 
+                   to a random sieve value is thus
+                   log(p) * (1/p + 1/p^2 + 1/p^3 + ...)
+                   = (log(p) / p) * 1 / (1 - (1/p)) 
+                   = log(p) / (p-1)
+                   This contribution occurs once for each
+                   square root used for sieving. There are two
+                   roots for each factor base prime, unless
+                   the prime divides k*n. In that case there 
+                   is only one root """
+
+                if (knmodp == 0):
+                    scores[j] -= contrib
+                else:
+                    scores[j] -= 2 * contrib
+
+    # use the multiplier that generates the best score 
+
+    best_score = 1000.0
+    best_mult = 1
+
+    #print(scores)
+   
+    for i in range(0, num_multipliers):
+        score = scores[i];
+        if (score < best_score):
+            best_score = score
+            best_mult = mult_list[i]
+    return best_mult;
+
 def mod_sqrt(n, p):
     """ Tonelli shanks algorithm """
     a = n % p
@@ -308,7 +408,7 @@ def linear_algebra(Rels, P):
     return M.left_kernel().basis()
 
 
-def process_basis_vectors(N, basis, Rels):
+def process_basis_vectors(N, basis, Rels, multiplier = 1):
     """ Process each basis vector, construct (a^2)-(b^2) mod n congruence. """
     for K in basis:
         lhs = rhs = Ahs = 1
@@ -320,8 +420,17 @@ def process_basis_vectors(N, basis, Rels):
                 Ahs *= i[1][2] # A-term in poly
             LHS = Ahs * lhs
             g = gcd(isqrt(LHS)-rhs,N)
-            if N > g > 1:
-                return [int(g), int(N//g)]
+            if N > g > 1:     
+                factors = [int(g), int(N//g)]
+                if multiplier > 1:
+                    tmp2 = []
+                    for factor in factors
+                        ignoreme  = gcd(factor,multiplier)
+                        if r > ignoreme > 1:
+                            tmp2.append(factor//ignoreme)
+                    return tmp2
+                else:
+                    return tmp
 
 
 def find_primebase(n, bound):
@@ -393,6 +502,11 @@ def _MPQS(N, verbose=True, M = 1):
         sys.stderr.write("Found small factor: %d\n" % Prime_base[0])
         return Prime_base + _MPQS(N // Prime_base[0])
 
+    multiplier = choose_multiplier(N, Prime_base)
+    Nm = multiplier * N
+    
+    sys.stderr.write("Multiplier is: %d" % multiplier)
+
     x_max = B2 *60  # size of the sieve
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
@@ -407,10 +521,11 @@ def _MPQS(N, verbose=True, M = 1):
     start = 0
     stop = B1 # range to sieve
 
-    polys, early_factor = generate_polys(N, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
+    polys, early_factor = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
     if polys == None and early_factor != None:
-        sys.stderr.write("Found small factor: %d\n" % early_factor)
-        return [early_factor] + _MPQS(N // early_factor)
+        if multiplier > 1 and gcd(early_factor, multiplier) == 1:
+            sys.stderr.write("Found small factor: %d\n" % early_factor)
+            return [early_factor] + _MPQS(N // early_factor)
 
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
@@ -430,7 +545,7 @@ def _MPQS(N, verbose=True, M = 1):
         for poly in polys:
             s1 = min(poly.start_vals[0]) 
             s2 = max(poly.start_vals[0]) 
-            inputs += [(N, start + s1, stop + s2 , Prime_base, smooth_base, Rels, required_relations,  poly)]
+            inputs += [(Nm, start + s1, stop + s2 , Prime_base, smooth_base, Rels, required_relations,  poly)]
 
         # deploy tasks to every cpu core.
         pols = []
@@ -454,7 +569,7 @@ def _MPQS(N, verbose=True, M = 1):
             t3 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
             sys.stderr.write("Matrix reduction...\n")
-            result = process_basis_vectors(N,basis, Rels)
+            result = process_basis_vectors(Nm, basis, Rels, multiplier)
             t4 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t4-t3))
             

From 649811a6b785d7ae9336ec3180a5c2b136032e35 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 14:00:16 -0300
Subject: [PATCH 04/56] Update docstrings

---
 MPQS.sage | 24 ++++++++++++++++++++----
 1 file changed, 20 insertions(+), 4 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index d3ea029..7e87641 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -272,7 +272,11 @@ class Poly:
         self.solve_for_x()
  
     def create(self):
-        """ Construct our polynomial. """
+        """ 
+        Construct our polynomial. 
+        Code borrowed from 
+        https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
+        """
         n = self.n
         x_max = self.x_max
         root_2n = isqrt(2*n)
@@ -313,7 +317,11 @@ class Poly:
 
 
     def solve_for_x(self):
-        """ Get the mininmum x value of the polynomial where y=f(x) is a guaranteed relation. """
+        """ 
+        Get the mininmum x value of the polynomial where y=f(x) is a guaranteed relation. 
+        Code borrowed https://github.com/cramppet/quadratic-sieve/blob/master/QS.py
+        """
+        
         A = self.A
         B = self.B
         C = self.C
@@ -434,7 +442,11 @@ def process_basis_vectors(N, basis, Rels, multiplier = 1):
 
 
 def find_primebase(n, bound):
-    """ Finds the base prime for given n and bound. """
+    """ 
+    Finds the base prime for given n and bound. 
+    https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
+    
+    """
     primes, mod_root, log_p, num_prime, p = [], [], [], 0, 3
     while p < bound or num_prime < 3:
         leg = legendre(n % p, p)
@@ -450,7 +462,11 @@ def find_primebase(n, bound):
 
 
 def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
-    """ Recalculates the min efective prime base bound. """
+    """ 
+    Recalculates the min efective prime base bound. 
+    code borrowed from:
+    https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
+    """
     min_prime = int(thresh * 3)
     fudge = sum(log_p[i] for i, p in enumerate(Prime_base) if p < min_prime)
     fudge = fudge // 4

From e83c319ddd1b2118b87a474e9cc436f50a516041 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 14:28:18 -0300
Subject: [PATCH 05/56] better handling of early factorizations

---
 MPQS.sage | 53 ++++++++++++++++++++++++++++-------------------------
 1 file changed, 28 insertions(+), 25 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 7e87641..f2e7312 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -268,6 +268,7 @@ class Poly:
         self.x_max = x_max
         self.search = search
         self.verbose = verbose
+        self.early_factors = []
         self.create()
         self.solve_for_x()
  
@@ -283,18 +284,15 @@ class Poly:
         root_A = next_prime(isqrt(root_2n // x_max))
         self.root_A = root_A
         s=0
-        while True:
+        leg = 2
+        while leg != 1:
             root_A = next_prime(root_A)
             leg = legendre(n, root_A)
             if leg == 1:
                 if s > self.search:
                     break
             elif leg == 0:
-                self.early_factor = root_A
-                self.A = None
-                self.B = None
-                self.C = None
-                return root_A, None, None
+                self.early_factors.append(root_A)
             s+=1
 
         self.s = s
@@ -481,22 +479,23 @@ def generate_polys(N, Prime_base, x_max, needed):
     cpolys = 0
     polys = []
     polys_ABC = []
-    early_factor = None
+    early_factors = []
     while cpolys <= needed:
         pol = Poly(N, Prime_base, x_max, search=n, verbose=False)
-        if pol.A == None and pol.B == None and pol.C == None:
-            early_factor = pol.early_factor
-            polys = None
-            break
+        if len(pol.early_factors) > 0:
+            for early_factor in pol.early_factors:
+                if early_factor not in early_factors:
+                    early_factors.append(early_factor)
+        else:
+            pol_ABC = (pol.A,pol.B,pol.C)
+            if pol_ABC not in polys_ABC:
+                m = "Found poly: f(x) = %dx^2 + %dx + %d\n" % (pol_ABC)
+                sys.stderr.write(m.replace("+ -","- "))
+                polys_ABC.append(pol_ABC)
+                polys.append(pol)
+                cpolys += 1
         n += 1
-        pol_ABC = (pol.A,pol.B,pol.C)
-        if pol_ABC not in polys_ABC:
-            m = "Found poly: f(x) = %dx^2 + %dx + %d\n" % (pol_ABC)
-            sys.stderr.write(m.replace("+ -","- "))
-            polys_ABC.append(pol_ABC)
-            polys.append(pol)
-            cpolys += 1
-    return polys, early_factor
+    return polys, early_factors
 
 
 def _MPQS(N, verbose=True, M = 1):
@@ -537,12 +536,16 @@ def _MPQS(N, verbose=True, M = 1):
     start = 0
     stop = B1 # range to sieve
 
-    polys, early_factor = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
-    if polys == None and early_factor != None:
-        if multiplier > 1 and gcd(early_factor, multiplier) == 1:
-            sys.stderr.write("Found small factor: %d\n" % early_factor)
-            return [early_factor] + _MPQS(N // early_factor)
-
+    polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
+    if (early_factors) > 0:
+        tmp = 1
+        for early_factor in early_factors:
+            g = gcd(early_factor, N)    
+            if N > g > 1: 
+                tmp *= g
+         if tmp > 1:
+             return [tmp, N // tmp]
+    
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
 

From 6360a1b012bf058fa1b09a983185fdc1fb8bacf2 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 14:31:42 -0300
Subject: [PATCH 06/56] better handling of early factorizations

---
 MPQS.sage | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index f2e7312..da1cb33 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -539,12 +539,14 @@ def _MPQS(N, verbose=True, M = 1):
     polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
     if (early_factors) > 0:
         tmp = 1
+        small = []
         for early_factor in early_factors:
             g = gcd(early_factor, N)    
             if N > g > 1: 
                 tmp *= g
+                small.append(g)
          if tmp > 1:
-             return [tmp, N // tmp]
+             return small + MPQS(N // tmp)
     
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.

From 10f5ca05c16e5ecf18f5ec6d9edb44501093cf4f Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 14:33:10 -0300
Subject: [PATCH 07/56] better handling of early factorizations

---
 MPQS.sage | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index da1cb33..6e132f8 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -436,7 +436,7 @@ def process_basis_vectors(N, basis, Rels, multiplier = 1):
                             tmp2.append(factor//ignoreme)
                     return tmp2
                 else:
-                    return tmp
+                    return factors
 
 
 def find_primebase(n, bound):

From f10db1a0e0eedc20bfc02387a070c94331ea373e Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 14:34:27 -0300
Subject: [PATCH 08/56] some fixes

---
 MPQS.sage | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 6e132f8..b6f5d58 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -430,7 +430,7 @@ def process_basis_vectors(N, basis, Rels, multiplier = 1):
                 factors = [int(g), int(N//g)]
                 if multiplier > 1:
                     tmp2 = []
-                    for factor in factors
+                    for factor in factors:
                         ignoreme  = gcd(factor,multiplier)
                         if r > ignoreme > 1:
                             tmp2.append(factor//ignoreme)
@@ -539,14 +539,14 @@ def _MPQS(N, verbose=True, M = 1):
     polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
     if (early_factors) > 0:
         tmp = 1
-        small = []
+         small = []
         for early_factor in early_factors:
             g = gcd(early_factor, N)    
             if N > g > 1: 
                 tmp *= g
                 small.append(g)
-         if tmp > 1:
-             return small + MPQS(N // tmp)
+        if tmp > 1:
+            return small + MPQS(N // tmp)
     
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.

From a50f44f5fe774476dbe4bc5a9e417a0cf687fc6a Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 15:00:26 -0300
Subject: [PATCH 09/56] Update, some fixes

---
 MPQS.sage | 84 +++++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 60 insertions(+), 24 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index b6f5d58..6e96553 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -68,7 +68,7 @@ def choose_multiplier(n, prime_list):
     # for the rest of the small factor base primes 
 
     for i in range (1, num_primes):
-        prime = prime_list[i]
+        prime = int(prime_list[i])
         contrib = log(prime) / (prime - 1);
         modp = n % prime
 
@@ -113,7 +113,9 @@ def choose_multiplier(n, prime_list):
     return best_mult;
 
 def mod_sqrt(n, p):
-    """ Tonelli shanks algorithm """
+    """ 
+    Tonelli shanks algorithm 
+    """
     a = n % p
     if p % 4 == 3:
         return pow(a, (p+1) >> 2, p)
@@ -147,7 +149,9 @@ def mod_sqrt(n, p):
 
 
 def trial_division(n, P):
-    """ A simple trial division, factors are given by P-list. """
+    """ 
+    A simple trial division, factors are given by P-list. 
+    """
     a = []
     r = n
     l = len(P)
@@ -166,7 +170,9 @@ def trial_division(n, P):
 
 
 def filter_out_even_powers(ppws):
-    """ Filter out even powers. """
+    """ 
+    Filter out even powers. 
+    """
     d = {}
     for p,pw in ppws:
         if p not in d:
@@ -181,13 +187,18 @@ def filter_out_even_powers(ppws):
 
 
 def trial_division_minus_even_powers(n, P):
+    """
+    Factor a composite n returning only odd power primes.
+    """
     a, r, n = trial_division(n, P)
     a = filter_out_even_powers(a)
     return [a,r,n]
 
 
 def is_smooth(x, P):
-    """ check if n is B-smooth """
+    """ 
+    check if n is B-smooth. 
+    """
     y = x
     for p in P:
         while y % p ==0:
@@ -196,7 +207,9 @@ def is_smooth(x, P):
 
 
 def minifactor2(x, P, D = 1):
-    """ tries to factor out common primes with gcd and filter out even power primes"""
+    """ 
+    Tries to factor out common primes with gcd and filter out even power primes. 
+    """
     tmp = []
     R = 1
     g = gcd(x, D) 
@@ -239,7 +252,9 @@ def minifactor2(x, P, D = 1):
 
 
 def minifactor3(x, P, smooth_base):
-    """ Minifactor """
+    """ 
+    Minifactor. 
+    """
     x = abs(x)
     smooth = gcd(x, smooth_base)
     if smooth > 1:
@@ -253,15 +268,19 @@ def minifactor3(x, P, smooth_base):
         
                          
 def minifactor(x, P):
-    """ Minifactor algo, finds odd-power primes in composites """
+    """ 
+    Minifactor algo, finds odd-power primes in composites. 
+    """
     p = trial_division_minus_even_powers(x,P)
     if p[1] == 1: # if 1 x is B-smooth
         return p
 
 
 class Poly:
-    """ Quadratic polynomial helper class: 
-    type Ax + 2Bx + C """
+    """ 
+    Quadratic polynomial helper class: 
+    type Ax + 2Bx + C 
+    """
     def __init__(self, n, P, x_max, search = 0, verbose = None):
         self.n = int(n)
         self.P = P
@@ -285,7 +304,7 @@ class Poly:
         self.root_A = root_A
         s=0
         leg = 2
-        while leg != 1:
+        while True:
             root_A = next_prime(root_A)
             leg = legendre(n, root_A)
             if leg == 1:
@@ -293,7 +312,10 @@ class Poly:
                     break
             elif leg == 0:
                 self.early_factors.append(root_A)
-            s+=1
+                if s > self.search:
+                    break
+                #break
+            s += 1
 
         self.s = s
         A = int(pow(root_A, 2))
@@ -345,7 +367,9 @@ class Poly:
 
 
     def eval(self, x):
-        """ Eval the poly: y=f(x), return y and radical (Ax+B)"""
+        """ 
+        Eval the poly: y=f(x), return y and radical (Ax+B).
+        """
         A = self.A
         B = self.B
         C = self.C
@@ -355,7 +379,9 @@ class Poly:
         
 
 def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = None):
-    """ relations search funcion """
+    """ 
+    Relations search funcion 
+    """
     sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
 
     #if (stop-start) < 0:
@@ -409,13 +435,17 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
 
 
 def linear_algebra(Rels, P):
-    """ Linear algebra, it generates a matrix in GF(2) then computes it's left null-space. """
+    """ 
+    Linear algebra, it generates a matrix in GF(2) then computes it's left null-space. 
+    """
     M = matrix(GF(2), len(Rels), len(P), lambda i, j:P[j] in Rels[i][0][0])
     return M.left_kernel().basis()
 
 
 def process_basis_vectors(N, basis, Rels, multiplier = 1):
-    """ Process each basis vector, construct (a^2)-(b^2) mod n congruence. """
+    """ 
+    Process each basis vector, construct (a^2)-(b^2) mod n congruence. 
+    """
     for K in basis:
         lhs = rhs = Ahs = 1
         I = [f for f, k in zip(Rels, K) if k == 1]
@@ -474,14 +504,17 @@ def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
 
 
 def generate_polys(N, Prime_base, x_max, needed):
-    """ It searchs for distinct needed polys congruent to n. """
+    """ 
+    It searchs for distinct needed polys congruent to n. 
+    """
     n=1
     cpolys = 0
     polys = []
     polys_ABC = []
     early_factors = []
     while cpolys <= needed:
-        pol = Poly(N, Prime_base, x_max, search=n, verbose=False)
+        pol = Poly(N, Prime_base, x_max, search = n, verbose=False)
+        #print(pol.early_factors)
         if len(pol.early_factors) > 0:
             for early_factor in pol.early_factors:
                 if early_factor not in early_factors:
@@ -499,7 +532,9 @@ def generate_polys(N, Prime_base, x_max, needed):
 
 
 def _MPQS(N, verbose=True, M = 1):
-    """ Main MPQS function. """
+    """ 
+    Main MPQS function. 
+    """
     bN, lN = int(log2(N)), len(str(N))
     i2N = isqrt(N) 
     i2Np1 = i2N + 1 
@@ -520,8 +555,6 @@ def _MPQS(N, verbose=True, M = 1):
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
     
-    sys.stderr.write("Multiplier is: %d" % multiplier)
-
     x_max = B2 *60  # size of the sieve
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
@@ -531,15 +564,16 @@ def _MPQS(N, verbose=True, M = 1):
     required_relations = round(len(Prime_base) * required_relations_ratio) 
 
     sys.stderr.write("Factoring N: %d, bits: %d, digits: %d, B1: %d, B2: %d\n" % (N,bN,lN,B1,B2))
+    sys.stderr.write("Multiplier is: %d\n" % multiplier)
     sys.stderr.write("Need %d relations\n" % (required_relations))
 
     start = 0
     stop = B1 # range to sieve
 
     polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
-    if (early_factors) > 0:
+    if len(early_factors) > 0:
         tmp = 1
-         small = []
+        small = []
         for early_factor in early_factors:
             g = gcd(early_factor, N)    
             if N > g > 1: 
@@ -604,7 +638,9 @@ def _MPQS(N, verbose=True, M = 1):
       
 
 def MPQS(N):
-    """ Iterative version of MPQS. """
+    """ 
+    Iterative version of MPQS. 
+    """
     result = _MPQS(N)
     if result != None:
         R= []

From 5814d44fb8f403b343b24c4050eca9e5427ae06b Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 17:03:23 -0300
Subject: [PATCH 10/56] add linear algebra, Gauss_elimination_GF2,
 matrix_transpose

---
 MPQS.sage | 36 +++++++++++++++++++++++++++++++++++-
 1 file changed, 35 insertions(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index 6e96553..ddab7e0 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -433,6 +433,40 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
     #Rels += Found_Rels
     return Found_Rels
 
+def transpose(A):
+    """
+    Transpose matrix so columns become rows
+    """
+    new_A = []
+    for i in range(len(A[0])):
+        new_row = []
+        for row in A:
+            new_row.append(row[i])
+        new_A.append(new_row)
+    return(new_A)
+
+def Gaussian_elimination_GF2(A):
+  h = len(A)
+  m = len(A[0])
+  marks = [False] * h
+  for j in range(0,m):
+    for i in range(0,h):
+      if A[i][j] == 1:
+          marks[i] = True
+          for k in range(j+1,m):
+            if A[i][k] == 1:
+              A[i][k] = (A[i][j] ^ A[i][k]) 
+          break
+  return marks, A
+
+def left_nullspace(A):
+    """
+    Compute left null space:
+    x such that xTA = 0T
+    """
+    A = transpose(A)
+    marks, A = Gaussian_elimination_GF2(A)
+    return marks, A
 
 def linear_algebra(Rels, P):
     """ 
@@ -473,7 +507,7 @@ def find_primebase(n, bound):
     """ 
     Finds the base prime for given n and bound. 
     https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
-    
+    Same as the sieve of erathostenes.
     """
     primes, mod_root, log_p, num_prime, p = [], [], [], 0, 3
     while p < bound or num_prime < 3:

From c8690491bb225930953c44453737b772639e0cf9 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 17:14:29 -0300
Subject: [PATCH 11/56] update docstring

---
 MPQS.sage | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index ddab7e0..c8233f3 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -446,6 +446,10 @@ def transpose(A):
     return(new_A)
 
 def Gaussian_elimination_GF2(A):
+  """
+  Gaussian elimination over GF2
+  https://www.cs.umd.edu/~gasarch/TOPICS/factoring/fastgauss.pdf, page2
+  """
   h = len(A)
   m = len(A[0])
   marks = [False] * h
@@ -634,7 +638,7 @@ def _MPQS(N, verbose=True, M = 1):
         for poly in polys:
             s1 = min(poly.start_vals[0]) 
             s2 = max(poly.start_vals[0]) 
-            inputs += [(Nm, start + s1, stop + s2 , Prime_base, smooth_base, Rels, required_relations,  poly)]
+            inputs += [(Nm, start + s1, stop + s1 , Prime_base, smooth_base, Rels, required_relations,  poly)]
 
         # deploy tasks to every cpu core.
         pols = []

From b53779c3efc77044f8f6a2625dc0971bfd6498c1 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Tue, 8 Jun 2021 19:56:47 -0300
Subject: [PATCH 12/56] add partial relations merging

---
 MPQS.sage | 67 ++++++++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 51 insertions(+), 16 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index c8233f3..ef62101 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -168,6 +168,17 @@ def trial_division(n, P):
             i += 1
     return a,r,n
 
+def merge_powers(ppws):
+    d = {}
+    for p,pw in ppws:
+        if p not in d:
+            d[p] = pw
+        else:
+            d[p] += pw
+    tmp2 = []
+    for p in d:
+       tmp2.append((p,d[p]))
+    return tmp2
 
 def filter_out_even_powers(ppws):
     """ 
@@ -265,6 +276,23 @@ def minifactor3(x, P, smooth_base):
             a = a1 + a2
             a3 = filter_out_even_powers(a)
             return a3, r2, x
+
+
+def minifactor4(x, P, smooth_base):
+    """ 
+    Minifactor. 
+    """
+    x = abs(x)
+    a1 = []
+    smooth = gcd(x, smooth_base)
+    if smooth > 1:
+        not_smooth = x // smooth
+        a1, r1, n1 = trial_division(smooth, P)
+        a2, r2, n2 = trial_division(not_smooth, P)
+    else:
+        a2, r2, n2 = trial_division(x, P) 
+    return merge_powers(a1 + a2), r2, x
+
         
                          
 def minifactor(x, P):
@@ -404,21 +432,26 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
     ld = len(Diffs)
     m = 1000
     msg = ""
+    partials = {}
+    merged = 0
     for i in range(ld):
         if len(Rels) > required_relations:
             break
         yRad, x = Diffs[i]
         y, Rad = yRad
-        #r = minifactor2(y, P, D//y)
-        #f = minifactor(y, P)
-        f = minifactor3(y, P, smooth_base)
-        #print(r,trial_division(y,P),y)
+        f = minifactor4(y, P, smooth_base)
         if 1:
-        #if r != None:
-            #f, D = r
-            if f != None and f[1] == 1:
-                Rels.append((f,(y,Rad,A,x)))
-
+            if f != None:
+                filtered = (filter_out_even_powers(f[0]),f[1],y)
+                if f[1] == 1:    
+                    Rels.append([f[0], y, Rad, A])
+                elif f[1] in partials:
+                    a = partials[f[1]]
+                    p = filter_out_even_powers(f[0] + a[0])
+                    Rels.append([p, y * a[1], Rad * a[2], A* a[3]])
+                    merged += 1
+                else:
+                    partials[f[1]] = [f[0], y, Rad, A]
         if i % m == 0:
             lt = time.time()
             td = lt - ltd
@@ -426,10 +459,12 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
             eta = td * (ld/m)
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
-            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),required_relations,tds,etas)
+            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged,tds,etas)
             sys.stderr.write(msg)
+    
+    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged)
     sys.stderr.write(msg)
-     
+    
     #Rels += Found_Rels
     return Found_Rels
 
@@ -476,7 +511,7 @@ def linear_algebra(Rels, P):
     """ 
     Linear algebra, it generates a matrix in GF(2) then computes it's left null-space. 
     """
-    M = matrix(GF(2), len(Rels), len(P), lambda i, j:P[j] in Rels[i][0][0])
+    M = matrix(GF(2), len(Rels), len(P), lambda i, j:P[j] in Rels[i][0])
     return M.left_kernel().basis()
 
 
@@ -489,11 +524,11 @@ def process_basis_vectors(N, basis, Rels, multiplier = 1):
         I = [f for f, k in zip(Rels, K) if k == 1]
         if len(I) > 0:
             for i in I:
-                lhs *= i[1][0] # left-hand side
-                rhs *= i[1][1] # rigth-hand side
-                Ahs *= i[1][2] # A-term in poly
+                lhs *= i[1] # left-hand side
+                rhs *= i[2] # rigth-hand side
+                Ahs *= i[3] # A-term in poly
             LHS = Ahs * lhs
-            g = gcd(isqrt(LHS)-rhs,N)
+            g = gcd(isqrt(abs(LHS))-rhs,N)
             if N > g > 1:     
                 factors = [int(g), int(N//g)]
                 if multiplier > 1:

From 970000663b14c2da1c35104ac192deaa746c92ea Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Tue, 8 Jun 2021 20:01:18 -0300
Subject: [PATCH 13/56] Update MPQS.sage

---
 MPQS.sage | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index ef62101..9ffc294 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -3,6 +3,16 @@
 # Author Dario Clavijo 2021
 # License: GPLv3
 
+"""
+The quadratic sieve algorithm (QS) is an integer factorization algorithm and, in practice, 
+the second fastest method known (after the general number field sieve). 
+It is still the fastest for integers under 100 decimal digits or so, and is considerably simpler than the number field sieve. 
+It is a general-purpose factorization algorithm, meaning that its running time depends solely on the size of the integer to be factored,
+and not on special structure or properties. 
+It was invented by Carl Pomerance in 1981 as an improvement to Schroeppel's linear sieve
+https://en.wikipedia.org/wiki/Quadratic_sieve
+"""
+
 import time
 import sys
 from gmpy2 import gcd, gcdext, isqrt, is_prime, next_prime, log2, log10, legendre, powmod, invert
@@ -16,23 +26,18 @@ def choose_multiplier(n, prime_list):
     """
     Code borrowed from msieve/mpqs.c
     """
-    mult_list = [1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 26, 29, 30, 31, 33, 34, 35, 37, 38,39, 41, 42, 43, 46, 47, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73]
-    #mult_list = [1] + prime_list
-    #fb_size = len(P)
+    mult_list = [1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 26, 29, 30, 31, 33, 34, 35, 37, 38]
+    mult_list += [39, 41, 42, 43, 46, 47, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73]
+
     MAX_MP_WORDS = 30
     NUM_MULTIPLIERS=len(mult_list)
     NUM_TEST_PRIMES = 300
-
     i =  j = 0
-    #num_primes = min(2 * fb_size, NUM_TEST_PRIMES);
     best_score = 0.0
     best_mult = 0
     scores = [0.0] * NUM_MULTIPLIERS
     num_multipliers = 0
-    #log2n = log(n)
     M_LN2 = log(2)
-
-    #num_primes = min(num_primes, len(prime_list))
     num_primes = len(prime_list)
 
     """ measure the contribution of 2 as a factor of sieve
@@ -49,9 +54,6 @@ def choose_multiplier(n, prime_list):
         # only consider multipliers k such than
         #   k*n will not overflow an mp_t */
 
-        #if (log2n + logmult > (32 * MAX_MP_WORDS - 2) * M_LN2)
-        #    break;
-
         scores[i] = 0.5 * logmult;
         if knmod8 == 1:
             scores[i] -= 2 * M_LN2;

From 57c7de454cdb87213e4b59271ca56a99f8fa4f7b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Tue, 8 Jun 2021 20:18:49 -0300
Subject: [PATCH 14/56] Fix potential bug

---
 MPQS.sage | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index 9ffc294..5f8bce7 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -446,7 +446,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
             if f != None:
                 filtered = (filter_out_even_powers(f[0]),f[1],y)
                 if f[1] == 1:    
-                    Rels.append([f[0], y, Rad, A])
+                    Rels.append([filtered, y, Rad, A])
                 elif f[1] in partials:
                     a = partials[f[1]]
                     p = filter_out_even_powers(f[0] + a[0])

From b19d1b8faf81e5a9e2ce5ac2c5370f67925bd036 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 11:13:37 -0300
Subject: [PATCH 15/56] update MPQS.sage

---
 MPQS.sage | 26 +++++++++++++++-----------
 1 file changed, 15 insertions(+), 11 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 5f8bce7..7ac9cbc 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -408,7 +408,7 @@ class Poly:
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
         
 
-def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol = None):
+def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_relations, pol = None):
     """ 
     Relations search funcion 
     """
@@ -425,17 +425,16 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
     #D = reduce(lambda x, y: x * y, I)
     Diffs = [(pol.eval(x),x) for x in I]
     Found_Rels = []
-
     A = pol.A
     #B = pol.B
     #C = pol.C
-
-    ltd = time.time()
+    st = time.time
+    ltd = st
     ld = len(Diffs)
     m = 1000
     msg = ""
     partials = {}
-    merged = 0
+
     for i in range(ld):
         if len(Rels) > required_relations:
             break
@@ -451,7 +450,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
                     a = partials[f[1]]
                     p = filter_out_even_powers(f[0] + a[0])
                     Rels.append([p, y * a[1], Rad * a[2], A* a[3]])
-                    merged += 1
+                    merged_count += 1
                 else:
                     partials[f[1]] = [f[0], y, Rad, A]
         if i % m == 0:
@@ -461,10 +460,12 @@ def relations_find(N, start, stop, P, smooth_base, Rels, required_relations, pol
             eta = td * (ld/m)
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
-            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged,tds,etas)
+            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged_count,tds,etas)
             sys.stderr.write(msg)
     
-    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged)
+    td = time.time() - st
+    tds = humanfriendly.format_timespan(td)
+    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, time elapsed: %s\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged_count, tds)
     sys.stderr.write(msg)
     
     #Rels += Found_Rels
@@ -660,6 +661,8 @@ def _MPQS(N, verbose=True, M = 1):
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
 
+    merged_count = manager.Value("i", 0)
+
     while True:
         # trim primes, recalc min
         Prime_base = [p for p in Prime_base if p > min_prime]
@@ -675,7 +678,7 @@ def _MPQS(N, verbose=True, M = 1):
         for poly in polys:
             s1 = min(poly.start_vals[0]) 
             s2 = max(poly.start_vals[0]) 
-            inputs += [(Nm, start + s1, stop + s1 , Prime_base, smooth_base, Rels, required_relations,  poly)]
+            inputs += [(Nm, start + s1, stop + s1 , Prime_base, smooth_base, Rels, merged_count, required_relations,  poly)]
 
         # deploy tasks to every cpu core.
         pols = []
@@ -733,5 +736,6 @@ if __name__ == "__main__":
     ts = time.time()
     r = MPQS(N)
     print(r)
-    td = time.time()
-    sys.stderr.write("All done in: %f secs.\n" % (td-ts))
+    td = time.time() - ts
+    tds = humanfriendly.format_timespan(td) 
+    sys.stderr.write("All done in: %s.\n" % (tds))

From 9ec8253886c8c32023912882ec94705893f7f2be Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 11:42:05 -0300
Subject: [PATCH 16/56] lock interprocess counter, better estimation of ETA

---
 MPQS.sage | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 7ac9cbc..66d9626 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -450,17 +450,20 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
                     a = partials[f[1]]
                     p = filter_out_even_powers(f[0] + a[0])
                     Rels.append([p, y * a[1], Rad * a[2], A* a[3]])
-                    merged_count += 1
+                    with merged_count.get_lock():
+                        merged_count += 1
+                    del partials[f[1]]
                 else:
                     partials[f[1]] = [f[0], y, Rad, A]
         if i % m == 0:
+            lRels = len(Rels)
             lt = time.time()
             td = lt - ltd
             ltd = lt
-            eta = td * (ld/m)
+            eta = td * (((ld / m) + (required_relations / lRels)) / 2)
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
-            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged_count,tds,etas)
+            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),lRels,required_relations, merged_count,tds,etas)
             sys.stderr.write(msg)
     
     td = time.time() - st

From c814bd63aa775a53de69c339eb6d169a2be0e2ae Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 11:44:48 -0300
Subject: [PATCH 17/56] fix potential div by 0 bug

---
 MPQS.sage | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index 66d9626..a20f6b4 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -460,7 +460,10 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
             lt = time.time()
             td = lt - ltd
             ltd = lt
-            eta = td * (((ld / m) + (required_relations / lRels)) / 2)
+            if lRels > 0:
+                eta = td * (((ld / m) + (required_relations / lRels)) / 2)
+            else:
+                eta = td * (ld / m)
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
             msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),lRels,required_relations, merged_count,tds,etas)

From 19a0aa42041c9b6b3bdd8f347c83e5c1db2fd8bd Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 12:15:01 -0300
Subject: [PATCH 18/56] add new minifactor strategy

---
 MPQS.sage | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 54 insertions(+), 3 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index a20f6b4..8e65c01 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -264,6 +264,43 @@ def minifactor2(x, P, D = 1):
         return (tmp, R, x), D                 
 
 
+def pollard_rho(n, seed=2, p=2, c=1):
+        f = lambda x: x ** p + c
+        x, y, d = seed, seed, 1
+        while d == 1:
+            x = f(x) % n
+            y = f(f(y)) % n
+            d = gcd((x - y) % n, n)
+            if N > d > 1:
+                return d
+                
+def pollard_rho_iter(n, P):
+    tmp = []
+    r = [n]
+    R = n
+    while True:
+        x = r.pop()
+        if x < P[-1]
+            for p in P:
+                if R % p == 0:
+                    while x % p == 0:
+                        R // = p            
+                        pw += 1
+                    tmp.append((p,pw))
+        else:
+            p = pollard_rho_iter(x)
+            if is_prime(p):
+                if p in P:
+                    if R % p == 0:
+                        while R > 1 and R % p == 0:
+                            R // = p
+                            pw += 1
+                        tmp.append((p,pw))
+            #else:   
+            else:
+                r.append(p)
+    return merge_powers(tmp), R, n        
+                
 def minifactor3(x, P, smooth_base):
     """ 
     Minifactor. 
@@ -289,13 +326,27 @@ def minifactor4(x, P, smooth_base):
     smooth = gcd(x, smooth_base)
     if smooth > 1:
         not_smooth = x // smooth
-        a1, r1, n1 = trial_division(smooth, P)
+        a1, r1, n1 = pollard_rho_iter(smooth, P)
         a2, r2, n2 = trial_division(not_smooth, P)
     else:
         a2, r2, n2 = trial_division(x, P) 
     return merge_powers(a1 + a2), r2, x
 
-        
+
+def minifactor5(x, P, smooth_base):
+    """ 
+    Minifactor. 
+    """
+    x = abs(x)
+    a1 = []
+    smooth = gcd(x, smooth_base)
+    if smooth > 1:
+        not_smooth = x // smooth
+        a1, r1, n1 = trial_division(smooth, P)
+        a2, r2, n2 = pollard_rho_iter(not_smooth, P)
+    else:
+        a2, r2, n2 = pollard_rho_iter(x, P) 
+    return merge_powers(a1 + a2), r2, x
                          
 def minifactor(x, P):
     """ 
@@ -440,7 +491,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
             break
         yRad, x = Diffs[i]
         y, Rad = yRad
-        f = minifactor4(y, P, smooth_base)
+        f = minifactor5(y, P, smooth_base)
         if 1:
             if f != None:
                 filtered = (filter_out_even_powers(f[0]),f[1],y)

From cd35e048f848d626097f1c3b5f1a9d3d935aabe3 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 13:05:25 -0300
Subject: [PATCH 19/56] add new bounds estimator borrowed from msieve/mpqs.c

---
 MPQS.sage | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 53 insertions(+), 2 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 8e65c01..da3fdf2 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -664,6 +664,56 @@ def generate_polys(N, Prime_base, x_max, needed):
     return polys, early_factors
 
 
+def prebuilt_params(bits):
+    """
+    Bounds estimation
+    borrowed from msieve/mpqs.c
+    """
+    if bits <= 64:
+	    return [100, 40, 1 * 65536]
+	if bits <= 128: 
+	    return [450, 40, 1 * 65536]
+	if bits <= 183:
+	    return [2000, 40, 1 * 65536]
+	if bits <= 200: 
+	    return [3000, 50, 1 * 65536]
+	if bits <= 212: 
+	    return [5400, 50, 3 * 65536]
+	if bits <= 233
+	    return [10000, 100, 3 * 65536]
+	if bits <= 249:
+	    return [27000, 100, 3 * 65536]
+	if bits <= 266:
+	    return [50000, 100, 3 * 65536]
+	if bits <= 283:
+	    return [55000, 80, 3 * 65536]
+	if bits <= 298:
+	    return [60000, 80, 9 * 65536]
+	if bits <= 315:
+	    return [80000, 150, 9 * 65536]
+	if bits <= 332:
+	    return [100000, 150, 9 * 65536]
+	if bits <= 348:
+	    return [140000, 150, 9 * 65536]
+	if bits <= 363:
+	    return [210000, 150, 13 * 65536]
+	if bits <= 379:
+	    return [300000, 150, 17 * 65536]
+	if bits <= 395:
+	    return [400000, 150, 21 * 65536]
+	if bits <= 415:
+	    return [500000, 150, 25 * 65536] # beyond this point you're crazy 
+	if bits <= 440:
+	    return [700000, 150, 33 * 65536]
+	if bits <= 465:
+	    return [900000, 150, 50 * 65536]
+	if bits <= 490:
+	    return [1100000, 150, 75 * 65536]
+	if bits <= 512:
+	    return [1300000, 150, 100 * 65536]
+	return [1300000, 150, 100 * 65536]
+
+
 def _MPQS(N, verbose=True, M = 1):
     """ 
     Main MPQS function. 
@@ -676,8 +726,9 @@ def _MPQS(N, verbose=True, M = 1):
 
     T = cpu_count()
 
-    B1 = pow(int(log10(pow(N, 6))),2) 
-    Prime_base, log_p = find_primebase(N, B1)
+    #B1 = pow(int(log10(pow(N, 6))),2) 
+    B2, _ , B1 = prebuilt_params(log2(N))
+    Prime_base, log_p = find_primebase(N, B2)
     B1 //= M
     B2 = len(Prime_base)
   

From bfb1aae615692249bc30f5b27085d1ad3bd97496 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 13:22:24 -0300
Subject: [PATCH 20/56] forgot various:, fix lines 287, 301, fix tabs in func
 prebuilt_params

---
 MPQS.sage | 92 +++++++++++++++++++++++++++----------------------------
 1 file changed, 46 insertions(+), 46 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index da3fdf2..07a9b17 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -280,20 +280,20 @@ def pollard_rho_iter(n, P):
     R = n
     while True:
         x = r.pop()
-        if x < P[-1]
+        if x < P[-1]:
             for p in P:
                 if R % p == 0:
                     while x % p == 0:
-                        R // = p            
+                        R //= p            
                         pw += 1
                     tmp.append((p,pw))
         else:
-            p = pollard_rho_iter(x)
+            p = pollard_rho_iter(x, P)
             if is_prime(p):
                 if p in P:
                     if R % p == 0:
                         while R > 1 and R % p == 0:
-                            R // = p
+                            R //= p
                             pw += 1
                         tmp.append((p,pw))
             #else:   
@@ -670,48 +670,48 @@ def prebuilt_params(bits):
     borrowed from msieve/mpqs.c
     """
     if bits <= 64:
-	    return [100, 40, 1 * 65536]
-	if bits <= 128: 
-	    return [450, 40, 1 * 65536]
-	if bits <= 183:
-	    return [2000, 40, 1 * 65536]
-	if bits <= 200: 
-	    return [3000, 50, 1 * 65536]
-	if bits <= 212: 
-	    return [5400, 50, 3 * 65536]
-	if bits <= 233
-	    return [10000, 100, 3 * 65536]
-	if bits <= 249:
-	    return [27000, 100, 3 * 65536]
-	if bits <= 266:
-	    return [50000, 100, 3 * 65536]
-	if bits <= 283:
-	    return [55000, 80, 3 * 65536]
-	if bits <= 298:
-	    return [60000, 80, 9 * 65536]
-	if bits <= 315:
-	    return [80000, 150, 9 * 65536]
-	if bits <= 332:
-	    return [100000, 150, 9 * 65536]
-	if bits <= 348:
-	    return [140000, 150, 9 * 65536]
-	if bits <= 363:
-	    return [210000, 150, 13 * 65536]
-	if bits <= 379:
-	    return [300000, 150, 17 * 65536]
-	if bits <= 395:
-	    return [400000, 150, 21 * 65536]
-	if bits <= 415:
-	    return [500000, 150, 25 * 65536] # beyond this point you're crazy 
-	if bits <= 440:
-	    return [700000, 150, 33 * 65536]
-	if bits <= 465:
-	    return [900000, 150, 50 * 65536]
-	if bits <= 490:
-	    return [1100000, 150, 75 * 65536]
-	if bits <= 512:
-	    return [1300000, 150, 100 * 65536]
-	return [1300000, 150, 100 * 65536]
+        return [100, 40, 1 * 65536]
+    if bits <= 128: 
+        return [450, 40, 1 * 65536]
+    if bits <= 183:
+        return [2000, 40, 1 * 65536]
+    if bits <= 200: 
+        return [3000, 50, 1 * 65536]
+    if bits <= 212: 
+        return [5400, 50, 3 * 65536]
+    if bits <= 233:
+        return [10000, 100, 3 * 65536]
+    if bits <= 249:
+        return [27000, 100, 3 * 65536]
+    if bits <= 266:
+        return [50000, 100, 3 * 65536]
+    if bits <= 283:
+        return [55000, 80, 3 * 65536]
+    if bits <= 298:
+        return [60000, 80, 9 * 65536]
+    if bits <= 315:
+        return [80000, 150, 9 * 65536]
+    if bits <= 332:
+        return [100000, 150, 9 * 65536]
+    if bits <= 348:
+        return [140000, 150, 9 * 65536]
+    if bits <= 363:
+        return [210000, 150, 13 * 65536]
+    if bits <= 379:
+        return [300000, 150, 17 * 65536]
+    if bits <= 395:
+        return [400000, 150, 21 * 65536]
+    if bits <= 415:
+        return [500000, 150, 25 * 65536] # beyond this point you're crazy 
+    if bits <= 440:
+        return [700000, 150, 33 * 65536]
+    if bits <= 465:
+        return [900000, 150, 50 * 65536]
+    if bits <= 490:
+        return [1100000, 150, 75 * 65536]
+    if bits <= 512:
+        return [1300000, 150, 100 * 65536]
+    return [1300000, 150, 100 * 65536]
 
 
 def _MPQS(N, verbose=True, M = 1):

From 0fbafe6bbfd7cb20b8259a2c2b146485725a25d1 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 13:28:31 -0300
Subject: [PATCH 21/56] fix bug

---
 MPQS.sage | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 07a9b17..4ee7113 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -278,7 +278,7 @@ def pollard_rho_iter(n, P):
     tmp = []
     r = [n]
     R = n
-    while True:
+    while len(r) > 0:
         x = r.pop()
         if x < P[-1]:
             for p in P:
@@ -288,7 +288,7 @@ def pollard_rho_iter(n, P):
                         pw += 1
                     tmp.append((p,pw))
         else:
-            p = pollard_rho_iter(x, P)
+            p = pollard_rho(x)
             if is_prime(p):
                 if p in P:
                     if R % p == 0:

From 775e4b2de643980ee66d339790ba9f6b2fb7cae7 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 13:29:58 -0300
Subject: [PATCH 22/56] forgot ()

---
 MPQS.sage | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index 4ee7113..ef2e47a 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -479,7 +479,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
     A = pol.A
     #B = pol.B
     #C = pol.C
-    st = time.time
+    st = time.time()
     ltd = st
     ld = len(Diffs)
     m = 1000

From 8099666c2fc64b60b6edbe9c09f92f8715d93595 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 13:34:31 -0300
Subject: [PATCH 23/56] fix merge_count

---
 MPQS.sage | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index ef2e47a..601e4cf 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -502,7 +502,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
                     p = filter_out_even_powers(f[0] + a[0])
                     Rels.append([p, y * a[1], Rad * a[2], A* a[3]])
                     with merged_count.get_lock():
-                        merged_count += 1
+                        merged_count.value += 1
                     del partials[f[1]]
                 else:
                     partials[f[1]] = [f[0], y, Rad, A]
@@ -517,12 +517,12 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
                 eta = td * (ld / m)
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
-            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),lRels,required_relations, merged_count,tds,etas)
+            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),lRels,required_relations, merged_count.value,tds,etas)
             sys.stderr.write(msg)
     
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
-    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, time elapsed: %s\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged_count, tds)
+    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, time elapsed: %s\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged_count.value, tds)
     sys.stderr.write(msg)
     
     #Rels += Found_Rels

From f58f321c077deff263dd818d7ddb7296a9029eea Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 14:16:33 -0300
Subject: [PATCH 24/56] fix find_primebase

---
 MPQS.sage | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index 601e4cf..bae4449 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -609,7 +609,7 @@ def find_primebase(n, bound):
     Same as the sieve of erathostenes.
     """
     primes, mod_root, log_p, num_prime, p = [], [], [], 0, 3
-    while p < bound or num_prime < 3:
+    while num_prime <= bound:
         leg = legendre(n % p, p)
         if leg == 1:
             primes += [p]

From f8a8fc58964abfc9cd4892e72eeecf036e5d2a9f Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 14:22:21 -0300
Subject: [PATCH 25/56] use euclidean distance to estimate ETA

---
 MPQS.sage | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/MPQS.sage b/MPQS.sage
index bae4449..301fc73 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -512,7 +512,8 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
             td = lt - ltd
             ltd = lt
             if lRels > 0:
-                eta = td * (((ld / m) + (required_relations / lRels)) / 2)
+                #eta = td * (((ld / m) + (required_relations / lRels)) / 2)
+                eta = td * isqrt((ld / m)**2 + (required_relations / lRels)**2)
             else:
                 eta = td * (ld / m)
             tds = humanfriendly.format_timespan(td)

From 9664b99e544cf2287ecc1ef1af729e08ad65fdf3 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 18:04:55 -0300
Subject: [PATCH 26/56] Update MPQS.sage

---
 MPQS.sage | 55 ++++++++++++++++++++++++++-----------------------------
 1 file changed, 26 insertions(+), 29 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 301fc73..3f6e613 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -15,7 +15,7 @@ https://en.wikipedia.org/wiki/Quadratic_sieve
 
 import time
 import sys
-from gmpy2 import gcd, gcdext, isqrt, is_prime, next_prime, log2, log10, legendre, powmod, invert
+from gmpy2 import gcd, gcdext, sqrt, isqrt, is_prime, next_prime, log2, log10, legendre, powmod, invert
 from sage.parallel.multiprocessing_sage import parallel_iter
 from multiprocessing import cpu_count, Pool, Manager
 from itertools import repeat
@@ -292,6 +292,7 @@ def pollard_rho_iter(n, P):
             if is_prime(p):
                 if p in P:
                     if R % p == 0:
+                        pw = 0
                         while R > 1 and R % p == 0:
                             R //= p
                             pw += 1
@@ -326,7 +327,7 @@ def minifactor4(x, P, smooth_base):
     smooth = gcd(x, smooth_base)
     if smooth > 1:
         not_smooth = x // smooth
-        a1, r1, n1 = pollard_rho_iter(smooth, P)
+        a1, r1, n1 = trial_division(smooth, P)
         a2, r2, n2 = trial_division(not_smooth, P)
     else:
         a2, r2, n2 = trial_division(x, P) 
@@ -427,7 +428,7 @@ class Poly:
         B = self.B
         C = self.C
         n = self.n
-        start_vals = []
+        start_vals = {}
         #for p in self.P:
         p = self.root_A
         g = 1
@@ -442,8 +443,8 @@ class Poly:
                 r2 = (-1 * r1) % p
                 start1 = (ainv * (r1 - B)) % p
                 start2 = (ainv * (r2 - B)) % p
-                start_vals.append([int(start1), int(start2)])
-                self.start_vals = start_vals
+                start_vals[p] = [int(start1), int(start2)]
+        self.start_vals = start_vals
         return start_vals
 
 
@@ -464,17 +465,8 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
     Relations search funcion 
     """
     sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
-
-    #if (stop-start) < 0:
-    #    return [] 
-
-    if (stop-start) > 0:
-        I = [abs(x) for x in range(start, stop) if not is_prime(x) and not is_square(x)]
-    else: # this may not happen
-        I = [abs(x) for x in range(start, stop, -1) if not is_prime(x) and not is_square(x)]
-
     #D = reduce(lambda x, y: x * y, I)
-    Diffs = [(pol.eval(x),x) for x in I]
+    Diffs = [(pol.eval(abs(x)),abs(x)) for x in range(start, stop)]
     Found_Rels = []
     A = pol.A
     #B = pol.B
@@ -491,19 +483,22 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
             break
         yRad, x = Diffs[i]
         y, Rad = yRad
-        f = minifactor5(y, P, smooth_base)
+        y = abs(y)
+        f = minifactor4(y, P, smooth_base)
         if 1:
             if f != None:
-                filtered = (filter_out_even_powers(f[0]),f[1],y)
-                if f[1] == 1:    
-                    Rels.append([filtered, y, Rad, A])
+                filtered = filter_out_even_powers(f[0])
+                if f[1] == 1:   
+                    rel = [filtered, y, Rad, A] 
+                    #print(rel)
+                    Rels.append(rel)
                 elif f[1] in partials:
                     a = partials[f[1]]
                     p = filter_out_even_powers(f[0] + a[0])
                     Rels.append([p, y * a[1], Rad * a[2], A* a[3]])
-                    with merged_count.get_lock():
-                        merged_count.value += 1
-                    del partials[f[1]]
+                    #with merged_count.value.get_lock():
+                    merged_count.value += 1
+                    #del partials[f[1]]
                 else:
                     partials[f[1]] = [f[0], y, Rad, A]
         if i % m == 0:
@@ -513,9 +508,9 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
             ltd = lt
             if lRels > 0:
                 #eta = td * (((ld / m) + (required_relations / lRels)) / 2)
-                eta = td * isqrt((ld / m)**2 + (required_relations / lRels)**2)
+                eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
             else:
-                eta = td * (ld / m)
+                eta = int(td * (ld / m))
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
             msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),lRels,required_relations, merged_count.value,tds,etas)
@@ -528,6 +523,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
     
     #Rels += Found_Rels
     return Found_Rels
+    
 
 def transpose(A):
     """
@@ -583,6 +579,8 @@ def process_basis_vectors(N, basis, Rels, multiplier = 1):
     for K in basis:
         lhs = rhs = Ahs = 1
         I = [f for f, k in zip(Rels, K) if k == 1]
+        #print(I)
+        #sys.exit(0)
         if len(I) > 0:
             for i in I:
                 lhs *= i[1] # left-hand side
@@ -609,12 +607,12 @@ def find_primebase(n, bound):
     https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
     Same as the sieve of erathostenes.
     """
-    primes, mod_root, log_p, num_prime, p = [], [], [], 0, 3
+    primes, mod_root, log_p, num_prime, p = [], {}, [], 0, 3
     while num_prime <= bound:
         leg = legendre(n % p, p)
         if leg == 1:
             primes += [p]
-            #mod_root += mod_sqrt(n ,p)
+            #mod_root[p] = mod_sqrt(n ,p)
             log_p += [log10(p)]
             num_prime += 1
         elif leg == 0:
@@ -785,9 +783,8 @@ def _MPQS(N, verbose=True, M = 1):
         
         # generate tasks parameters
         for poly in polys:
-            s1 = min(poly.start_vals[0]) 
-            s2 = max(poly.start_vals[0]) 
-            inputs += [(Nm, start + s1, stop + s1 , Prime_base, smooth_base, Rels, merged_count, required_relations,  poly)]
+            inputs += [(Nm, start, stop, Prime_base, smooth_base, Rels, merged_count, required_relations, poly)]
+
 
         # deploy tasks to every cpu core.
         pols = []

From 44bb8989d393377965eb3046da347340535d01b1 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Wed, 9 Jun 2021 19:12:57 -0300
Subject: [PATCH 27/56] Update MPQS.sage

---
 MPQS.sage | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 3f6e613..b83f09d 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -428,12 +428,12 @@ class Poly:
         B = self.B
         C = self.C
         n = self.n
-        start_vals = {}
+        start_vals = []
         #for p in self.P:
         p = self.root_A
         g = 1
-        while g == 1:
         #for p in self.P:
+        while g == 1:
             p = next_prime(p)
             ainv = 1
             if A != 1:
@@ -443,7 +443,7 @@ class Poly:
                 r2 = (-1 * r1) % p
                 start1 = (ainv * (r1 - B)) % p
                 start2 = (ainv * (r2 - B)) % p
-                start_vals[p] = [int(start1), int(start2)]
+                start_vals = [[int(start1), int(start2)],p]
         self.start_vals = start_vals
         return start_vals
 
@@ -464,9 +464,11 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
     """ 
     Relations search funcion 
     """
-    sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
+    s = min(pol.start_vals[0])
+    #s = 0
+    sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start + s, stop + s, (stop-start)))
     #D = reduce(lambda x, y: x * y, I)
-    Diffs = [(pol.eval(abs(x)),abs(x)) for x in range(start, stop)]
+    Diffs = [(pol.eval(abs(x)),abs(x)) for x in range(start + s, stop + s)]
     Found_Rels = []
     A = pol.A
     #B = pol.B
@@ -486,9 +488,9 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
         y = abs(y)
         f = minifactor4(y, P, smooth_base)
         if 1:
-            if f != None:
-                filtered = filter_out_even_powers(f[0])
+            if f != None:            
                 if f[1] == 1:   
+                    filtered = filter_out_even_powers(f[0])
                     rel = [filtered, y, Rad, A] 
                     #print(rel)
                     Rels.append(rel)

From 54c4aff1f5c2104ab09f7d0d32e3e42034776ad5 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Thu, 10 Jun 2021 18:54:51 -0300
Subject: [PATCH 28/56] Update MPQS.sage

---
 MPQS.sage | 141 ++++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 100 insertions(+), 41 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index b83f09d..c57570e 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -363,9 +363,10 @@ class Poly:
     Quadratic polynomial helper class: 
     type Ax + 2Bx + C 
     """
-    def __init__(self, n, P, x_max, search = 0, verbose = None):
+    def __init__(self, n, P, x_max, search = 0, verbose = None, logp = None):
         self.n = int(n)
         self.P = P
+        self.logp = logp
         self.x_max = x_max
         self.search = search
         self.verbose = verbose
@@ -445,8 +446,30 @@ class Poly:
                 start2 = (ainv * (r2 - B)) % p
                 start_vals = [[int(start1), int(start2)],p]
         self.start_vals = start_vals
-        return start_vals
 
+    def calc_sums(self):    
+        x_max = self.x_max
+        a,b = self.start_vals[0]
+        sums = [0] * x_max * 2
+        logp = {}
+        for p in self.P:
+            if self.logp == None:
+                logp[p] = log10(p)
+            else:                 
+                logp[p] = self.logp[p]
+            k = 0
+            while k < x_max:
+                if k+a < x_max:
+                    sums[k+a] += logp[p]
+                if k+b < x_max:
+                    sums[k+b] += logp[p]
+                if k:
+                    sums[k-a+x_max] += logp[p]
+                    sums[k-b+x_max] += logp[p]
+                k += p
+        self.sums = sums
+        #print(sums)
+        return sums
 
     def eval(self, x):
         """ 
@@ -460,17 +483,18 @@ class Poly:
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
         
 
-def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_relations, pol = None):
+def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, poly = None):
     """ 
     Relations search funcion 
     """
-    s = min(pol.start_vals[0])
-    #s = 0
-    sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start + s, stop + s, (stop-start)))
+    sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
     #D = reduce(lambda x, y: x * y, I)
-    Diffs = [(pol.eval(abs(x)),abs(x)) for x in range(start + s, stop + s)]
+   
+    #sums = poly.calc_sums()
+    Diffs = [(poly.eval(abs(x)),abs(x)) for x in range(start, stop)]
+
     Found_Rels = []
-    A = pol.A
+    A = poly.A
     #B = pol.B
     #C = pol.C
     st = time.time()
@@ -479,15 +503,19 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
     m = 1000
     msg = ""
     partials = {}
-
+    i= 0
     for i in range(ld):
-        if len(Rels) > required_relations:
+        if len(Rels) > required_relations or len(cycleFactors) > 0:
             break
         yRad, x = Diffs[i]
         y, Rad = yRad
         y = abs(y)
-        f = minifactor4(y, P, smooth_base)
-        if 1:
+        l10 = log10(y) 
+        if l10 < thresh:
+            print("Value %d: %d < %d" % (y,l10,thresh))
+        else:
+            f = minifactor4(y, P, smooth_base)
+        #if 1:
             if f != None:            
                 if f[1] == 1:   
                     filtered = filter_out_even_powers(f[0])
@@ -497,7 +525,16 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
                 elif f[1] in partials:
                     a = partials[f[1]]
                     p = filter_out_even_powers(f[0] + a[0])
-                    Rels.append([p, y * a[1], Rad * a[2], A* a[3]])
+                    Ahs = A*a[3]
+                    lhs = Rad * a[2]
+                    rhs = y * a[1]
+                    LHS = Ahs + lhs
+                    g = gcd(isqrt(LHS) - rhs, N)
+                    if N > g > 1:
+                        cycleFactors.append([g , N // g])
+                        sys.stderr.write("Found cycle with partial\n")
+                    else:
+                        Rels.append([p, rhs, lhs, Ahs])
                     #with merged_count.value.get_lock():
                     merged_count.value += 1
                     #del partials[f[1]]
@@ -515,12 +552,12 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
                 eta = int(td * (ld / m))
             tds = humanfriendly.format_timespan(td)
             etas = humanfriendly.format_timespan(eta)
-            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,(stop-start),lRels,required_relations, merged_count.value,tds,etas)
+            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,ld,lRels,required_relations, merged_count.value,tds,etas)
             sys.stderr.write(msg)
     
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
-    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, time elapsed: %s\n" % (start,stop,i,(stop-start),len(Rels),required_relations, merged_count.value, tds)
+    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, time elapsed: %s\n" % (start,stop,i,ld,len(Rels),required_relations, merged_count.value, tds)
     sys.stderr.write(msg)
     
     #Rels += Found_Rels
@@ -609,13 +646,13 @@ def find_primebase(n, bound):
     https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
     Same as the sieve of erathostenes.
     """
-    primes, mod_root, log_p, num_prime, p = [], {}, [], 0, 3
+    primes, mod_root, log_p, num_prime, p = [], {}, {}, 0, 3
     while num_prime <= bound:
         leg = legendre(n % p, p)
         if leg == 1:
             primes += [p]
             #mod_root[p] = mod_sqrt(n ,p)
-            log_p += [log10(p)]
+            log_p[p] = log10(p)
             num_prime += 1
         elif leg == 0:
             return [p], None
@@ -632,19 +669,19 @@ def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
     min_prime = int(thresh * 3)
     fudge = sum(log_p[i] for i, p in enumerate(Prime_base) if p < min_prime)
     fudge = fudge // 4
-    sys.stderr.write("min_prime: %d, fudge: %f\n" % (min_prime,fudge))
+    sys.stderr.write("min_prime: %d, fudge: %f, thresh: %f\n" % (min_prime,fudge, thresh))
     thresh -= fudge
     return min_prime, thresh, fudge
 
 
-def generate_polys(N, Prime_base, x_max, needed):
+def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
     """ 
     It searchs for distinct needed polys congruent to n. 
     """
-    n=1
+    n = min_search
     cpolys = 0
     polys = []
-    polys_ABC = []
+    #polys_ABC = []
     early_factors = []
     while cpolys <= needed:
         pol = Poly(N, Prime_base, x_max, search = n, verbose=False)
@@ -662,7 +699,7 @@ def generate_polys(N, Prime_base, x_max, needed):
                 polys.append(pol)
                 cpolys += 1
         n += 1
-    return polys, early_factors
+    return polys, polys_ABC, early_factors
 
 
 def prebuilt_params(bits):
@@ -740,7 +777,8 @@ def _MPQS(N, verbose=True, M = 1):
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
     
-    x_max = B2 *60  # size of the sieve
+    #x_max = B2 *60  # size of the sieve
+    x_max = B1
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
     min_prime = int(thresh * 3)
@@ -755,23 +793,27 @@ def _MPQS(N, verbose=True, M = 1):
     start = 0
     stop = B1 # range to sieve
 
-    polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
-    if len(early_factors) > 0:
-        tmp = 1
-        small = []
-        for early_factor in early_factors:
-            g = gcd(early_factor, N)    
-            if N > g > 1: 
-                tmp *= g
-                small.append(g)
-        if tmp > 1:
-            return small + MPQS(N // tmp)
+    #polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
+    #if len(early_factors) > 0:
+    #    tmp = 1
+    #    small = []
+    #    for early_factor in early_factors:
+    #        g = gcd(early_factor, N)    
+    #        if N > g > 1: 
+    #            tmp *= g
+    #            small.append(g)
+    #    if tmp > 1:
+    #        return small + MPQS(N // tmp)
     
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
+    cycleFactors = manager.list()
 
     merged_count = manager.Value("i", 0)
 
+    min_search = 0
+    polys_ABC = []
+
     while True:
         # trim primes, recalc min
         Prime_base = [p for p in Prime_base if p > min_prime]
@@ -780,13 +822,26 @@ def _MPQS(N, verbose=True, M = 1):
 
         t1 = time.time()
         sys.stderr.write("Data collection with %d threads...\n" % T)
-
-        inputs = []
-        
+  
+        polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T, min_search, polys_ABC) # generate n distinct polys one for each cpu core.
+        if len(early_factors) > 0:
+            tmp = 1
+            small = []
+            for early_factor in early_factors:
+                g = gcd(early_factor, N)
+                if N > g > 1:
+                    tmp *= g
+                    small.append(g)
+            if tmp > 1:
+                return small + MPQS(N // tmp)
+
+        inputs = [] 
         # generate tasks parameters
         for poly in polys:
-            inputs += [(Nm, start, stop, Prime_base, smooth_base, Rels, merged_count, required_relations, poly)]
-
+            #print(poly.sums)
+            #calc_sums
+            inputs += [(Nm, start, stop, Prime_base, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, poly)]
+        #sys.exit(0)
 
         # deploy tasks to every cpu core.
         pols = []
@@ -800,6 +855,9 @@ def _MPQS(N, verbose=True, M = 1):
         sys.stderr.write("Done in: %f secs.\n" % (t2-t1))
         sys.stderr.write("Found %d rels with %d base primes.\n" % (len(Rels),len(Prime_base)))
 
+        if len(cycleFactors) > 0:
+            return cycleFactors
+
         # when needed relations is reached proceed to do linear algebra
         if len(Rels) > required_relations:
             sys.stderr.write("Found %d enough relations of %d needed relations...\n" % (len(Rels),required_relations))
@@ -817,8 +875,9 @@ def _MPQS(N, verbose=True, M = 1):
             if result != None:
                 return result
 
-        start = stop
-        stop *= 2
+        #start = stop
+        #stop += B1
+        min_search += T
 
         sys.stderr.write("Need to sieve %d differences more...\nNew sieving range %d:%d\n" % (required_relations - len(Rels),start,stop))
       

From 265ae7a184bdaad0cbf3ff3a40747f3060b8ae50 Mon Sep 17 00:00:00 2001
From: Dario Clavijo <daedalus2027@gmail.com>
Date: Fri, 11 Jun 2021 19:42:41 -0300
Subject: [PATCH 29/56] cleanup code, add test for perfect powers, primes,
 squares before minifactor

---
 MPQS.sage | 290 +++++++++++++++++++-----------------------------------
 1 file changed, 99 insertions(+), 191 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index c57570e..ab9d0a3 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -13,9 +13,10 @@ It was invented by Carl Pomerance in 1981 as an improvement to Schroeppel's line
 https://en.wikipedia.org/wiki/Quadratic_sieve
 """
 
+import os
 import time
 import sys
-from gmpy2 import gcd, gcdext, sqrt, isqrt, is_prime, next_prime, log2, log10, legendre, powmod, invert
+from gmpy2 import gcd, gcdext, sqrt, isqrt, is_prime, next_prime, log2, log10, log, legendre, powmod, invert
 from sage.parallel.multiprocessing_sage import parallel_iter
 from multiprocessing import cpu_count, Pool, Manager
 from itertools import repeat
@@ -105,8 +106,6 @@ def choose_multiplier(n, prime_list):
     best_score = 1000.0
     best_mult = 1
 
-    #print(scores)
-   
     for i in range(0, num_multipliers):
         score = scores[i];
         if (score < best_score):
@@ -150,6 +149,16 @@ def mod_sqrt(n, p):
         return a  # p == 2
 
 
+def is_power_logprime(n, min_log_primes):
+    ispow = False
+    for p in min_log_primes:
+        a = log(n) / p
+        b = int(a)
+        if a == b:
+            iwpow = True
+    return ispow
+
+
 def trial_division(n, P):
     """ 
     A simple trial division, factors are given by P-list. 
@@ -158,16 +167,20 @@ def trial_division(n, P):
     r = n
     l = len(P)
     i = 0
-    pi2 = pow(P[i],2)
-    while (i < l) and (pi2 <= n):
-        if r % P[i] == 0:
-            pw = 0
-            while r % P[i] == 0:
-                pw += 1
-                r //= P[i]
-            a.append((int(P[i]),int(pw)))
-        else:
-            i += 1
+    if n > P[0]:
+        while (i < l):
+            if r >= P[i]:
+            #if 1:
+                if r % P[i] == 0:
+                    pw = 0
+                    while r % P[i] == 0:
+                        pw += 1
+                        r //= P[i]
+                    a.append((int(P[i]),int(pw)))
+                else:
+                    i += 1
+            else:
+                break
     return a,r,n
 
 def merge_powers(ppws):
@@ -219,106 +232,7 @@ def is_smooth(x, P):
     return abs(y) == 1 
 
 
-def minifactor2(x, P, D = 1):
-    """ 
-    Tries to factor out common primes with gcd and filter out even power primes. 
-    """
-    tmp = []
-    R = 1
-    g = gcd(x, D) 
-    if x > g > 1: # if there is a common factor try to factor primes from that factor
-        gtmp = [g, x//g]
-        for g1 in gtmp:
-            if not is_prime(g1):
-                i = isqrt(abs(g1))
-                if i**2 != abs(g1): # g1 should not be square
-                    a, r, n = trial_division(g1,P)
-                    tmp += a
-                    R *= r
-                else:
-                    if is_prime(i):
-                        if i > P[-1] or i < P[0]:  
-                            R *= i
-                        else:
-                        #a, r, n = [i], 1, g1
-                            tmp += [(i,2)]
-                            #R //= g1
-                    else:
-                        #R *= r
-                        gtmp.append(i)    
-            else:
-                if g1 > P[-1] or g1 < P[0]:
-                    R *= g1
-                else:
-                    tmp += [(g1,1)]
-                    #R //= g1
-        D//=g # reduce D
-    else: # if not proceed as always
-        a, r, n = trial_division(x,P)
-        R *= r
-        tmp += a
-        D//=x # reduce D
-    #print(tmp)
-    if R == 1: # if remainder is 1 then we found a B-smooth number then relation.
-        tmp = filter_out_even_powers(a)  
-        return (tmp, R, x), D                 
-
-
-def pollard_rho(n, seed=2, p=2, c=1):
-        f = lambda x: x ** p + c
-        x, y, d = seed, seed, 1
-        while d == 1:
-            x = f(x) % n
-            y = f(f(y)) % n
-            d = gcd((x - y) % n, n)
-            if N > d > 1:
-                return d
-                
-def pollard_rho_iter(n, P):
-    tmp = []
-    r = [n]
-    R = n
-    while len(r) > 0:
-        x = r.pop()
-        if x < P[-1]:
-            for p in P:
-                if R % p == 0:
-                    while x % p == 0:
-                        R //= p            
-                        pw += 1
-                    tmp.append((p,pw))
-        else:
-            p = pollard_rho(x)
-            if is_prime(p):
-                if p in P:
-                    if R % p == 0:
-                        pw = 0
-                        while R > 1 and R % p == 0:
-                            R //= p
-                            pw += 1
-                        tmp.append((p,pw))
-            #else:   
-            else:
-                r.append(p)
-    return merge_powers(tmp), R, n        
-                
-def minifactor3(x, P, smooth_base):
-    """ 
-    Minifactor. 
-    """
-    x = abs(x)
-    smooth = gcd(x, smooth_base)
-    if smooth > 1:
-        not_smooth = x // smooth
-        a1, r1, n1 = trial_division(smooth, P)
-        a2, r2, n2 = trial_division(not_smooth, P)
-        if r2 == 1:
-            a = a1 + a2
-            a3 = filter_out_even_powers(a)
-            return a3, r2, x
-
-
-def minifactor4(x, P, smooth_base):
+def minifactor(x, P, smooth_base):
     """ 
     Minifactor. 
     """
@@ -334,30 +248,6 @@ def minifactor4(x, P, smooth_base):
     return merge_powers(a1 + a2), r2, x
 
 
-def minifactor5(x, P, smooth_base):
-    """ 
-    Minifactor. 
-    """
-    x = abs(x)
-    a1 = []
-    smooth = gcd(x, smooth_base)
-    if smooth > 1:
-        not_smooth = x // smooth
-        a1, r1, n1 = trial_division(smooth, P)
-        a2, r2, n2 = pollard_rho_iter(not_smooth, P)
-    else:
-        a2, r2, n2 = pollard_rho_iter(x, P) 
-    return merge_powers(a1 + a2), r2, x
-                         
-def minifactor(x, P):
-    """ 
-    Minifactor algo, finds odd-power primes in composites. 
-    """
-    p = trial_division_minus_even_powers(x,P)
-    if p[1] == 1: # if 1 x is B-smooth
-        return p
-
-
 class Poly:
     """ 
     Quadratic polynomial helper class: 
@@ -388,7 +278,9 @@ class Poly:
         s=0
         leg = 2
         while True:
-            root_A = next_prime(root_A)
+            root_A = next_prime(root_A + s)
+            #root_A = next_prime(root_A)
+
             leg = legendre(n, root_A)
             if leg == 1:
                 if s > self.search:
@@ -483,20 +375,39 @@ class Poly:
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
         
 
-def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, poly = None):
+def is_power(n, minprimes):
+    ispow = False
+    c = 0
+    for p in minprimes:
+        a = log(n) / log(p)
+        b = int(a)
+        if a == b:
+            ispow = True
+            c = b
+    return ispow
+
+
+def is_power_logprime(n, min_log_primes):
+    ispow = False
+    for p in min_log_primes:
+        a = log(n) / p
+        b = int(a)
+        if a == b:
+            iwpow = True
+    return ispow
+
+def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, poly = None):
     """ 
     Relations search funcion 
     """
-    sys.stderr.write("relations_find: range(%d, %d), interval: %d sieving start\n" % (start, stop, (stop-start)))
-    #D = reduce(lambda x, y: x * y, I)
-   
-    #sums = poly.calc_sums()
+    #pid = os.getpid()
+    pid = taskid
+    tasks.value += 1
+    sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d sieving start\n" % (pid,start, stop, (stop-start)))
     Diffs = [(poly.eval(abs(x)),abs(x)) for x in range(start, stop)]
 
     Found_Rels = []
     A = poly.A
-    #B = pol.B
-    #C = pol.C
     st = time.time()
     ltd = st
     ld = len(Diffs)
@@ -510,17 +421,12 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
         yRad, x = Diffs[i]
         y, Rad = yRad
         y = abs(y)
-        l10 = log10(y) 
-        if l10 < thresh:
-            print("Value %d: %d < %d" % (y,l10,thresh))
-        else:
-            f = minifactor4(y, P, smooth_base)
-        #if 1:
+        if not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes):
+            f = minifactor(y, P, smooth_base)
             if f != None:            
                 if f[1] == 1:   
                     filtered = filter_out_even_powers(f[0])
                     rel = [filtered, y, Rad, A] 
-                    #print(rel)
                     Rels.append(rel)
                 elif f[1] in partials:
                     a = partials[f[1]]
@@ -535,9 +441,7 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
                         sys.stderr.write("Found cycle with partial\n")
                     else:
                         Rels.append([p, rhs, lhs, Ahs])
-                    #with merged_count.value.get_lock():
                     merged_count.value += 1
-                    #del partials[f[1]]
                 else:
                     partials[f[1]] = [f[0], y, Rad, A]
         if i % m == 0:
@@ -545,22 +449,23 @@ def relations_find(N, start, stop, P, smooth_base, Rels, merged_count, required_
             lt = time.time()
             td = lt - ltd
             ltd = lt
-            if lRels > 0:
-                #eta = td * (((ld / m) + (required_relations / lRels)) / 2)
-                eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
+            if i > 0:
+                if lRels > 0:
+                    eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
+                else:
+                    eta = int(td * (ld / m))
+                etas = humanfriendly.format_timespan(eta)
             else:
-                eta = int(td * (ld / m))
+                etas = "No ETA estimated"
             tds = humanfriendly.format_timespan(td)
-            etas = humanfriendly.format_timespan(eta)
-            msg = "relations_find: range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s.\n" % (start,stop,i,ld,lRels,required_relations, merged_count.value,tds,etas)
+            msg = "[%d] relations_find: range (%d, %d , %d), found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s, alive tasks: %d.\n" % (pid, start,i,ld,lRels,required_relations, merged_count.value,tds,etas,tasks.value)
             sys.stderr.write(msg)
     
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
-    msg = "relations_find: Ended range(%d, %d), inverval: %d of %d, found: %d of %d, merged: %d, time elapsed: %s\n" % (start,stop,i,ld,len(Rels),required_relations, merged_count.value, tds)
+    msg = "[%d] relations_find: Ended range(%d, %d, %d), found: %d of %d, merged: %d, time elapsed: %s, alive tasks: %d.\n" % (pid, start,i,ld,len(Rels),required_relations, merged_count.value, tds, tasks.value)
     sys.stderr.write(msg)
-    
-    #Rels += Found_Rels
+    tasks.value -= 1
     return Found_Rels
     
 
@@ -576,6 +481,7 @@ def transpose(A):
         new_A.append(new_row)
     return(new_A)
 
+
 def Gaussian_elimination_GF2(A):
   """
   Gaussian elimination over GF2
@@ -594,6 +500,7 @@ def Gaussian_elimination_GF2(A):
           break
   return marks, A
 
+
 def left_nullspace(A):
     """
     Compute left null space:
@@ -601,7 +508,13 @@ def left_nullspace(A):
     """
     A = transpose(A)
     marks, A = Gaussian_elimination_GF2(A)
-    return marks, A
+
+    B = []
+    for row in range(len(A)):
+        if marks[row] == False:
+            B.append(A[row])
+    return B
+
 
 def linear_algebra(Rels, P):
     """ 
@@ -618,8 +531,6 @@ def process_basis_vectors(N, basis, Rels, multiplier = 1):
     for K in basis:
         lhs = rhs = Ahs = 1
         I = [f for f, k in zip(Rels, K) if k == 1]
-        #print(I)
-        #sys.exit(0)
         if len(I) > 0:
             for i in I:
                 lhs *= i[1] # left-hand side
@@ -679,13 +590,11 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
     It searchs for distinct needed polys congruent to n. 
     """
     n = min_search
-    cpolys = 0
+    cpolys = 1
     polys = []
-    #polys_ABC = []
     early_factors = []
     while cpolys <= needed:
         pol = Poly(N, Prime_base, x_max, search = n, verbose=False)
-        #print(pol.early_factors)
         if len(pol.early_factors) > 0:
             for early_factor in pol.early_factors:
                 if early_factor not in early_factors:
@@ -756,15 +665,20 @@ def _MPQS(N, verbose=True, M = 1):
     """ 
     Main MPQS function. 
     """
+    if is_prime(N):
+        return [N]
+
     bN, lN = int(log2(N)), len(str(N))
-    i2N = isqrt(N) 
+    i2N = isqrt(N)
+
+    if pow(i2N,2) == N:
+        return [i2N, i2N]
+
     i2Np1 = i2N + 1 
     root_2n = isqrt(2*N)
     Rels = []
 
     T = cpu_count()
-
-    #B1 = pow(int(log10(pow(N, 6))),2) 
     B2, _ , B1 = prebuilt_params(log2(N))
     Prime_base, log_p = find_primebase(N, B2)
     B1 //= M
@@ -777,7 +691,6 @@ def _MPQS(N, verbose=True, M = 1):
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
     
-    #x_max = B2 *60  # size of the sieve
     x_max = B1
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
@@ -793,37 +706,28 @@ def _MPQS(N, verbose=True, M = 1):
     start = 0
     stop = B1 # range to sieve
 
-    #polys, early_factors = generate_polys(Nm, Prime_base, x_max, T) # generate n distinct polys one for each cpu core.
-    #if len(early_factors) > 0:
-    #    tmp = 1
-    #    small = []
-    #    for early_factor in early_factors:
-    #        g = gcd(early_factor, N)    
-    #        if N > g > 1: 
-    #            tmp *= g
-    #            small.append(g)
-    #    if tmp > 1:
-    #        return small + MPQS(N // tmp)
-    
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
     cycleFactors = manager.list()
 
     merged_count = manager.Value("i", 0)
+    tasks = manager.Value("i", 0)
 
     min_search = 0
     polys_ABC = []
 
     while True:
         # trim primes, recalc min
+        min_log_primes = [log(p) for p in Prime_base if p <= min_prime]
         Prime_base = [p for p in Prime_base if p > min_prime]
+        log_primes = [log(p) for p in Prime_base]
         smooth_base = prod(Prime_base)
         min_prime, thresh, fudge = recalculate_min_prime_thresh(thresh, Prime_base, log_p)
 
         t1 = time.time()
         sys.stderr.write("Data collection with %d threads...\n" % T)
   
-        polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T, min_search, polys_ABC) # generate n distinct polys one for each cpu core.
+        polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_search, polys_ABC) # generate n distinct polys one for each cpu core.
         if len(early_factors) > 0:
             tmp = 1
             small = []
@@ -836,12 +740,11 @@ def _MPQS(N, verbose=True, M = 1):
                 return small + MPQS(N // tmp)
 
         inputs = [] 
+        taskid = 1
         # generate tasks parameters
         for poly in polys:
-            #print(poly.sums)
-            #calc_sums
-            inputs += [(Nm, start, stop, Prime_base, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, poly)]
-        #sys.exit(0)
+            inputs += [(taskid, Nm, start, stop, Prime_base, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, poly)]
+            taskid += 1
 
         # deploy tasks to every cpu core.
         pols = []
@@ -902,7 +805,12 @@ if __name__ == "__main__":
     N = Integer(sys.argv[1])
     ts = time.time()
     r = MPQS(N)
-    print(r)
+    Fp = trial_division(N,sorted(r))[0] 
+    tmp = []
+    for f,p in Fp:
+        tmp += ["%d ^ %d" % (f,p)]
+    msg = "The factorization for N = %d is: %s" % (N," * ".join(tmp))
+    print(msg.replace(" ^ 1",""))
     td = time.time() - ts
     tds = humanfriendly.format_timespan(td) 
     sys.stderr.write("All done in: %s.\n" % (tds))

From ca6667a461dffc248226a94a2f78352b87e81bfc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sat, 12 Jun 2021 01:42:01 -0300
Subject: [PATCH 30/56] clean up code

---
 MPQS.sage | 168 ++++++++++++++++++++++++++++--------------------------
 1 file changed, 87 insertions(+), 81 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index ab9d0a3..4527842 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -106,6 +106,8 @@ def choose_multiplier(n, prime_list):
     best_score = 1000.0
     best_mult = 1
 
+    #print(scores)
+   
     for i in range(0, num_multipliers):
         score = scores[i];
         if (score < best_score):
@@ -167,22 +169,19 @@ def trial_division(n, P):
     r = n
     l = len(P)
     i = 0
-    if n > P[0]:
-        while (i < l):
-            if r >= P[i]:
-            #if 1:
-                if r % P[i] == 0:
-                    pw = 0
-                    while r % P[i] == 0:
-                        pw += 1
-                        r //= P[i]
-                    a.append((int(P[i]),int(pw)))
-                else:
-                    i += 1
-            else:
-                break
+    for p in P:
+        if r >= p:
+            if r % p == 0:
+                pw = 0
+                while r % p == 0:
+                    pw += 1
+                    r //= p
+               a.append((int(p),int(pw)))
+        else:
+            break
     return a,r,n
 
+
 def merge_powers(ppws):
     d = {}
     for p,pw in ppws:
@@ -195,6 +194,7 @@ def merge_powers(ppws):
        tmp2.append((p,d[p]))
     return tmp2
 
+
 def filter_out_even_powers(ppws):
     """ 
     Filter out even powers. 
@@ -232,7 +232,7 @@ def is_smooth(x, P):
     return abs(y) == 1 
 
 
-def minifactor(x, P, smooth_base):
+def minifactor4(x, P, smooth_base):
     """ 
     Minifactor. 
     """
@@ -315,8 +315,7 @@ class Poly:
         """ 
         Get the mininmum x value of the polynomial where y=f(x) is a guaranteed relation. 
         Code borrowed https://github.com/cramppet/quadratic-sieve/blob/master/QS.py
-        """
-        
+        """        
         A = self.A
         B = self.B
         C = self.C
@@ -330,38 +329,15 @@ class Poly:
             p = next_prime(p)
             ainv = 1
             if A != 1:
-                 g, inv, _ = gcdext(A, p)
-            if g != 1:
-                r1 = mod_sqrt(C, p)
-                r2 = (-1 * r1) % p
-                start1 = (ainv * (r1 - B)) % p
-                start2 = (ainv * (r2 - B)) % p
-                start_vals = [[int(start1), int(start2)],p]
-        self.start_vals = start_vals
-
-    def calc_sums(self):    
-        x_max = self.x_max
-        a,b = self.start_vals[0]
-        sums = [0] * x_max * 2
-        logp = {}
-        for p in self.P:
-            if self.logp == None:
-                logp[p] = log10(p)
-            else:                 
-                logp[p] = self.logp[p]
-            k = 0
-            while k < x_max:
-                if k+a < x_max:
-                    sums[k+a] += logp[p]
-                if k+b < x_max:
-                    sums[k+b] += logp[p]
-                if k:
-                    sums[k-a+x_max] += logp[p]
-                    sums[k-b+x_max] += logp[p]
-                k += p
-        self.sums = sums
-        #print(sums)
-        return sums
+                g, inv, _ = gcdext(A, p)
+                if g != 1:
+                    r1 = mod_sqrt(C, p)
+                    r2 = (-1 * r1) % p
+                    start1 = (ainv * (r1 - B)) % p
+                    start2 = (ainv * (r2 - B)) % p
+                    self.X = [int(start1), int(start2)]
+                    self.minima = sum(self.X)//2
+
 
     def eval(self, x):
         """ 
@@ -396,6 +372,7 @@ def is_power_logprime(n, min_log_primes):
             iwpow = True
     return ispow
 
+
 def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, poly = None):
     """ 
     Relations search funcion 
@@ -404,8 +381,9 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
     pid = taskid
     tasks.value += 1
     sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d sieving start\n" % (pid,start, stop, (stop-start)))
-    Diffs = [(poly.eval(abs(x)),abs(x)) for x in range(start, stop)]
-
+    m = poly.minima
+    Diffs = [(poly.eval(abs(x)),abs(x)) for x in range(start + m, stop +m)]
+    proc = noproc = 0
     Found_Rels = []
     A = poly.A
     st = time.time()
@@ -421,12 +399,14 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
         yRad, x = Diffs[i]
         y, Rad = yRad
         y = abs(y)
-        if not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes):
-            f = minifactor(y, P, smooth_base)
+        if y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes):
+            proc += 1 
+            f = minifactor4(y, P, smooth_base)
             if f != None:            
                 if f[1] == 1:   
                     filtered = filter_out_even_powers(f[0])
                     rel = [filtered, y, Rad, A] 
+                    #print(rel)
                     Rels.append(rel)
                 elif f[1] in partials:
                     a = partials[f[1]]
@@ -441,9 +421,13 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
                         sys.stderr.write("Found cycle with partial\n")
                     else:
                         Rels.append([p, rhs, lhs, Ahs])
+                    #with merged_count.value.get_lock():
                     merged_count.value += 1
                 else:
                     partials[f[1]] = [f[0], y, Rad, A]
+        else:
+            noproc += 1
+
         if i % m == 0:
             lRels = len(Rels)
             lt = time.time()
@@ -451,6 +435,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             ltd = lt
             if i > 0:
                 if lRels > 0:
+                    #eta = td * (((ld / m) + (required_relations / lRels)) / 2)
                     eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
                 else:
                     eta = int(td * (ld / m))
@@ -458,14 +443,17 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             else:
                 etas = "No ETA estimated"
             tds = humanfriendly.format_timespan(td)
-            msg = "[%d] relations_find: range (%d, %d , %d), found: %d of %d, merged: %d, iter_elapsed: %s, eta: %s, alive tasks: %d.\n" % (pid, start,i,ld,lRels,required_relations, merged_count.value,tds,etas,tasks.value)
+            msg = "[%d] relations_find: range (%d, %d , %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, iter_elapsed: %s, eta: %s, alive tasks: %d.\n" 
+            msg = msg % (pid, start, i, ld, lRels, required_relations, merged_count.value, proc, noproc , tds, etas, tasks.value)
             sys.stderr.write(msg)
     
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
-    msg = "[%d] relations_find: Ended range(%d, %d, %d), found: %d of %d, merged: %d, time elapsed: %s, alive tasks: %d.\n" % (pid, start,i,ld,len(Rels),required_relations, merged_count.value, tds, tasks.value)
+    msg = "[%d] relations_find: Ended range(%d, %d, %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, time elapsed: %s, alive tasks: %d.\n" 
+    msg = msg % (pid, start,i,ld,len(Rels),required_relations, merged_count.value, proc, noproc, tds, tasks.value)
     sys.stderr.write(msg)
     tasks.value -= 1
+    #Rels += Found_Rels
     return Found_Rels
     
 
@@ -500,7 +488,6 @@ def Gaussian_elimination_GF2(A):
           break
   return marks, A
 
-
 def left_nullspace(A):
     """
     Compute left null space:
@@ -508,6 +495,7 @@ def left_nullspace(A):
     """
     A = transpose(A)
     marks, A = Gaussian_elimination_GF2(A)
+    #marks, A = gauss(A)
 
     B = []
     for row in range(len(A)):
@@ -679,18 +667,20 @@ def _MPQS(N, verbose=True, M = 1):
     Rels = []
 
     T = cpu_count()
+
     B2, _ , B1 = prebuilt_params(log2(N))
     Prime_base, log_p = find_primebase(N, B2)
-    B1 //= M
+    #B1 //= M
     B2 = len(Prime_base)
   
     if B2 == 1:
         sys.stderr.write("Found small factor: %d\n" % Prime_base[0])
-        return Prime_base + _MPQS(N // Prime_base[0])
+        return Prime_base + _MPQS(N // Prime_base[0]), 0
 
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
     
+    #x_max = B2 *60  # size of the sieve
     x_max = B1
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
@@ -705,7 +695,7 @@ def _MPQS(N, verbose=True, M = 1):
 
     start = 0
     stop = B1 # range to sieve
-
+    
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
     cycleFactors = manager.list()
@@ -713,7 +703,7 @@ def _MPQS(N, verbose=True, M = 1):
     merged_count = manager.Value("i", 0)
     tasks = manager.Value("i", 0)
 
-    min_search = 0
+    min_poly_search = 0
     polys_ABC = []
 
     while True:
@@ -726,29 +716,37 @@ def _MPQS(N, verbose=True, M = 1):
 
         t1 = time.time()
         sys.stderr.write("Data collection with %d threads...\n" % T)
-  
-        polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_search, polys_ABC) # generate n distinct polys one for each cpu core.
-        if len(early_factors) > 0:
-            tmp = 1
-            small = []
-            for early_factor in early_factors:
-                g = gcd(early_factor, N)
-                if N > g > 1:
-                    tmp *= g
-                    small.append(g)
-            if tmp > 1:
-                return small + MPQS(N // tmp)
 
+        need_more_polys = (len(Rels) == 0) # if we havent still got any relation we need another set of polys.
+
+        if need_more_polys == True:  
+            polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys_ABC) # generate n distinct polys one for each cpu core.
+            if len(early_factors) > 0:
+                tmp = 1
+                small = []
+                for early_factor in early_factors:
+                    g = gcd(early_factor, N)
+                    if N > g > 1:
+                        tmp *= g
+                        small.append(g)
+                if tmp > 1:
+                    return small + MPQS(N // tmp), 0
+             
         inputs = [] 
         taskid = 1
         # generate tasks parameters
-        for poly in polys:
+        for poly in polys[0 - (T * 2):len(polys)]:
+            #print(poly.sums)
+            #calc_sums
             inputs += [(taskid, Nm, start, stop, Prime_base, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, poly)]
+        #sys.exit(0)
             taskid += 1
 
         # deploy tasks to every cpu core.
         pols = []
         workpool = Pool(T)
+        #print(T,len(inputs))
+        #sys.exit(0)
         with workpool:
             R = workpool.starmap(relations_find, inputs)  
         workpool.close()
@@ -776,11 +774,13 @@ def _MPQS(N, verbose=True, M = 1):
             sys.stderr.write("Done in: %f secs.\n" % (t4-t3))
             
             if result != None:
-                return result
+                return result, len(polys)
 
-        #start = stop
-        #stop += B1
-        min_search += T
+        if not need_more_polys:
+            start = stop
+            stop += B1
+        else:
+            min_poly_search += (T * 2)
 
         sys.stderr.write("Need to sieve %d differences more...\nNew sieving range %d:%d\n" % (required_relations - len(Rels),start,stop))
       
@@ -789,7 +789,10 @@ def MPQS(N):
     """ 
     Iterative version of MPQS. 
     """
-    result = _MPQS(N)
+
+    result, polys = _MPQS(N)
+    polycount = polys
+
     if result != None:
         R= []
         for r in result:
@@ -797,20 +800,23 @@ def MPQS(N):
                 sys.stderr.write("Found factor: %d\n" % r)
                 R.append(r)
             else:
-                R += _MPQS(r)
-        return R
+                r, polys = _MPQS(r)
+                R+=r
+                polycount += polys
+        return R, polycount
 
 
 if __name__ == "__main__":
     N = Integer(sys.argv[1])
     ts = time.time()
-    r = MPQS(N)
+    r, polycount = MPQS(N)
     Fp = trial_division(N,sorted(r))[0] 
     tmp = []
     for f,p in Fp:
         tmp += ["%d ^ %d" % (f,p)]
     msg = "The factorization for N = %d is: %s" % (N," * ".join(tmp))
     print(msg.replace(" ^ 1",""))
+    print("With %d distinct generator polys" % polycount)
     td = time.time() - ts
     tds = humanfriendly.format_timespan(td) 
     sys.stderr.write("All done in: %s.\n" % (tds))

From eeb7937a353b02296ff797938f32ac6eb38335a6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sat, 12 Jun 2021 01:49:26 -0300
Subject: [PATCH 31/56] Update

---
 MPQS.sage | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 4527842..2703557 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -364,6 +364,10 @@ def is_power(n, minprimes):
 
 
 def is_power_logprime(n, min_log_primes):
+    """
+    Given the precomputed logs of a set of primes 
+    we can test if a number n is a perfect power of that prime.
+    """
     ispow = False
     for p in min_log_primes:
         a = log(n) / p
@@ -403,12 +407,12 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             proc += 1 
             f = minifactor4(y, P, smooth_base)
             if f != None:            
-                if f[1] == 1:   
+                if f[1] == 1:  # found a relation 
                     filtered = filter_out_even_powers(f[0])
                     rel = [filtered, y, Rad, A] 
                     #print(rel)
                     Rels.append(rel)
-                elif f[1] in partials:
+                elif f[1] in partials: # found a partial and try to merge
                     a = partials[f[1]]
                     p = filter_out_even_powers(f[0] + a[0])
                     Ahs = A*a[3]
@@ -419,7 +423,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
                     if N > g > 1:
                         cycleFactors.append([g , N // g])
                         sys.stderr.write("Found cycle with partial\n")
-                    else:
+                    else: # didnt merge so store the partial for later merging
                         Rels.append([p, rhs, lhs, Ahs])
                     #with merged_count.value.get_lock():
                     merged_count.value += 1

From 58009825893b017e46c4cada70a9fab55901d2ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sun, 13 Jun 2021 00:15:03 -0300
Subject: [PATCH 32/56] add poly repr, new poly generation heuristic

---
 MPQS.sage | 41 +++++++++++++++++++++++++++++++----------
 1 file changed, 31 insertions(+), 10 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 2703557..8450cc6 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -176,7 +176,7 @@ def trial_division(n, P):
                 while r % p == 0:
                     pw += 1
                     r //= p
-               a.append((int(p),int(pw)))
+                a.append((int(p),int(pw)))
         else:
             break
     return a,r,n
@@ -260,10 +260,12 @@ class Poly:
         self.x_max = x_max
         self.search = search
         self.verbose = verbose
+        self.counter = 0
         self.early_factors = []
         self.create()
         self.solve_for_x()
  
+ 
     def create(self):
         """ 
         Construct our polynomial. 
@@ -313,18 +315,15 @@ class Poly:
 
     def solve_for_x(self):
         """ 
-        Get the mininmum x value of the polynomial where y=f(x) is a guaranteed relation. 
+        Get the minima of the polynomial where y=f(x) is a guaranteed relation. 
         Code borrowed https://github.com/cramppet/quadratic-sieve/blob/master/QS.py
         """        
         A = self.A
         B = self.B
         C = self.C
         n = self.n
-        start_vals = []
-        #for p in self.P:
         p = self.root_A
         g = 1
-        #for p in self.P:
         while g == 1:
             p = next_prime(p)
             ainv = 1
@@ -349,7 +348,14 @@ class Poly:
         Ax = A * x
         Rad = Ax + B 
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
-        
+      
+    def __repr__():
+        """
+        Return the string representation of the constructed polynomial
+        """
+        m = "F(X) = %d X ^ 2 + %d X + %d with minima: %d" % (self.A, self.B, self.C, self.minima) 
+        m = m.replace("+ -","- ")
+        return m
 
 def is_power(n, minprimes):
     ispow = False
@@ -397,6 +403,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
     msg = ""
     partials = {}
     i= 0
+    rels_found = 0
     for i in range(ld):
         if len(Rels) > required_relations or len(cycleFactors) > 0:
             break
@@ -408,6 +415,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             f = minifactor4(y, P, smooth_base)
             if f != None:            
                 if f[1] == 1:  # found a relation 
+                    rels_found += 1
                     filtered = filter_out_even_powers(f[0])
                     rel = [filtered, y, Rad, A] 
                     #print(rel)
@@ -433,6 +441,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             noproc += 1
 
         if i % m == 0:
+            poly.counter = rels_found
             lRels = len(Rels)
             lt = time.time()
             td = lt - ltd
@@ -451,6 +460,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             msg = msg % (pid, start, i, ld, lRels, required_relations, merged_count.value, proc, noproc , tds, etas, tasks.value)
             sys.stderr.write(msg)
     
+    poly.counter = rels_found
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
     msg = "[%d] relations_find: Ended range(%d, %d, %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, time elapsed: %s, alive tasks: %d.\n" 
@@ -594,8 +604,9 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
         else:
             pol_ABC = (pol.A,pol.B,pol.C)
             if pol_ABC not in polys_ABC:
-                m = "Found poly: f(x) = %dx^2 + %dx + %d\n" % (pol_ABC)
-                sys.stderr.write(m.replace("+ -","- "))
+                #m = "Found poly: f(x) = %d X ^ 2 + %d X + %d, minima: %d\n" % (pol_ABC[0], pol_ABC[1], pol_ABC[2], pol.minima)
+                m = repr(pol) 
+                sys.stderr.write(m)
                 polys_ABC.append(pol_ABC)
                 polys.append(pol)
                 cpolys += 1
@@ -603,6 +614,11 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
     return polys, polys_ABC, early_factors
 
 
+def poly_stats(polys):
+    for poly in polys:
+        sys.stderr.write("For poly: %s the count is: %d\n" % (repr(poly),poly.count))
+
+
 def prebuilt_params(bits):
     """
     Bounds estimation
@@ -709,6 +725,8 @@ def _MPQS(N, verbose=True, M = 1):
 
     min_poly_search = 0
     polys_ABC = []
+    last_lRels = 0
+    force_new_polys = False
 
     while True:
         # trim primes, recalc min
@@ -721,9 +739,12 @@ def _MPQS(N, verbose=True, M = 1):
         t1 = time.time()
         sys.stderr.write("Data collection with %d threads...\n" % T)
 
-        need_more_polys = (len(Rels) == 0) # if we havent still got any relation we need another set of polys.
+        lRels = len(Rels)
+        need_more_polys = not (lRels > last_lRels) or force_new_polys
+        last_lRels = lRels
 
-        if need_more_polys == True:  
+        if need_more_polys == True: 
+            sys.stderr.write("Need more polys...") 
             polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys_ABC) # generate n distinct polys one for each cpu core.
             if len(early_factors) > 0:
                 tmp = 1

From 554b0b857bb01eb9891749e56db84f6b6d23a23f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sun, 13 Jun 2021 02:34:12 -0300
Subject: [PATCH 33/56] hashable class polys, new generation algorithm

---
 MPQS.sage | 79 ++++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 55 insertions(+), 24 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 8450cc6..77ea5f1 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -260,7 +260,7 @@ class Poly:
         self.x_max = x_max
         self.search = search
         self.verbose = verbose
-        self.counter = 0
+        #self.counter = 0
         self.early_factors = []
         self.create()
         self.solve_for_x()
@@ -349,7 +349,7 @@ class Poly:
         Rad = Ax + B 
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
       
-    def __repr__():
+    def __repr__(self):
         """
         Return the string representation of the constructed polynomial
         """
@@ -357,6 +357,17 @@ class Poly:
         m = m.replace("+ -","- ")
         return m
 
+    def __hash__(self):
+        h = hash("%d-%d-%d" % (self.A,self.B,self.C))
+        #print("hash: %s" % h)
+        return h
+
+    def __eq__(self, other):
+        if isinstance(other, Poly):
+            return ((self.A == other.A) and (self.B == other.B) and (self.C == other.C))
+        else:
+            return NotImplemented
+
 def is_power(n, minprimes):
     ispow = False
     c = 0
@@ -383,7 +394,7 @@ def is_power_logprime(n, min_log_primes):
     return ispow
 
 
-def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, poly = None):
+def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
     """ 
     Relations search funcion 
     """
@@ -433,6 +444,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
                         sys.stderr.write("Found cycle with partial\n")
                     else: # didnt merge so store the partial for later merging
                         Rels.append([p, rhs, lhs, Ahs])
+                        rels_found += 1
                     #with merged_count.value.get_lock():
                     merged_count.value += 1
                 else:
@@ -441,7 +453,9 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             noproc += 1
 
         if i % m == 0:
-            poly.counter = rels_found
+            #print(repr(poly))
+            #poly.counter = rels_found
+            polycounts[poly] = rels_found
             lRels = len(Rels)
             lt = time.time()
             td = lt - ltd
@@ -460,7 +474,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             msg = msg % (pid, start, i, ld, lRels, required_relations, merged_count.value, proc, noproc , tds, etas, tasks.value)
             sys.stderr.write(msg)
     
-    poly.counter = rels_found
+    polycounts[poly] = rels_found
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
     msg = "[%d] relations_find: Ended range(%d, %d, %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, time elapsed: %s, alive tasks: %d.\n" 
@@ -587,13 +601,15 @@ def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
     return min_prime, thresh, fudge
 
 
-def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
+#def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
+def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
+
     """ 
     It searchs for distinct needed polys congruent to n. 
     """
     n = min_search
-    cpolys = 1
-    polys = []
+    cpolys = 0
+    #polys = []
     early_factors = []
     while cpolys <= needed:
         pol = Poly(N, Prime_base, x_max, search = n, verbose=False)
@@ -602,21 +618,26 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
                 if early_factor not in early_factors:
                     early_factors.append(early_factor)
         else:
-            pol_ABC = (pol.A,pol.B,pol.C)
-            if pol_ABC not in polys_ABC:
+            #pol_ABC = (pol.A,pol.B,pol.C)
+            #if pol_ABC not in polys_ABC:
+            if pol not in polys:
                 #m = "Found poly: f(x) = %d X ^ 2 + %d X + %d, minima: %d\n" % (pol_ABC[0], pol_ABC[1], pol_ABC[2], pol.minima)
                 m = repr(pol) 
-                sys.stderr.write(m)
-                polys_ABC.append(pol_ABC)
+                sys.stderr.write("New Poly %d: %s\n" % (cpolys,m))
+                #polys_ABC.append(pol_ABC)
                 polys.append(pol)
                 cpolys += 1
         n += 1
-    return polys, polys_ABC, early_factors
+    #return polys, polys_ABC, early_factors
+    return polys, early_factors
 
 
-def poly_stats(polys):
+
+def poly_stats(polys, polycounts):
     for poly in polys:
-        sys.stderr.write("For poly: %s the count is: %d\n" % (repr(poly),poly.count))
+        if poly in polycounts:
+            poly_count = polycounts[poly]
+            sys.stderr.write("For poly: %s the count is: %d\n" % (repr(poly),poly_count))
 
 
 def prebuilt_params(bits):
@@ -706,12 +727,10 @@ def _MPQS(N, verbose=True, M = 1):
     thresh = log10(m_val) * 0.735
     min_prime = int(thresh * 3)
     
-    required_relations_ratio = 1.05
-    required_relations = round(len(Prime_base) * required_relations_ratio) 
+
 
     sys.stderr.write("Factoring N: %d, bits: %d, digits: %d, B1: %d, B2: %d\n" % (N,bN,lN,B1,B2))
     sys.stderr.write("Multiplier is: %d\n" % multiplier)
-    sys.stderr.write("Need %d relations\n" % (required_relations))
 
     start = 0
     stop = B1 # range to sieve
@@ -719,17 +738,25 @@ def _MPQS(N, verbose=True, M = 1):
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
     cycleFactors = manager.list()
+    polycounts = manager.dict()
 
     merged_count = manager.Value("i", 0)
     tasks = manager.Value("i", 0)
 
     min_poly_search = 0
-    polys_ABC = []
+    #polys_ABC = []
+    polys = []
+
     last_lRels = 0
     force_new_polys = False
 
+    required_relations_ratio = 1.05
+    required_relations = len(Prime_base)
+
     while True:
         # trim primes, recalc min
+        required_relations = int(required_relations * required_relations_ratio)
+        sys.stderr.write("Need %d relations\n" % (required_relations))
         min_log_primes = [log(p) for p in Prime_base if p <= min_prime]
         Prime_base = [p for p in Prime_base if p > min_prime]
         log_primes = [log(p) for p in Prime_base]
@@ -744,8 +771,10 @@ def _MPQS(N, verbose=True, M = 1):
         last_lRels = lRels
 
         if need_more_polys == True: 
-            sys.stderr.write("Need more polys...") 
-            polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys_ABC) # generate n distinct polys one for each cpu core.
+            sys.stderr.write("Need more polys...\n") 
+            #polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys_ABC) # generate n distinct polys one for each cpu core.
+            polys, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys) # generate n distinct polys one for each cpu core.
+
             if len(early_factors) > 0:
                 tmp = 1
                 small = []
@@ -763,9 +792,9 @@ def _MPQS(N, verbose=True, M = 1):
         for poly in polys[0 - (T * 2):len(polys)]:
             #print(poly.sums)
             #calc_sums
-            inputs += [(taskid, Nm, start, stop, Prime_base, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, poly)]
-        #sys.exit(0)
+            inputs += [(taskid, Nm, start, stop, Prime_base, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
             taskid += 1
+        #sys.exit(0)
 
         # deploy tasks to every cpu core.
         pols = []
@@ -797,7 +826,9 @@ def _MPQS(N, verbose=True, M = 1):
             result = process_basis_vectors(Nm, basis, Rels, multiplier)
             t4 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t4-t3))
-            
+             
+            poly_stats(polys, polycounts)
+
             if result != None:
                 return result, len(polys)
 

From 1df68e0fdb5dfc715a091b7bdf1b5b21e7aeae61 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sun, 13 Jun 2021 02:38:12 -0300
Subject: [PATCH 34/56] clean up code

---
 MPQS.sage | 29 ++++++-----------------------
 1 file changed, 6 insertions(+), 23 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 77ea5f1..2ed8e1a 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -115,6 +115,7 @@ def choose_multiplier(n, prime_list):
             best_mult = mult_list[i]
     return best_mult;
 
+
 def mod_sqrt(n, p):
     """ 
     Tonelli shanks algorithm 
@@ -349,6 +350,7 @@ class Poly:
         Rad = Ax + B 
         return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
       
+
     def __repr__(self):
         """
         Return the string representation of the constructed polynomial
@@ -357,17 +359,20 @@ class Poly:
         m = m.replace("+ -","- ")
         return m
 
+
     def __hash__(self):
         h = hash("%d-%d-%d" % (self.A,self.B,self.C))
         #print("hash: %s" % h)
         return h
 
+
     def __eq__(self, other):
         if isinstance(other, Poly):
             return ((self.A == other.A) and (self.B == other.B) and (self.C == other.C))
         else:
             return NotImplemented
 
+
 def is_power(n, minprimes):
     ispow = False
     c = 0
@@ -429,7 +434,6 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
                     rels_found += 1
                     filtered = filter_out_even_powers(f[0])
                     rel = [filtered, y, Rad, A] 
-                    #print(rel)
                     Rels.append(rel)
                 elif f[1] in partials: # found a partial and try to merge
                     a = partials[f[1]]
@@ -453,8 +457,6 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             noproc += 1
 
         if i % m == 0:
-            #print(repr(poly))
-            #poly.counter = rels_found
             polycounts[poly] = rels_found
             lRels = len(Rels)
             lt = time.time()
@@ -462,7 +464,6 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             ltd = lt
             if i > 0:
                 if lRels > 0:
-                    #eta = td * (((ld / m) + (required_relations / lRels)) / 2)
                     eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
                 else:
                     eta = int(td * (ld / m))
@@ -481,7 +482,6 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
     msg = msg % (pid, start,i,ld,len(Rels),required_relations, merged_count.value, proc, noproc, tds, tasks.value)
     sys.stderr.write(msg)
     tasks.value -= 1
-    #Rels += Found_Rels
     return Found_Rels
     
 
@@ -516,6 +516,7 @@ def Gaussian_elimination_GF2(A):
           break
   return marks, A
 
+
 def left_nullspace(A):
     """
     Compute left null space:
@@ -601,15 +602,12 @@ def recalculate_min_prime_thresh(thresh, Prime_base, log_p):
     return min_prime, thresh, fudge
 
 
-#def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys_ABC=[]):
 def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
-
     """ 
     It searchs for distinct needed polys congruent to n. 
     """
     n = min_search
     cpolys = 0
-    #polys = []
     early_factors = []
     while cpolys <= needed:
         pol = Poly(N, Prime_base, x_max, search = n, verbose=False)
@@ -618,17 +616,12 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
                 if early_factor not in early_factors:
                     early_factors.append(early_factor)
         else:
-            #pol_ABC = (pol.A,pol.B,pol.C)
-            #if pol_ABC not in polys_ABC:
             if pol not in polys:
-                #m = "Found poly: f(x) = %d X ^ 2 + %d X + %d, minima: %d\n" % (pol_ABC[0], pol_ABC[1], pol_ABC[2], pol.minima)
                 m = repr(pol) 
                 sys.stderr.write("New Poly %d: %s\n" % (cpolys,m))
-                #polys_ABC.append(pol_ABC)
                 polys.append(pol)
                 cpolys += 1
         n += 1
-    #return polys, polys_ABC, early_factors
     return polys, early_factors
 
 
@@ -721,13 +714,10 @@ def _MPQS(N, verbose=True, M = 1):
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
     
-    #x_max = B2 *60  # size of the sieve
     x_max = B1
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
     min_prime = int(thresh * 3)
-    
-
 
     sys.stderr.write("Factoring N: %d, bits: %d, digits: %d, B1: %d, B2: %d\n" % (N,bN,lN,B1,B2))
     sys.stderr.write("Multiplier is: %d\n" % multiplier)
@@ -744,7 +734,6 @@ def _MPQS(N, verbose=True, M = 1):
     tasks = manager.Value("i", 0)
 
     min_poly_search = 0
-    #polys_ABC = []
     polys = []
 
     last_lRels = 0
@@ -772,7 +761,6 @@ def _MPQS(N, verbose=True, M = 1):
 
         if need_more_polys == True: 
             sys.stderr.write("Need more polys...\n") 
-            #polys, polys_ABC, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys_ABC) # generate n distinct polys one for each cpu core.
             polys, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys) # generate n distinct polys one for each cpu core.
 
             if len(early_factors) > 0:
@@ -790,17 +778,12 @@ def _MPQS(N, verbose=True, M = 1):
         taskid = 1
         # generate tasks parameters
         for poly in polys[0 - (T * 2):len(polys)]:
-            #print(poly.sums)
-            #calc_sums
             inputs += [(taskid, Nm, start, stop, Prime_base, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
             taskid += 1
-        #sys.exit(0)
 
         # deploy tasks to every cpu core.
         pols = []
         workpool = Pool(T)
-        #print(T,len(inputs))
-        #sys.exit(0)
         with workpool:
             R = workpool.starmap(relations_find, inputs)  
         workpool.close()

From dbbad02600acc001a84726a87632fe26e9c3a320 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sun, 13 Jun 2021 02:48:44 -0300
Subject: [PATCH 35/56] Rearrange functions, add missing docstrings

---
 MPQS.sage | 182 ++++++++++++++++++++++++++++--------------------------
 1 file changed, 93 insertions(+), 89 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 2ed8e1a..ff94dd7 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -23,6 +23,56 @@ from itertools import repeat
 import humanfriendly
 
 
+def prebuilt_params(bits):
+    """
+    Bounds estimation
+    borrowed from msieve/mpqs.c
+    """
+    if bits <= 64:
+        return [100, 40, 1 * 65536]
+    if bits <= 128: 
+        return [450, 40, 1 * 65536]
+    if bits <= 183:
+        return [2000, 40, 1 * 65536]
+    if bits <= 200: 
+        return [3000, 50, 1 * 65536]
+    if bits <= 212: 
+        return [5400, 50, 3 * 65536]
+    if bits <= 233:
+        return [10000, 100, 3 * 65536]
+    if bits <= 249:
+        return [27000, 100, 3 * 65536]
+    if bits <= 266:
+        return [50000, 100, 3 * 65536]
+    if bits <= 283:
+        return [55000, 80, 3 * 65536]
+    if bits <= 298:
+        return [60000, 80, 9 * 65536]
+    if bits <= 315:
+        return [80000, 150, 9 * 65536]
+    if bits <= 332:
+        return [100000, 150, 9 * 65536]
+    if bits <= 348:
+        return [140000, 150, 9 * 65536]
+    if bits <= 363:
+        return [210000, 150, 13 * 65536]
+    if bits <= 379:
+        return [300000, 150, 17 * 65536]
+    if bits <= 395:
+        return [400000, 150, 21 * 65536]
+    if bits <= 415:
+        return [500000, 150, 25 * 65536] # beyond this point you're crazy 
+    if bits <= 440:
+        return [700000, 150, 33 * 65536]
+    if bits <= 465:
+        return [900000, 150, 50 * 65536]
+    if bits <= 490:
+        return [1100000, 150, 75 * 65536]
+    if bits <= 512:
+        return [1300000, 150, 100 * 65536]
+    return [1300000, 150, 100 * 65536]
+
+
 def choose_multiplier(n, prime_list):
     """
     Code borrowed from msieve/mpqs.c
@@ -152,16 +202,6 @@ def mod_sqrt(n, p):
         return a  # p == 2
 
 
-def is_power_logprime(n, min_log_primes):
-    ispow = False
-    for p in min_log_primes:
-        a = log(n) / p
-        b = int(a)
-        if a == b:
-            iwpow = True
-    return ispow
-
-
 def trial_division(n, P):
     """ 
     A simple trial division, factors are given by P-list. 
@@ -184,6 +224,9 @@ def trial_division(n, P):
 
 
 def merge_powers(ppws):
+    """
+    Merge powers such that: (a^b) * (a^c) == a^(b+c)
+    """
     d = {}
     for p,pw in ppws:
         if p not in d:
@@ -198,7 +241,7 @@ def merge_powers(ppws):
 
 def filter_out_even_powers(ppws):
     """ 
-    Filter out even powers. 
+    Same as merge but Filter out even powers. 
     """
     d = {}
     for p,pw in ppws:
@@ -233,6 +276,35 @@ def is_smooth(x, P):
     return abs(y) == 1 
 
 
+def is_power(n, primes):
+    """
+    Given a interger n and a set of primes it checks if n is a perfect power.
+    """
+    ispow = False
+    c = 0
+    for p in primes:
+        a = log(n) / log(p)
+        b = int(a)
+        if a == b:
+            ispow = True
+            c = b
+    return ispow
+
+
+def is_power_logprime(n, log_primes):
+    """
+    Given the precomputed logs of a set of primes 
+    we can test if a number n is a perfect power of that prime.
+    """
+    ispow = False
+    for p in log_primes:
+        a = log(n) / p
+        b = int(a)
+        if a == b:
+            iwpow = True
+    return ispow
+
+
 def minifactor4(x, P, smooth_base):
     """ 
     Minifactor. 
@@ -361,44 +433,24 @@ class Poly:
 
 
     def __hash__(self):
+        """
+        Hashes unique values of the polynomial to construct an python internal id.
+        """
         h = hash("%d-%d-%d" % (self.A,self.B,self.C))
         #print("hash: %s" % h)
         return h
 
 
     def __eq__(self, other):
+        """
+        Internal python facility needed to check if two polynomials are the same.
+        """
         if isinstance(other, Poly):
             return ((self.A == other.A) and (self.B == other.B) and (self.C == other.C))
         else:
             return NotImplemented
 
 
-def is_power(n, minprimes):
-    ispow = False
-    c = 0
-    for p in minprimes:
-        a = log(n) / log(p)
-        b = int(a)
-        if a == b:
-            ispow = True
-            c = b
-    return ispow
-
-
-def is_power_logprime(n, min_log_primes):
-    """
-    Given the precomputed logs of a set of primes 
-    we can test if a number n is a perfect power of that prime.
-    """
-    ispow = False
-    for p in min_log_primes:
-        a = log(n) / p
-        b = int(a)
-        if a == b:
-            iwpow = True
-    return ispow
-
-
 def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
     """ 
     Relations search funcion 
@@ -625,62 +677,14 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
     return polys, early_factors
 
 
-
 def poly_stats(polys, polycounts):
+    """
+    Prints polynomials relations stats
+    """
     for poly in polys:
         if poly in polycounts:
             poly_count = polycounts[poly]
-            sys.stderr.write("For poly: %s the count is: %d\n" % (repr(poly),poly_count))
-
-
-def prebuilt_params(bits):
-    """
-    Bounds estimation
-    borrowed from msieve/mpqs.c
-    """
-    if bits <= 64:
-        return [100, 40, 1 * 65536]
-    if bits <= 128: 
-        return [450, 40, 1 * 65536]
-    if bits <= 183:
-        return [2000, 40, 1 * 65536]
-    if bits <= 200: 
-        return [3000, 50, 1 * 65536]
-    if bits <= 212: 
-        return [5400, 50, 3 * 65536]
-    if bits <= 233:
-        return [10000, 100, 3 * 65536]
-    if bits <= 249:
-        return [27000, 100, 3 * 65536]
-    if bits <= 266:
-        return [50000, 100, 3 * 65536]
-    if bits <= 283:
-        return [55000, 80, 3 * 65536]
-    if bits <= 298:
-        return [60000, 80, 9 * 65536]
-    if bits <= 315:
-        return [80000, 150, 9 * 65536]
-    if bits <= 332:
-        return [100000, 150, 9 * 65536]
-    if bits <= 348:
-        return [140000, 150, 9 * 65536]
-    if bits <= 363:
-        return [210000, 150, 13 * 65536]
-    if bits <= 379:
-        return [300000, 150, 17 * 65536]
-    if bits <= 395:
-        return [400000, 150, 21 * 65536]
-    if bits <= 415:
-        return [500000, 150, 25 * 65536] # beyond this point you're crazy 
-    if bits <= 440:
-        return [700000, 150, 33 * 65536]
-    if bits <= 465:
-        return [900000, 150, 50 * 65536]
-    if bits <= 490:
-        return [1100000, 150, 75 * 65536]
-    if bits <= 512:
-        return [1300000, 150, 100 * 65536]
-    return [1300000, 150, 100 * 65536]
+            sys.stderr.write("For polynomial: %s the relations count is: %d\n" % (repr(poly),poly_count))
 
 
 def _MPQS(N, verbose=True, M = 1):

From 25c3e3e97265a879280f1fd54538c14286ed0d96 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Mon, 14 Jun 2021 00:34:36 -0300
Subject: [PATCH 36/56] Complete Murphy(1999) plog check threshold, add to poly
 class __add__, __sub__, __mul__, add create_matrix

---
 MPQS.sage | 208 ++++++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 164 insertions(+), 44 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index ff94dd7..9f36c0e 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -21,7 +21,7 @@ from sage.parallel.multiprocessing_sage import parallel_iter
 from multiprocessing import cpu_count, Pool, Manager
 from itertools import repeat
 import humanfriendly
-
+from copy import *
 
 def prebuilt_params(bits):
     """
@@ -238,7 +238,6 @@ def merge_powers(ppws):
        tmp2.append((p,d[p]))
     return tmp2
 
-
 def filter_out_even_powers(ppws):
     """ 
     Same as merge but Filter out even powers. 
@@ -326,7 +325,7 @@ class Poly:
     Quadratic polynomial helper class: 
     type Ax + 2Bx + C 
     """
-    def __init__(self, n, P, x_max, search = 0, verbose = None, logp = None):
+    def __init__(self, n, P, x_max, search = 0, verbose = None, logp = None, A = None, B = None, C = None):
         self.n = int(n)
         self.P = P
         self.logp = logp
@@ -335,7 +334,18 @@ class Poly:
         self.verbose = verbose
         #self.counter = 0
         self.early_factors = []
-        self.create()
+
+        if A == None and B == None and C == None:
+            self.create()
+        else:
+            print("given poly parameters: %d %d %d" % (A,B,C))
+            self.A = A
+            self.B = B
+            self.C = C
+            root_2n = isqrt(2 * self.n)
+            root_A = next_prime(isqrt(root_2n // x_max))
+            self.root_A = root_A
+
         self.solve_for_x()
  
  
@@ -396,19 +406,24 @@ class Poly:
         C = self.C
         n = self.n
         p = self.root_A
-        g = 1
-        while g == 1:
-            p = next_prime(p)
-            ainv = 1
-            if A != 1:
-                g, inv, _ = gcdext(A, p)
-                if g != 1:
-                    r1 = mod_sqrt(C, p)
-                    r2 = (-1 * r1) % p
-                    start1 = (ainv * (r1 - B)) % p
-                    start2 = (ainv * (r2 - B)) % p
-                    self.X = [int(start1), int(start2)]
-                    self.minima = sum(self.X)//2
+        if A > 1:
+            g = 1
+            while g == 1:
+                p = next_prime(p)
+                ainv = 1
+                if A != 1:
+                    g, inv, _ = gcdext(A, p)
+                    if g != 1:
+                        r1 = mod_sqrt(C, p)
+                        r2 = (-1 * r1) % p
+                        start1 = (ainv * (r1 - B)) % p
+                        start2 = (ainv * (r2 - B)) % p
+                        self.X = [int(start1), int(start2)]
+                        self.minima = sum(self.X)//2
+        elif A == 1:
+            i = isqrt(abs(C))
+            self.X = [-i,i]
+            self.minima = 0
 
 
     def eval(self, x):
@@ -420,7 +435,8 @@ class Poly:
         C = self.C
         Ax = A * x
         Rad = Ax + B 
-        return (Rad + B) * x + C, Rad # same as (Ax + 2B)x + C
+        y = (Rad + B) * x + C
+        return y, Rad, x # same as (Ax + 2B)x + C
       
 
     def __repr__(self):
@@ -431,6 +447,37 @@ class Poly:
         m = m.replace("+ -","- ")
         return m
 
+    def __add__(self, other):
+        """
+        Adds one polynomial to other
+        """
+        if isinstance(other, Poly):
+            if self.n == other.n and self.x_max == olther.x_max:
+                A, B, C = self.A + other.A, self.B + other.B, self.C + other.C
+                return Poly(self.n, self.x_max, search = None, A = A, B = B, C = C)
+        else:
+            return NotImplemented
+
+    def __sub__(self, other):
+        """
+        Substracts one polynomial to other
+        """
+        if isinstance(other, Poly):
+            if self.n == other.n and self.x_max == olther.x_max:
+                A, B, C = self.A - other.A, self.B - other.B, self.C - other.C
+                return Poly(self.n, self.x_max, search = None, A = A, B = B, C = C)
+        else:
+            return NotImplemented
+
+    def __mul__(self, other):
+        """
+        Scalar multiply: scalar*(F(x))
+        """
+        if isinstance(other, int):
+            A, B, C = self.A * other, self.B * other, self.C * other
+            return Poly(self.n, self.x_max, search = None, A = A, B = B, C = C)
+        else:
+            return NotImplemented
 
     def __hash__(self):
         """
@@ -451,22 +498,39 @@ class Poly:
             return NotImplemented
 
 
-def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
+def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
     """ 
     Relations search funcion 
     """
     #pid = os.getpid()
     pid = taskid
     tasks.value += 1
-    sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d sieving start\n" % (pid,start, stop, (stop-start)))
+    sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d sieving start, tresh: %f.\n" % (pid,start, stop, (stop-start), thresh))
     m = poly.minima
-    Diffs = [(poly.eval(abs(x)),abs(x)) for x in range(start + m, stop +m)]
+    Diffs = [poly.eval(abs(x)) for x in range(start + m, stop +m)]
+    
+
+    ld = len(Diffs)
+    pre_log_filter = True
+    if pre_log_filter:
+        """
+        In theory this should make computations faster as we skeep to factor y=f(x) that does not shield any useful relation,
+        but in the practice is not working.
+        """
+        logs_y = [0] * (stop - start) * 2
+        for j in range(len(P)):
+            p = P[j]
+            log_p = log_primes[j]
+            for i in range(start + m , stop + m, p):
+                logs_y[i - m] += log_p
+        sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d log_filtering ended\n" % (pid,start, stop, (stop-start)))    
+
     proc = noproc = 0
     Found_Rels = []
     A = poly.A
     st = time.time()
     ltd = st
-    ld = len(Diffs)
+
     m = 1000
     msg = ""
     partials = {}
@@ -474,12 +538,19 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
     rels_found = 0
     for i in range(ld):
         if len(Rels) > required_relations or len(cycleFactors) > 0:
+            polycounts[poly] = rels_found
             break
-        yRad, x = Diffs[i]
-        y, Rad = yRad
+        y, Rad, x = Diffs[i]
         y = abs(y)
-        if y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes):
-            proc += 1 
+
+        if pre_log_filter:
+           candidate = (0 < logs_y[i - m] <= thresh and not is_prime(y))
+           a, b, c = is_square(y), is_prime(y), is_power_logprime(y, min_log_primes)
+        else:
+           candidate = (y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes))
+
+        if candidate:
+            proc += 1
             f = minifactor4(y, P, smooth_base)
             if f != None:            
                 if f[1] == 1:  # found a relation 
@@ -509,7 +580,6 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, smooth_base, Rels,
             noproc += 1
 
         if i % m == 0:
-            polycounts[poly] = rels_found
             lRels = len(Rels)
             lt = time.time()
             td = lt - ltd
@@ -564,11 +634,27 @@ def Gaussian_elimination_GF2(A):
           marks[i] = True
           for k in range(j+1,m):
             if A[i][k] == 1:
-              A[i][k] = (A[i][j] ^ A[i][k]) 
+              A[i][k] = (A[i][j] + A[i][k]) % 2 
           break
   return marks, A
 
 
+def create_matrix(Rels, P):
+    """
+    Create a matrix with the relations.
+    """
+    M1 = []
+    for i in range(0,len(Rels)):
+        row = []
+        for j in range(0,len(P)): 
+            if P[j] in Rels[i][0]:
+                row.append(1)
+            else:
+                row.append(0)
+        M1.append(row)
+    return M1
+
+
 def left_nullspace(A):
     """
     Compute left null space:
@@ -579,9 +665,9 @@ def left_nullspace(A):
     #marks, A = gauss(A)
 
     B = []
-    for row in range(len(A)):
-        if marks[row] == False:
-            B.append(A[row])
+    for row_index in range(0,len(A)):
+        if marks[row_index] == False:
+            B.append(A[row_index])
     return B
 
 
@@ -676,6 +762,27 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
         n += 1
     return polys, early_factors
 
+def getbestpolys(polys, polycounts, T):
+    """
+    Experimental function, break things.
+    """
+    tmp_polys = copy(polys)
+    new_polys = []
+    print("have polys: %d" % len(tmp_polys))
+    print(polys)
+    while len(new_polys) < T:
+        max_poly_count = 0
+        for i in range(0,len(tmp_polys)):
+            poly = tmp_polys[i]
+            if poly in polycounts:
+                poly_count = polycounts[poly]
+                if poly_count > max_poly_count:
+                    max_poly_count = poly_count
+                    best_poly = poly
+                    del polys[i]
+        print("new best poly: %s" % repr(poly))
+        new_polys.append(poly)
+    return new_polys
 
 def poly_stats(polys, polycounts):
     """
@@ -692,19 +799,20 @@ def _MPQS(N, verbose=True, M = 1):
     Main MPQS function. 
     """
     if is_prime(N):
-        return [N]
+        return [N], 0 
 
     bN, lN = int(log2(N)), len(str(N))
     i2N = isqrt(N)
 
     if pow(i2N,2) == N:
-        return [i2N, i2N]
+        return [i2N, i2N],0
 
     i2Np1 = i2N + 1 
     root_2n = isqrt(2*N)
     Rels = []
 
     T = cpu_count()
+    T *= M
 
     B2, _ , B1 = prebuilt_params(log2(N))
     Prime_base, log_p = find_primebase(N, B2)
@@ -713,7 +821,8 @@ def _MPQS(N, verbose=True, M = 1):
   
     if B2 == 1:
         sys.stderr.write("Found small factor: %d\n" % Prime_base[0])
-        return Prime_base + _MPQS(N // Prime_base[0]), 0
+        r, polycount = _MPQS(N // Prime_base[0])
+        return Prime_base + r, polycount
 
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
@@ -751,10 +860,10 @@ def _MPQS(N, verbose=True, M = 1):
         required_relations = int(required_relations * required_relations_ratio)
         sys.stderr.write("Need %d relations\n" % (required_relations))
         min_log_primes = [log(p) for p in Prime_base if p <= min_prime]
-        Prime_base = [p for p in Prime_base if p > min_prime]
+        filtered_Prime_base = [p for p in Prime_base if p > min_prime]
         log_primes = [log(p) for p in Prime_base]
         smooth_base = prod(Prime_base)
-        min_prime, thresh, fudge = recalculate_min_prime_thresh(thresh, Prime_base, log_p)
+        min_prime, thresh, fudge = recalculate_min_prime_thresh(thresh, filtered_Prime_base, log_p)
 
         t1 = time.time()
         sys.stderr.write("Data collection with %d threads...\n" % T)
@@ -765,7 +874,7 @@ def _MPQS(N, verbose=True, M = 1):
 
         if need_more_polys == True: 
             sys.stderr.write("Need more polys...\n") 
-            polys, early_factors = generate_polys(Nm, Prime_base, x_max, T * 2, min_poly_search, polys) # generate n distinct polys one for each cpu core.
+            polys, early_factors = generate_polys(Nm, Prime_base, x_max, T, min_poly_search, polys) # generate n distinct polys one for each cpu core.
 
             if len(early_factors) > 0:
                 tmp = 1
@@ -780,9 +889,13 @@ def _MPQS(N, verbose=True, M = 1):
              
         inputs = [] 
         taskid = 1
+        #polys += [Poly(N, Prime_base, x_max, search = None, verbose=False, A=1, B=0, C=-N)] # lower shield base poly for QS.
+
         # generate tasks parameters
-        for poly in polys[0 - (T * 2):len(polys)]:
-            inputs += [(taskid, Nm, start, stop, Prime_base, min_log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
+        #polys = getbestpolys(polys, polycounts, T)
+        for poly in polys[0 - (T * M):len(polys)]:
+        #for poly in polys:
+            inputs += [(taskid, Nm, start, stop, filtered_Prime_base, min_log_primes, log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
             taskid += 1
 
         # deploy tasks to every cpu core.
@@ -796,10 +909,17 @@ def _MPQS(N, verbose=True, M = 1):
         t2 = time.time()
         sys.stderr.write("Done in: %f secs.\n" % (t2-t1))
         sys.stderr.write("Found %d rels with %d base primes.\n" % (len(Rels),len(Prime_base)))
+        #sys.stderr.write("Found %d rels with %d base primes.\n Sorting..." % (len(Rels),len(Prime_base)))
+
 
         if len(cycleFactors) > 0:
-            return cycleFactors
+            return cycleFactors, len(polys)
+
+        #Rels = list(set(Rels))
+        t3 = time.time()
+        #sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
 
+        
         # when needed relations is reached proceed to do linear algebra
         if len(Rels) > required_relations:
             sys.stderr.write("Found %d enough relations of %d needed relations...\n" % (len(Rels),required_relations))
@@ -807,12 +927,12 @@ def _MPQS(N, verbose=True, M = 1):
             
             basis = linear_algebra(Rels, Prime_base) # calculate left nullspace of a GF(2) matrix.
 
-            t3 = time.time()
-            sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
-            sys.stderr.write("Matrix reduction...\n")
-            result = process_basis_vectors(Nm, basis, Rels, multiplier)
             t4 = time.time()
             sys.stderr.write("Done in: %f secs.\n" % (t4-t3))
+            sys.stderr.write("Matrix reduction...")
+            result = process_basis_vectors(Nm, basis, Rels, multiplier)
+            t5 = time.time()
+            sys.stderr.write("Done in: %f secs.\n" % (t5-t4))
              
             poly_stats(polys, polycounts)
 

From 9b484e341d57eb2b2df615da87e6a1fd6d914f1e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Mon, 14 Jun 2021 18:59:29 -0300
Subject: [PATCH 37/56] take logp table to outside relations_find, cleanup code

---
 MPQS.sage | 72 +++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 46 insertions(+), 26 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 9f36c0e..2cad2e8 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -23,6 +23,7 @@ from itertools import repeat
 import humanfriendly
 from copy import *
 
+
 def prebuilt_params(bits):
     """
     Bounds estimation
@@ -238,6 +239,7 @@ def merge_powers(ppws):
        tmp2.append((p,d[p]))
     return tmp2
 
+
 def filter_out_even_powers(ppws):
     """ 
     Same as merge but Filter out even powers. 
@@ -447,6 +449,7 @@ class Poly:
         m = m.replace("+ -","- ")
         return m
 
+
     def __add__(self, other):
         """
         Adds one polynomial to other
@@ -458,6 +461,7 @@ class Poly:
         else:
             return NotImplemented
 
+
     def __sub__(self, other):
         """
         Substracts one polynomial to other
@@ -469,6 +473,7 @@ class Poly:
         else:
             return NotImplemented
 
+
     def __mul__(self, other):
         """
         Scalar multiply: scalar*(F(x))
@@ -479,6 +484,7 @@ class Poly:
         else:
             return NotImplemented
 
+
     def __hash__(self):
         """
         Hashes unique values of the polynomial to construct an python internal id.
@@ -498,7 +504,22 @@ class Poly:
             return NotImplemented
 
 
-def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
+def compute_logs_y_sums(start, stop, P, log_primes):
+     """
+     Compute log table
+     """
+     logs_y = [0] * ((stop - start) + 1)
+     for j in range(len(P)):
+         p = P[j]
+         log_p = log_primes[j]
+         #for i in range(start + m, stop + m, p):
+         #    logs_y[i - m] += log_p
+         for i in range(start , stop , p):
+             logs_y[i] += log_p
+     return logs_y
+
+
+def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
     """ 
     Relations search funcion 
     """
@@ -512,18 +533,6 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, smooth
 
     ld = len(Diffs)
     pre_log_filter = True
-    if pre_log_filter:
-        """
-        In theory this should make computations faster as we skeep to factor y=f(x) that does not shield any useful relation,
-        but in the practice is not working.
-        """
-        logs_y = [0] * (stop - start) * 2
-        for j in range(len(P)):
-            p = P[j]
-            log_p = log_primes[j]
-            for i in range(start + m , stop + m, p):
-                logs_y[i - m] += log_p
-        sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d log_filtering ended\n" % (pid,start, stop, (stop-start)))    
 
     proc = noproc = 0
     Found_Rels = []
@@ -544,7 +553,9 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, smooth
         y = abs(y)
 
         if pre_log_filter:
-           candidate = (0 < logs_y[i - m] <= thresh and not is_prime(y))
+           #x = i
+           #print(x)
+           candidate = (0 < logs_y[i] <= thresh and not is_prime(y))
            a, b, c = is_square(y), is_prime(y), is_power_logprime(y, min_log_primes)
         else:
            candidate = (y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes))
@@ -754,7 +765,7 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
                 if early_factor not in early_factors:
                     early_factors.append(early_factor)
         else:
-            if pol not in polys:
+            if pol not in polys and pol.minima > 0:
                 m = repr(pol) 
                 sys.stderr.write("New Poly %d: %s\n" % (cpolys,m))
                 polys.append(pol)
@@ -762,6 +773,7 @@ def generate_polys(N, Prime_base, x_max, needed, min_search = 0, polys=[]):
         n += 1
     return polys, early_factors
 
+
 def getbestpolys(polys, polycounts, T):
     """
     Experimental function, break things.
@@ -784,6 +796,7 @@ def getbestpolys(polys, polycounts, T):
         new_polys.append(poly)
     return new_polys
 
+
 def poly_stats(polys, polycounts):
     """
     Prints polynomials relations stats
@@ -794,7 +807,7 @@ def poly_stats(polys, polycounts):
             sys.stderr.write("For polynomial: %s the relations count is: %d\n" % (repr(poly),poly_count))
 
 
-def _MPQS(N, verbose=True, M = 1):
+def _MPQS(N, verbose=True, M = 2):
     """ 
     Main MPQS function. 
     """
@@ -840,6 +853,7 @@ def _MPQS(N, verbose=True, M = 1):
     
     manager = Manager()
     Rels = manager.list() # placeholder list for relations shareable between child processes.
+    newlist = manager.list()
     cycleFactors = manager.list()
     polycounts = manager.dict()
 
@@ -855,13 +869,23 @@ def _MPQS(N, verbose=True, M = 1):
     required_relations_ratio = 1.05
     required_relations = len(Prime_base)
 
+    logs_y = []
+
     while True:
         # trim primes, recalc min
-        required_relations = int(required_relations * required_relations_ratio)
-        sys.stderr.write("Need %d relations\n" % (required_relations))
+
+        sys.stderr.write("Generating primebase and logprime table...\n")
         min_log_primes = [log(p) for p in Prime_base if p <= min_prime]
         filtered_Prime_base = [p for p in Prime_base if p > min_prime]
         log_primes = [log(p) for p in Prime_base]
+        filtered_log_primes = [log(p) for p in filtered_Prime_base]
+      
+        if (stop-start) > len(logs_y):
+            logs_y = compute_logs_y_sums(0, stop-start, filtered_Prime_base, filtered_log_primes)
+
+        required_relations = int(required_relations * required_relations_ratio)
+        sys.stderr.write("Need %d relations\n" % (required_relations))
+        
         smooth_base = prod(Prime_base)
         min_prime, thresh, fudge = recalculate_min_prime_thresh(thresh, filtered_Prime_base, log_p)
 
@@ -887,15 +911,16 @@ def _MPQS(N, verbose=True, M = 1):
                 if tmp > 1:
                     return small + MPQS(N // tmp), 0
              
+        logs_y = compute_logs_y_sums(0, (stop-start), filtered_Prime_base, filtered_log_primes)
+
         inputs = [] 
         taskid = 1
         #polys += [Poly(N, Prime_base, x_max, search = None, verbose=False, A=1, B=0, C=-N)] # lower shield base poly for QS.
 
         # generate tasks parameters
-        #polys = getbestpolys(polys, polycounts, T)
         for poly in polys[0 - (T * M):len(polys)]:
         #for poly in polys:
-            inputs += [(taskid, Nm, start, stop, filtered_Prime_base, min_log_primes, log_primes, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
+            inputs += [(taskid, Nm, start, stop, filtered_Prime_base, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
             taskid += 1
 
         # deploy tasks to every cpu core.
@@ -909,17 +934,12 @@ def _MPQS(N, verbose=True, M = 1):
         t2 = time.time()
         sys.stderr.write("Done in: %f secs.\n" % (t2-t1))
         sys.stderr.write("Found %d rels with %d base primes.\n" % (len(Rels),len(Prime_base)))
-        #sys.stderr.write("Found %d rels with %d base primes.\n Sorting..." % (len(Rels),len(Prime_base)))
-
 
         if len(cycleFactors) > 0:
             return cycleFactors, len(polys)
 
-        #Rels = list(set(Rels))
         t3 = time.time()
-        #sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
 
-        
         # when needed relations is reached proceed to do linear algebra
         if len(Rels) > required_relations:
             sys.stderr.write("Found %d enough relations of %d needed relations...\n" % (len(Rels),required_relations))
@@ -943,7 +963,7 @@ def _MPQS(N, verbose=True, M = 1):
             start = stop
             stop += B1
         else:
-            min_poly_search += (T * 2)
+            min_poly_search += T
 
         sys.stderr.write("Need to sieve %d differences more...\nNew sieving range %d:%d\n" % (required_relations - len(Rels),start,stop))
       

From e75d0be0bc609d124a172ad97e5c45488610afc2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Mon, 14 Jun 2021 19:28:14 -0300
Subject: [PATCH 38/56] Update MPQS.sage

---
 MPQS.sage | 13 ++++---------
 1 file changed, 4 insertions(+), 9 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 2cad2e8..4f1cc19 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -504,17 +504,15 @@ class Poly:
             return NotImplemented
 
 
-def compute_logs_y_sums(start, stop, P, log_primes):
+def compute_logs_y_sums(l, P, log_primes):
      """
      Compute log table
      """
-     logs_y = [0] * ((stop - start) + 1)
+     logs_y = [0] * (l + 1)
      for j in range(len(P)):
          p = P[j]
          log_p = log_primes[j]
-         #for i in range(start + m, stop + m, p):
-         #    logs_y[i - m] += log_p
-         for i in range(start , stop , p):
+         for i in range(0, l, p):
              logs_y[i] += log_p
      return logs_y
 
@@ -553,10 +551,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y
         y = abs(y)
 
         if pre_log_filter:
-           #x = i
-           #print(x)
            candidate = (0 < logs_y[i] <= thresh and not is_prime(y))
-           a, b, c = is_square(y), is_prime(y), is_power_logprime(y, min_log_primes)
         else:
            candidate = (y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes))
 
@@ -881,7 +876,7 @@ def _MPQS(N, verbose=True, M = 2):
         filtered_log_primes = [log(p) for p in filtered_Prime_base]
       
         if (stop-start) > len(logs_y):
-            logs_y = compute_logs_y_sums(0, stop-start, filtered_Prime_base, filtered_log_primes)
+            logs_y = compute_logs_y_sums((stop-start), filtered_Prime_base, filtered_log_primes)
 
         required_relations = int(required_relations * required_relations_ratio)
         sys.stderr.write("Need %d relations\n" % (required_relations))

From c3d849e8c64b9ffb4009874fa3212ef858ca312f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Wed, 16 Jun 2021 16:05:35 -0300
Subject: [PATCH 39/56] better logp table generation

---
 MPQS.sage | 88 +++++++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 69 insertions(+), 19 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 4f1cc19..4426df0 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -504,19 +504,34 @@ class Poly:
             return NotImplemented
 
 
-def compute_logs_y_sums(l, P, log_primes):
+def compute_logs_y_sums(start, stop, P, log_primes, last_logp_table = []):
      """
      Compute log table
      """
-     logs_y = [0] * (l + 1)
+     logs_y = copy(last_logp_table)
+     logs_y += [0] * (stop-start)
      for j in range(len(P)):
          p = P[j]
          log_p = log_primes[j]
-         for i in range(0, l, p):
-             logs_y[i] += log_p
+         for i in range(0 , stop , p):
+             if i >= start:
+                 logs_y[i] += log_p
      return logs_y
 
 
+def logptable(start, stop, primes, logp, last_logp_table = []):
+    logp_table = copy(last_logp_table)
+    logp_table += [0] * (stop-start)
+    for j in range(0,len(primes)):
+        p = primes[j]
+        log_p = logp[j]
+        for i in range(0, stop, p):
+            if i >= start:
+                logp_table[i] += log_p
+    return logp_table
+
+
+
 def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
     """ 
     Relations search funcion 
@@ -551,7 +566,10 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y
         y = abs(y)
 
         if pre_log_filter:
+           #x = i
+           #print(x)
            candidate = (0 < logs_y[i] <= thresh and not is_prime(y))
+           a, b, c = is_square(y), is_prime(y), is_power_logprime(y, min_log_primes)
         else:
            candidate = (y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes))
 
@@ -802,6 +820,18 @@ def poly_stats(polys, polycounts):
             sys.stderr.write("For polynomial: %s the relations count is: %d\n" % (repr(poly),poly_count))
 
 
+def unique(List, newlist):
+    """
+    Creates a unique list of elements
+    """
+    for element in List:
+        if element not in tmp:
+            newlist.append(element)
+        else:
+            print("repeated")
+    return newlist
+    
+
 def _MPQS(N, verbose=True, M = 2):
     """ 
     Main MPQS function. 
@@ -823,23 +853,25 @@ def _MPQS(N, verbose=True, M = 2):
     T *= M
 
     B2, _ , B1 = prebuilt_params(log2(N))
-    Prime_base, log_p = find_primebase(N, B2)
-    #B1 //= M
-    B2 = len(Prime_base)
-  
-    if B2 == 1:
-        sys.stderr.write("Found small factor: %d\n" % Prime_base[0])
-        r, polycount = _MPQS(N // Prime_base[0])
-        return Prime_base + r, polycount
-
-    multiplier = choose_multiplier(N, Prime_base)
-    Nm = multiplier * N
     
     x_max = B1
     m_val = (x_max * root_2n) >> 1
     thresh = log10(m_val) * 0.735
     min_prime = int(thresh * 3)
 
+    Prime_base, log_p = find_primebase(N, B2 + min_prime)
+    
+    #print(B2,len(Prime_base), min_prime)
+    #sys.exit(0)
+
+    multiplier = choose_multiplier(N, Prime_base)
+    Nm = multiplier * N
+
+    if B2 == 1:
+        sys.stderr.write("Found small factor: %d\n" % Prime_base[0])
+        r, polycount = _MPQS(N // Prime_base[0])
+        return Prime_base + r, polycount
+
     sys.stderr.write("Factoring N: %d, bits: %d, digits: %d, B1: %d, B2: %d\n" % (N,bN,lN,B1,B2))
     sys.stderr.write("Multiplier is: %d\n" % multiplier)
 
@@ -859,7 +891,7 @@ def _MPQS(N, verbose=True, M = 2):
     polys = []
 
     last_lRels = 0
-    force_new_polys = False
+    force_new_polys = True
 
     required_relations_ratio = 1.05
     required_relations = len(Prime_base)
@@ -875,8 +907,8 @@ def _MPQS(N, verbose=True, M = 2):
         log_primes = [log(p) for p in Prime_base]
         filtered_log_primes = [log(p) for p in filtered_Prime_base]
       
-        if (stop-start) > len(logs_y):
-            logs_y = compute_logs_y_sums((stop-start), filtered_Prime_base, filtered_log_primes)
+        if stop > len(logs_y):
+            logs_y = compute_logs_y_sums(start, stop, filtered_Prime_base, filtered_log_primes, logs_y)
 
         required_relations = int(required_relations * required_relations_ratio)
         sys.stderr.write("Need %d relations\n" % (required_relations))
@@ -888,7 +920,9 @@ def _MPQS(N, verbose=True, M = 2):
         sys.stderr.write("Data collection with %d threads...\n" % T)
 
         lRels = len(Rels)
+
         need_more_polys = not (lRels > last_lRels) or force_new_polys
+
         last_lRels = lRels
 
         if need_more_polys == True: 
@@ -913,10 +947,18 @@ def _MPQS(N, verbose=True, M = 2):
         #polys += [Poly(N, Prime_base, x_max, search = None, verbose=False, A=1, B=0, C=-N)] # lower shield base poly for QS.
 
         # generate tasks parameters
-        for poly in polys[0 - (T * M):len(polys)]:
+        if not need_more_polys:
+            sys.stderr.write("Using the best polys to get more relations...\n")
+            polys_new = getbestpolys(polys, polycounts, T)
+        else:
+            polys_new = polys[0 - (T * M):len(polys)]
+
+        for poly in polys_new:
         #for poly in polys:
             inputs += [(taskid, Nm, start, stop, filtered_Prime_base, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
             taskid += 1
+            if poly not in polys:
+                poly.append(poly) 
 
         # deploy tasks to every cpu core.
         pols = []
@@ -928,13 +970,21 @@ def _MPQS(N, verbose=True, M = 2):
         
         t2 = time.time()
         sys.stderr.write("Done in: %f secs.\n" % (t2-t1))
+        #sys.stderr.write("Found %d rels with %d base primes.\nSorting..." % (len(Rels),len(Prime_base)))
         sys.stderr.write("Found %d rels with %d base primes.\n" % (len(Rels),len(Prime_base)))
 
+
         if len(cycleFactors) > 0:
             return cycleFactors, len(polys)
 
+        #unique(Rels, newlist)
+        #Rels = newlist
+        
         t3 = time.time()
 
+        #sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
+        #sys.stderr.write("Unique rels after sort: %d \n" % (len(Rels)))
+        
         # when needed relations is reached proceed to do linear algebra
         if len(Rels) > required_relations:
             sys.stderr.write("Found %d enough relations of %d needed relations...\n" % (len(Rels),required_relations))

From 0e2ac8249004c6c57e38bc7fdbbc677bf9fc8fdb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Wed, 16 Jun 2021 18:47:52 -0300
Subject: [PATCH 40/56] add multitasking

---
 PisanoPeriod.py | 219 ++++++++++++++++++++++++++----------------------
 1 file changed, 121 insertions(+), 98 deletions(-)

diff --git a/PisanoPeriod.py b/PisanoPeriod.py
index 3e4b867..2f69abc 100644
--- a/PisanoPeriod.py
+++ b/PisanoPeriod.py
@@ -7,20 +7,17 @@
 
 import random
 import time
-from gmpy2 import sqrt,isqrt, gcd, fib, fib2, f_mod as mod, powmod, is_prime, get_context, log2, log10
+from gmpy2 import *
 import sys
 sys.setrecursionlimit(5000)
+from multiprocessing import Pool, cpu_count, Manager
 
-ctx = get_context()
-ctx.precision += 1000
-sqrt5 = sqrt(5)
-
-import bitarray
 
 class Fibonacci:
-    def __init__(self):
-        pass
-    
+    def __init__(self, progress = False, verbose = True, multitasked = True):
+        self.progress = progress
+        self.verbose = verbose
+        self.multitasked = multitasked
 
     def _fib_res(self,n,p):
         """ fibonacci sequence nth item modulo p """
@@ -34,57 +31,37 @@ def _fib_res(self,n,p):
                 return (c, d)
             else:
                 return (d, mod((c + d), p))
-    
-    def _fib_eig(self,n,p):
-        lambda1 = (1 + sqrt5) / 2
-        lambda2 = (1 - sqrt5) / 2
-        fk = (pow(lambda1, n) - pow(lambda2, n)) / (lambda1 - lambda2)
-        return mod(int(fk), p)
 
 
-    def get_n_mod_d(self,n,d, use='mersenne'):
-    #def get_n_mod_d(self,n,d, use=''):
+    def get_n_mod_d(self,n,d, use = 'mersenne'):
         if n < 0:
             ValueError("Negative arguments not implemented")
         if use == 'gmpy':
             return mod(fib(n), d)
-        elif use == 'eig':
-            return self._fib_eig(n,d) # takes forever
         elif use == 'mersenne':
-            return mod(powmod(2,n, d)-1, d)
+            return powmod(2, n, d) - 1
         else:
             return self._fib_res(n,d)[0]
 
+    def ranged_period(self, N, start, stop, look_up_dest):
+        print("ranged_period (%d,%d) start" % (start,stop))
+        tmp_look_up = {}
+        for x in range(start, stop):
+            tmp_look_up[self.get_n_mod_d(x, N)] = x
+        look_up_dest.update(tmp_look_up)
+        #look_up_dest |= tmp_look_up
+        print("ranged_period (%d,%d) end" % (start,stop))
+        return 1
     
-    def binary_search(self,L,n):
-        """ Finds item index in O(log2(N)) """ 
-        left = 0
-        right = len(L) - 1
-        while left <= right:
-            mid=(left + right) >> 1
-            if n == L[mid]:
-                return mid
-            elif n < L[mid]:
-                right = mid - 1
-            else:
-                left = mid + 1
-        return -1
-    
-    
-    def sort_list(self,L):
-        from operator import itemgetter
-        indices, L_sorted = zip(*sorted(enumerate(L), key=itemgetter(1)))
-        return list(L_sorted),list(indices)
-   
- 
-    def get_period_bigint(self, N, min_accept, xdiff, verbose = False, use_qbf = False):            
+
+    def get_period_bigint(self, N, min_accept, xdiff, verbose = False):            
         search_len = int(pow(N, (1.0 / 6) / 100))
         
         if search_len < min_accept:
             search_len = min_accept
   
-        if verbose:
-            print('search_len:%d'%(search_len))
+        if self.verbose:
+            print('Search_len: %d, log2(N): %d' % (search_len,int(log2(N))))
         
         starttime = time.time()
         diff = xdiff 
@@ -94,61 +71,72 @@ def get_period_bigint(self, N, min_accept, xdiff, verbose = False, use_qbf = Fal
             begin = 1
         end = N + int('9' * p_len)
     
-        if verbose:    
-            print('search begin:%d,end:%d'%(begin, end))
-               
-        rs = [self.get_n_mod_d(x, N) for x in range(search_len)]
-        rs_sort, rs_indices = self.sort_list(rs)
-
-        if use_qbf:
-            QBF = bitarray.bitarray(rs_sort[-1]+1) 
-            QBF.setall(0)
-
-            for index in rs_sort:
-                QBF[index] = 1
-
-        if verbose:    
-            #print(rs, rs_sort, rs_indices)        
-            print('sort complete! time used: %f secs' % (time.time() - starttime))
-                
-        T = 0
-        has_checked_list = []
+        if self.verbose:    
+            print('Search begin: %d, end: %d'%(begin, end))
+        
+        if self.multitasked:
+            C = cpu_count() * 2
+            search_work = search_len // C
+
+            manager =  Manager()
+            look_up = manager.dict()
+
+            if self.verbose:
+                print("Precompute LUT with %d tasks..." % C)
+
+            inputs = []
+            for x in range(0, search_len, search_work):
+                inputs += [(N, x, x + search_work, look_up)]
+
+            workpool = Pool(C)
+
+            with workpool:
+                results = workpool.starmap(self.ranged_period,inputs)
+                print(results)
+                workpool.close()
+                workpool.join()
+
+        else:
+            look_up = {}
+            self.ranged_period(N, 0, search_len, look_up)
+
+        if self.verbose:
+            print("LUT creation ended size: %d..." % len(look_up))
+            print("Searching...")
+        
 
         while True:       
             randi = random.randint(begin,end)            
-            if (not randi in has_checked_list):
-                has_checked_list.append(randi)
-            
-                res = self.get_n_mod_d(randi, N)
-             
-                if use_qbf: 
-                    if res < len(QBF) and QBF[res] == 1:
-                        inx = self.binary_search(rs_sort, res)
-                        print("res:",res,"inx:",inx)
-                    else:
-                        inx = -1
-                else:
-                    inx = self.binary_search(rs_sort, res)
-                                    
-                if inx > -1:                
-                    res_n = rs_indices[inx]
-                    T = randi - res_n
-                     
-                    if self.get_n_mod_d(T, N) == 0:
-                        td = int(time.time() - starttime)
-                        if verbose:
-                            print('For N = %d Found T:%d, randi: %d, time used %f secs.' % (N , T, randi, td))
-                        return td, T, randi
-                    else:
-                        if verbose:
-                            print('For N = %d\n Found res: %d, inx: %d, res_n: %d , T: %d\n but failed!' % (N, res, inx, res_n, T))
+            res = self.get_n_mod_d(randi, N)
+            if res > 0:
+                #print(res, res in look_up)
+                if res in look_up:
+                    res_n = look_up[res]
+                    T = randi - res_n   
+                    if T & 1 == 0: 
+                        if self.get_n_mod_d(T, N) == 0:
+                            td = int(time.time() - starttime)
+                            if self.verbose:
+                                print('For N = %d Found T:%d, randi: %d, time used %f secs.' % (N , T, randi, td))
+                            return td, T, randi
+                        else:
+                            if self.verbose:
+                                print('For N = %d\n Found res: %d, res_n: %d , T: %d\n but failed!' % (N, res, res_n, T))
+            else:
+                if randi & 1 == 0:
+                    T = randi
+                    td = int(time.time() - starttime)
+                    if self.verbose:
+                        print('First shot, For N = %d Found T:%d, randi: %d, time used %f secs.' % (N , T, randi, td))
+                    return td, T, randi
 
   
-    def _trivial_factorization_with_n_phi(self, N, phi):
+    def _trivial_factorization_with_n_phi(self, N, T):
+        phi = abs(N - T) + 1
         p1 = []
         d2 = N << 2
 
-        phi2 = pow(phi,2)
+        phi2 = pow(phi, 2)
         phi2p4d = phi2 + d2
         phi2m4d = phi2 - d2
 
@@ -181,8 +169,7 @@ def _trivial_factorization_with_n_phi(self, N, phi):
     def factorization(self, N, min_accept, xdiff):
         res = self.get_period_bigint(N, min_accept, xdiff, verbose=True) 
         if res != None:
-            t, T, r = res
-            phi = abs(N - T) + 1
+            t, phi, r = res 
             return self._trivial_factorization_with_n_phi(N, phi)
 
 """
@@ -228,9 +215,11 @@ def test(Ns,B2=0):
     n = 1
     tti = time.time()
     for N in Ns:
+        if int(log2(N)) < 38:
+            continue
         ti = time.time()
-        B1, B2 = pow(10, int((log10(N)) // 2)-1), B2
-        print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,log2(N),B1,B2))
+        B1, B2 = pow(10, int((log10(N)) // 2)-0), B2
+        print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,int(log2(N)),B1,B2))
         P = Fib.factorization(N,B1,B2)
         if P != None:
             phi = (P[0]-1) * (P[1]-1)
@@ -255,7 +244,7 @@ def test2():
   for B1 in B1s:
       for B2 in B2s:
           ti = time.time()
-          print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,log2(N),B1,B2))
+          print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,int(log2(N)),B1,B2))
           P = Fib.factorization(N,B1,B2)
           if P != None:
               phi = (P[0]-1) * (P[1]-1)
@@ -273,14 +262,48 @@ def test3(N, B2 = 0):
   Fib = Fibonacci()
   ti = time.time()
   B1, B2 = pow(10, int((log10(N)) // 2)), B2
-  print("Test: N: %d, log2(N): %d, B1: %d, B2: %d" % (N,log2(N),B1,B2))
+  print("Test: N: %d, log2(N): %d, B1: %d, B2: %d" % (N,int(log2(N)),B1,B2))
   P = Fib.factorization(N,B1,B2)
   print(P)
   td = time.time() - ti
   print("Runtime: %f" % (td))
 
+def test4(l,B2=0):
+    Fib = Fibonacci(multitasked = True)
+    #times = []
+    ff = 0
+    n = 1
+    tti = time.time()
+    for n in range(10,l):
+        ti = time.time()
+        n2 = 2**(n//2)
+        n21 = (2**(n//2)-1)
+        a = random.randint(n21,n2)
+        p = next_prime(a)
+        q = p
+        while q == p:
+           b = random.randint(n21,n2)
+           q = next_prime(next_prime(b))
+        N = p * q
+        B1, B2 = pow(10, int((log10(N)) // 2)-0), B2
+        print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,int(log2(N)),B1,B2))
+        P = Fib.factorization(N,B1,B2)
+        if P != None:
+            phi = (P[0]-1) * (P[1]-1)
+            print("phi(N): %d" % phi)
+            print("factors: %s" % str(P))
+            ff += 1
+        td = time.time() - ti
+        ttd = time.time() - tti
+        print("Runtime: %f\nFully factored:%d of %d" % (td,ff,l))
+        print("Total Runtime: %f" % (ttd))
+        n += 1
+        print('------------------------------------------------------------------')
+
+
 
 if __name__=='__main__':
-   test(Ns0+Ns1+Ns2,B2=int(sys.argv[1]))
+   #test(Ns0+Ns1+Ns2,B2=int(sys.argv[1]))
    #test2()
    #test3(int(sys.argv[1]))
+   test4(100)

From 2be9010fb98a95f450fb19991b7656110810f95e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Wed, 16 Jun 2021 18:52:49 -0300
Subject: [PATCH 41/56] Update PisanoPeriod.py

---
 PisanoPeriod.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PisanoPeriod.py b/PisanoPeriod.py
index 2f69abc..d6631a7 100644
--- a/PisanoPeriod.py
+++ b/PisanoPeriod.py
@@ -74,7 +74,7 @@ def get_period_bigint(self, N, min_accept, xdiff, verbose = False):
         if self.verbose:    
             print('Search begin: %d, end: %d'%(begin, end))
         
-        if self.multitasked:
+        if self.multitasked and search_len > 1000:
             C = cpu_count() * 2
             search_work = search_len // C
 

From acc3bca1a985914eefe4350859466e6ea8f34020 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sat, 19 Jun 2021 02:25:51 -0300
Subject: [PATCH 42/56] Working 99 bit

---
 MPQS.sage | 248 ++++++++++++++++++++++++++++--------------------------
 1 file changed, 128 insertions(+), 120 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 4426df0..96a86af 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -365,9 +365,6 @@ class Poly:
         s=0
         leg = 2
         while True:
-            root_A = next_prime(root_A + s)
-            #root_A = next_prime(root_A)
-
             leg = legendre(n, root_A)
             if leg == 1:
                 if s > self.search:
@@ -378,7 +375,7 @@ class Poly:
                     break
                 #break
             s += 1
-
+            root_A = next_prime(root_A + s)
         self.s = s
         A = int(pow(root_A, 2))
         # solve for an adequate B. B*B is a quadratic residue mod n,
@@ -532,102 +529,104 @@ def logptable(start, stop, primes, logp, last_logp_table = []):
 
 
 
-def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly = None):
+def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, polys = None):
     """ 
     Relations search funcion 
     """
     #pid = os.getpid()
     pid = taskid
     tasks.value += 1
-    sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d sieving start, tresh: %f.\n" % (pid,start, stop, (stop-start), thresh))
-    m = poly.minima
-    Diffs = [poly.eval(abs(x)) for x in range(start + m, stop +m)]
-    
-
-    ld = len(Diffs)
+    sys.stderr.write("[%d] relations_find: range(%d, %d), interval: %d sieving start, tresh: %f, pb: %d.\n" % (pid,start, stop, (stop-start), thresh,len(polys)))
     pre_log_filter = True
-
+    merge = True
     proc = noproc = 0
     Found_Rels = []
-    A = poly.A
     st = time.time()
     ltd = st
-
     m = 1000
     msg = ""
     partials = {}
-    i= 0
     rels_found = 0
-    for i in range(ld):
-        if len(Rels) > required_relations or len(cycleFactors) > 0:
-            polycounts[poly] = rels_found
-            break
-        y, Rad, x = Diffs[i]
-        y = abs(y)
-
-        if pre_log_filter:
-           #x = i
-           #print(x)
-           candidate = (0 < logs_y[i] <= thresh and not is_prime(y))
-           a, b, c = is_square(y), is_prime(y), is_power_logprime(y, min_log_primes)
-        else:
-           candidate = (y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes))
-
-        if candidate:
-            proc += 1
-            f = minifactor4(y, P, smooth_base)
-            if f != None:            
-                if f[1] == 1:  # found a relation 
-                    rels_found += 1
-                    filtered = filter_out_even_powers(f[0])
-                    rel = [filtered, y, Rad, A] 
-                    Rels.append(rel)
-                elif f[1] in partials: # found a partial and try to merge
-                    a = partials[f[1]]
-                    p = filter_out_even_powers(f[0] + a[0])
-                    Ahs = A*a[3]
-                    lhs = Rad * a[2]
-                    rhs = y * a[1]
-                    LHS = Ahs + lhs
-                    g = gcd(isqrt(LHS) - rhs, N)
-                    if N > g > 1:
-                        cycleFactors.append([g , N // g])
-                        sys.stderr.write("Found cycle with partial\n")
-                    else: # didnt merge so store the partial for later merging
-                        Rels.append([p, rhs, lhs, Ahs])
+
+    for poly in polys:
+        m = poly.minima
+        Diffs = [poly.eval(abs(x)) for x in range(start + m, stop +m)]
+        ld = len(Diffs)
+        A = poly.A
+        for i in range(ld):
+            if len(Rels) > required_relations or len(cycleFactors) > 0:
+                polycounts[poly] = (start,stop,rels_found)
+                break
+            y, Rad, x = Diffs[i]
+            y = abs(y)
+
+            if pre_log_filter:
+                candidate = (0 < logs_y[i] <= thresh and not is_prime(y))
+            else:
+                candidate = (y > P[0] and not is_prime(y) and not is_square(y) and not is_power_logprime(y, min_log_primes))
+
+            if candidate:
+                proc += 1
+                f = minifactor4(y, P, smooth_base)
+                if f != None:            
+                    if f[1] == 1:  # found a relation 
                         rels_found += 1
-                    #with merged_count.value.get_lock():
-                    merged_count.value += 1
-                else:
-                    partials[f[1]] = [f[0], y, Rad, A]
-        else:
-            noproc += 1
-
-        if i % m == 0:
-            lRels = len(Rels)
-            lt = time.time()
-            td = lt - ltd
-            ltd = lt
-            if i > 0:
-                if lRels > 0:
-                    eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
-                else:
-                    eta = int(td * (ld / m))
-                etas = humanfriendly.format_timespan(eta)
+                        filtered = filter_out_even_powers(f[0])
+                        rel = [filtered, y, Rad, A] 
+                        Rels.append(rel)
+                    elif merge and (f[1] in partials): # found a partial and try to merge
+                        a = partials[f[1]]
+                        p = filter_out_even_powers(f[0] + a[0])
+                        Ahs = A*a[3]
+                        lhs = Rad * a[2]
+                        rhs = y * a[1]
+                        LHS = Ahs + lhs
+                        g = gcd(isqrt(LHS) - rhs, N)
+                        if N > g > 1:
+                            cycleFactors.append([g , N // g])
+                            sys.stderr.write("Found cycle with partial\n")
+                        else: # didnt merge so store the partial for later merging
+                            if merge:
+                                rel = [p, rhs, lhs, Ahs]
+                                if rel not in Rels:
+                                    Rels.append(rel)
+                                    rels_found += 1
+                                    #with merged_count.value.get_lock():
+                                    merged_count.value += 1
+                                    del partials[f[1]]
+                    else:
+                        if merge:
+                            partials[f[1]] = [f[0], y, Rad, A]
             else:
-                etas = "No ETA estimated"
-            tds = humanfriendly.format_timespan(td)
-            msg = "[%d] relations_find: range (%d, %d , %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, iter_elapsed: %s, eta: %s, alive tasks: %d.\n" 
-            msg = msg % (pid, start, i, ld, lRels, required_relations, merged_count.value, proc, noproc , tds, etas, tasks.value)
-            sys.stderr.write(msg)
+                noproc += 1
+
+            if i % m == 0:
+                lRels = len(Rels)
+                lt = time.time()
+                td = lt - ltd
+                ltd = lt
+                if i > 0:
+                    if lRels > 0:
+                        eta = int(td * sqrt((ld / m)**2 + (required_relations / lRels)**2))
+                    else:
+                        eta = int(td * (ld / m))
+                    etas = humanfriendly.format_timespan(eta)
+                else:
+                    etas = "No ETA estimated"
+                tds = humanfriendly.format_timespan(td)
+                msg = "[%d] relations_find: range (%d, %d , %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, iter_elapsed: %s, eta: %s, alive tasks: %d.\n" 
+                msg = msg % (pid, start, i, ld, lRels, required_relations, merged_count.value, proc, noproc , tds, etas, tasks.value)
+                sys.stderr.write(msg)
     
-    polycounts[poly] = rels_found
+        polycounts[poly] = (start,stop,rels_found)
+
     td = time.time() - st
     tds = humanfriendly.format_timespan(td)
     msg = "[%d] relations_find: Ended range(%d, %d, %d), found: %d of %d, merged: %d, proc: %d, noproc: %d, time elapsed: %s, alive tasks: %d.\n" 
     msg = msg % (pid, start,i,ld,len(Rels),required_relations, merged_count.value, proc, noproc, tds, tasks.value)
     sys.stderr.write(msg)
     tasks.value -= 1
+
     return Found_Rels
     
 
@@ -736,6 +735,7 @@ def find_primebase(n, bound):
     https://github.com/elliptic-shiho/primefac-fork/blob/master/_primefac/_factor_algo/_mpqs.py
     Same as the sieve of erathostenes.
     """
+    print(n,bound)
     primes, mod_root, log_p, num_prime, p = [], {}, {}, 0, 3
     while num_prime <= bound:
         leg = legendre(n % p, p)
@@ -795,18 +795,27 @@ def getbestpolys(polys, polycounts, T):
     new_polys = []
     print("have polys: %d" % len(tmp_polys))
     print(polys)
+    j = -1
     while len(new_polys) < T:
         max_poly_count = 0
         for i in range(0,len(tmp_polys)):
             poly = tmp_polys[i]
             if poly in polycounts:
-                poly_count = polycounts[poly]
+                poly_count = polycounts[poly][2]
                 if poly_count > max_poly_count:
                     max_poly_count = poly_count
-                    best_poly = poly
-                    del polys[i]
-        print("new best poly: %s" % repr(poly))
-        new_polys.append(poly)
+                    best_poly = poly                
+                    j = i
+
+        #if -1 < j < len(tmp_polys):
+        if -1 < j < len(tmp_polys):
+            print("new best poly: %s, j: %d, len: %d, count: %d" % (repr(poly),j,len(tmp_polys),max_poly_count))
+            new_polys.append(best_poly)
+            del tmp_polys[j]
+        else:
+            print(tmp_polys)
+            break
+    sys.stderr.write("Requested best_polys: %d, got: %d\n" % (T,len(new_polys)))
     return new_polys
 
 
@@ -817,7 +826,7 @@ def poly_stats(polys, polycounts):
     for poly in polys:
         if poly in polycounts:
             poly_count = polycounts[poly]
-            sys.stderr.write("For polynomial: %s the relations count is: %d\n" % (repr(poly),poly_count))
+            sys.stderr.write("For polynomial: %s the relations count is: %s\n" % (repr(poly),poly_count))
 
 
 def unique(List, newlist):
@@ -832,7 +841,7 @@ def unique(List, newlist):
     return newlist
     
 
-def _MPQS(N, verbose=True, M = 2):
+def _MPQS(N, verbose=True, M = 5):
     """ 
     Main MPQS function. 
     """
@@ -849,8 +858,8 @@ def _MPQS(N, verbose=True, M = 2):
     root_2n = isqrt(2*N)
     Rels = []
 
-    T = cpu_count()
-    T *= M
+    T = cpu_count() * M
+    #T *= M
 
     B2, _ , B1 = prebuilt_params(log2(N))
     
@@ -860,14 +869,13 @@ def _MPQS(N, verbose=True, M = 2):
     min_prime = int(thresh * 3)
 
     Prime_base, log_p = find_primebase(N, B2 + min_prime)
-    
-    #print(B2,len(Prime_base), min_prime)
-    #sys.exit(0)
 
     multiplier = choose_multiplier(N, Prime_base)
     Nm = multiplier * N
+    
+    required_relations = len(Prime_base)
 
-    if B2 == 1:
+    if required_relations == 1:
         sys.stderr.write("Found small factor: %d\n" % Prime_base[0])
         r, polycount = _MPQS(N // Prime_base[0])
         return Prime_base + r, polycount
@@ -892,9 +900,10 @@ def _MPQS(N, verbose=True, M = 2):
 
     last_lRels = 0
     force_new_polys = True
+    use_best = True
 
     required_relations_ratio = 1.05
-    required_relations = len(Prime_base)
+    required_relations = int(required_relations * required_relations_ratio)
 
     logs_y = []
 
@@ -909,8 +918,8 @@ def _MPQS(N, verbose=True, M = 2):
       
         if stop > len(logs_y):
             logs_y = compute_logs_y_sums(start, stop, filtered_Prime_base, filtered_log_primes, logs_y)
+            required_relations = int(required_relations * required_relations_ratio)
 
-        required_relations = int(required_relations * required_relations_ratio)
         sys.stderr.write("Need %d relations\n" % (required_relations))
         
         smooth_base = prod(Prime_base)
@@ -922,43 +931,46 @@ def _MPQS(N, verbose=True, M = 2):
         lRels = len(Rels)
 
         need_more_polys = not (lRels > last_lRels) or force_new_polys
+        need_more_polys = False
 
         last_lRels = lRels
 
-        if need_more_polys == True: 
-            sys.stderr.write("Need more polys...\n") 
-            polys, early_factors = generate_polys(Nm, Prime_base, x_max, T, min_poly_search, polys) # generate n distinct polys one for each cpu core.
-
-            if len(early_factors) > 0:
-                tmp = 1
-                small = []
-                for early_factor in early_factors:
-                    g = gcd(early_factor, N)
-                    if N > g > 1:
-                        tmp *= g
-                        small.append(g)
-                if tmp > 1:
-                    return small + MPQS(N // tmp), 0
-             
-        logs_y = compute_logs_y_sums(0, (stop-start), filtered_Prime_base, filtered_log_primes)
-
         inputs = [] 
         taskid = 1
-        #polys += [Poly(N, Prime_base, x_max, search = None, verbose=False, A=1, B=0, C=-N)] # lower shield base poly for QS.
 
         # generate tasks parameters
         if not need_more_polys:
-            sys.stderr.write("Using the best polys to get more relations...\n")
-            polys_new = getbestpolys(polys, polycounts, T)
+            if len(polys) > 0 and use_best:
+                sys.stderr.write("Using the best polys to get more relations...\n")
+                polys_new = getbestpolys(polys, polycounts, T)
+                polys_new += polys[0 - T:len(polys)]
+            else:
+                #polys_new = generate_polys(Nm, Prime_base, x_max, T, min_poly_search, polys)
+                sys.stderr.write("Need more polys...\n")
+                polys_new, early_factors = generate_polys(Nm, Prime_base, x_max, T, min_poly_search, polys) # generate n distinct polys one for each cpu core.
+                if len(early_factors) > 0:
+                    tmp = 1
+                    small = []
+                    for early_factor in early_factors:
+                        g = gcd(early_factor, N)
+                        if N > g > 1:
+                            tmp *= g
+                            small.append(g)
+                    if tmp > 1:
+                        return small + MPQS(N // tmp), 0
+                min_poly_search += T
         else:
             polys_new = polys[0 - (T * M):len(polys)]
+            
+        pb = M
+        for i in range(0, len(polys_new), pb):
+            polys_batch = polys_new[i : i + pb]
 
-        for poly in polys_new:
-        #for poly in polys:
-            inputs += [(taskid, Nm, start, stop, filtered_Prime_base, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, poly)]
+            inputs += [(taskid, Nm, start, stop, filtered_Prime_base, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, polys_batch)]
             taskid += 1
-            if poly not in polys:
-                poly.append(poly) 
+            for poly in polys_batch:
+                if poly not in polys:
+                    poly.append(poly) 
 
         # deploy tasks to every cpu core.
         pols = []
@@ -977,9 +989,6 @@ def _MPQS(N, verbose=True, M = 2):
         if len(cycleFactors) > 0:
             return cycleFactors, len(polys)
 
-        #unique(Rels, newlist)
-        #Rels = newlist
-        
         t3 = time.time()
 
         #sys.stderr.write("Done in: %f secs.\n" % (t3-t2))
@@ -1017,7 +1026,6 @@ def MPQS(N):
     """ 
     Iterative version of MPQS. 
     """
-
     result, polys = _MPQS(N)
     polycount = polys
 

From f0b207046f700b536662025b879bd6e7f218e6c3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sat, 19 Jun 2021 02:27:56 -0300
Subject: [PATCH 43/56] get rid of duplicated funtion

---
 MPQS.sage | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 96a86af..60fe3ed 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -516,19 +516,6 @@ def compute_logs_y_sums(start, stop, P, log_primes, last_logp_table = []):
      return logs_y
 
 
-def logptable(start, stop, primes, logp, last_logp_table = []):
-    logp_table = copy(last_logp_table)
-    logp_table += [0] * (stop-start)
-    for j in range(0,len(primes)):
-        p = primes[j]
-        log_p = logp[j]
-        for i in range(0, stop, p):
-            if i >= start:
-                logp_table[i] += log_p
-    return logp_table
-
-
-
 def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y, smooth_base, Rels, merged_count, required_relations, cycleFactors, thresh, tasks, polycounts, polys = None):
     """ 
     Relations search funcion 

From 55718a17d051373248a2f03a8791932ae303ae1c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Sat, 19 Jun 2021 02:34:08 -0300
Subject: [PATCH 44/56] Update MPQS.sage

---
 MPQS.sage | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index 60fe3ed..db70e6b 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -556,15 +556,16 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y
                 proc += 1
                 f = minifactor4(y, P, smooth_base)
                 if f != None:            
-                    if f[1] == 1:  # found a relation 
+                    if f[1] == 1:  # Found a relation 
                         rels_found += 1
                         filtered = filter_out_even_powers(f[0])
-                        rel = [filtered, y, Rad, A] 
-                        Rels.append(rel)
-                    elif merge and (f[1] in partials): # found a partial and try to merge
+                        rel = [filtered, y, Rad, A]
+                        if rel not in Rels:
+                            Rels.append(rel)
+                    elif merge and (f[1] in partials): # Found a partial, try to find a cycle
                         a = partials[f[1]]
                         p = filter_out_even_powers(f[0] + a[0])
-                        Ahs = A*a[3]
+                        Ahs = A * a[3]
                         lhs = Rad * a[2]
                         rhs = y * a[1]
                         LHS = Ahs + lhs
@@ -572,7 +573,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y
                         if N > g > 1:
                             cycleFactors.append([g , N // g])
                             sys.stderr.write("Found cycle with partial\n")
-                        else: # didnt merge so store the partial for later merging
+                        else: # Cycle not found, merge
                             if merge:
                                 rel = [p, rhs, lhs, Ahs]
                                 if rel not in Rels:
@@ -582,7 +583,7 @@ def relations_find(taskid, N, start, stop, P, min_log_primes, log_primes, logs_y
                                     merged_count.value += 1
                                     del partials[f[1]]
                     else:
-                        if merge:
+                        if merge: # Didnt merge, store the partial for later merging
                             partials[f[1]] = [f[0], y, Rad, A]
             else:
                 noproc += 1

From 0c48a2a72ad8bfaf59f9a76ec3a9ead30707e5f5 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Mon, 19 Dec 2022 13:53:45 -0300
Subject: [PATCH 45/56] msieve_driver

---
 msieve_driver.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)
 create mode 100644 msieve_driver.py

diff --git a/msieve_driver.py b/msieve_driver.py
new file mode 100644
index 0000000..0fc41bb
--- /dev/null
+++ b/msieve_driver.py
@@ -0,0 +1,20 @@
+
+def msieve_factor(n):
+  import subprocess, re
+  BIN = "/home/dclavijo/code/FACTORING/msieve/msieve"
+  #ARGS="-v %d" % n
+  tmp = []
+  proc = subprocess.Popen([BIN,"-t","8","-v",str(n)],stdout=subprocess.PIPE)
+  for line in proc.stdout:
+    #the real code does filtering here
+    #print("test:", line.rstrip())
+    line = line.rstrip().decode("utf8")
+    #print(line)
+    if re.search("factor: ",line):
+      tmp += [int(line.split()[2])]
+  return tmp
+
+#n = 196616705776852626592749272255227899
+import sys
+n = int(sys.argv[1])
+print(msieve_factor(n))

From 2cd0915e2ffa08b95b8b806985a83eb54199ec8b Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Mon, 19 Dec 2022 13:54:18 -0300
Subject: [PATCH 46/56] Somsuk algorithm

---
 somsuk.py | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 61 insertions(+)
 create mode 100644 somsuk.py

diff --git a/somsuk.py b/somsuk.py
new file mode 100644
index 0000000..97bf194
--- /dev/null
+++ b/somsuk.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+# Author Dario Clavijo 2022
+
+from gmpy2 import isqrt, gcd, sqrt, powmod, isqrt_rem, is_square
+import random
+import time
+
+def fermat(n):
+    a, rem = isqrt_rem(n)
+    b2 = -rem
+    c0 = (a << 1) + 1
+    c = c0
+    while not is_square(b2):
+        b2 += c
+        c += 2
+    a = (c - 1) >> 1
+    b = isqrt(b2)
+    return a - b, a + b
+
+
+def somsuk(n):
+  """
+  Implementation based on Kritsanapong Somsuk paper: 
+  The new integer factorization algorithm based on fermat’s
+  factorization algorithm and euler’s theorem (2020).
+  ISSN: 2088-8708, DOI: 10.11591/ijece.v10i2.pp1469-1476
+  """
+  u = isqrt(n) << 1
+  c = n
+  while gcd(c, n) > 1:
+    c = random.randrange(2, n - 1)
+  a = powmod(c, -1, n)
+  s = powmod(c, 2, n)
+  phi_aprox = n - u  + 1
+  t = powmod(a, phi_aprox, n)
+  while True:
+    if t == 1:
+      x = u >> 1
+      y = isqrt(x * x - n)
+      return x - y, x + y
+    t = (t*s) % n
+    u += 2
+
+def timeit(f,n):
+  t0 = time.time()
+  r = f(n)
+  t1 = time.time()
+  print(str(f).split(" ")[1],"(", n , ") =", r,"time:",round(t1-t0,5))
+
+def test():
+  timeit(somsuk, 340213)
+  timeit(fermat, 340213)
+  timeit(somsuk, 788582867650121563)
+  timeit(fermat, 788582867650121563)
+  timeit(somsuk, 1047329636821139813)
+  timeit(fermat, 1047329636821139813)
+
+if __name__ == "__main__":
+  test()
+
+

From edd6238916987cd79459bbf1d3fd4ef8b821c156 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Wed, 21 Dec 2022 17:37:02 -0300
Subject: [PATCH 47/56] update qs.sage

---
 qs.sage | 76 ++++++++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 62 insertions(+), 14 deletions(-)

diff --git a/qs.sage b/qs.sage
index efc1ab6..af46b8d 100644
--- a/qs.sage
+++ b/qs.sage
@@ -6,6 +6,56 @@ from gmpy2 import gcd, isqrt, is_prime, next_prime
 import itertools
 
 
+def prebuilt_params(bits):
+    """
+    Bounds estimation
+    borrowed from msieve/mpqs.c
+    """
+    if bits <= 64:
+        return [100, 40, 1 * 65536]
+    if bits <= 128: 
+        return [450, 40, 1 * 65536]
+    if bits <= 183:
+        return [2000, 40, 1 * 65536]
+    if bits <= 200: 
+        return [3000, 50, 1 * 65536]
+    if bits <= 212: 
+        return [5400, 50, 3 * 65536]
+    if bits <= 233:
+        return [10000, 100, 3 * 65536]
+    if bits <= 249:
+        return [27000, 100, 3 * 65536]
+    if bits <= 266:
+        return [50000, 100, 3 * 65536]
+    if bits <= 283:
+        return [55000, 80, 3 * 65536]
+    if bits <= 298:
+        return [60000, 80, 9 * 65536]
+    if bits <= 315:
+        return [80000, 150, 9 * 65536]
+    if bits <= 332:
+        return [100000, 150, 9 * 65536]
+    if bits <= 348:
+        return [140000, 150, 9 * 65536]
+    if bits <= 363:
+        return [210000, 150, 13 * 65536]
+    if bits <= 379:
+        return [300000, 150, 17 * 65536]
+    if bits <= 395:
+        return [400000, 150, 21 * 65536]
+    if bits <= 415:
+        return [500000, 150, 25 * 65536] # beyond this point you're crazy 
+    if bits <= 440:
+        return [700000, 150, 33 * 65536]
+    if bits <= 465:
+        return [900000, 150, 50 * 65536]
+    if bits <= 490:
+        return [1100000, 150, 75 * 65536]
+    if bits <= 512:
+        return [1300000, 150, 100 * 65536]
+    return [1300000, 150, 100 * 65536]
+
+
 def trial_division_minus_even_powers(n,P):
     a = []
     r = n
@@ -33,20 +83,18 @@ def trial_division_minus_even_powers(n,P):
     return [a,r,n]
 
 
-def minifactor(x,P):
-    if not is_prime(x):
-        i = isqrt(x)
-        if pow(i,2) != x: # perfect squares are even powers
-            p = trial_division_minus_even_powers(x,P)
-            if p[1] == 1: # if 1 x is B-smooth
-                return p
+def minifactor(x, P):
+    if not is_prime(x) and not is_square(x):
+        p = trial_division_minus_even_powers(x, P)
+        if p[1] == 1: # if 1 x is B-smooth
+            return p
 
 
 def QS(N,differences,y):
     i2N = isqrt(N)
-    X = [Integer(a^2 - N) for a in range(i2N + 1,i2N + differences)]
+    X = [int(a * a - N) for a in range(i2N + 1, i2N + differences)]
     P = list(primes(2, y))
-    F = list(filter(None,map(minifactor, X, itertools.repeat(P, len(X)))))
+    F = list(filter(None, map(minifactor, X, itertools.repeat(P, len(X)))))
 
     M = matrix(GF(2), len(F), len(P), lambda i, j:P[j] in F[i][0])
     for K in M.left_kernel().basis():
@@ -54,14 +102,14 @@ def QS(N,differences,y):
         x = product([int(isqrt(f[2] + N)) for f in I])
         y = isqrt(product([f[2] for f in I]))
         g0, g1 = gcd(N, x - y), gcd(N, x + y)
-        if N > g0 > 1:
-            return g0, N // g0
-        if N > g1 > 1:
-            return g1, N // g1
+        if N > g0 > 1 or N > g1 > 1:
+            return g0, g1
 
 
 if __name__ == "__main__":
-    r = QS(Integer(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3]))
+    N = int(sys.argv[1])
+    y ,_ ,diff = prebuilt_params(N.bit_length())
+    r = QS(N, diff * 10, y * 10)
     if r != None:
         print("%d, %d" % r)
 

From c583c58eaa28c6217b5df6eafd4f430a50404cb7 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Thu, 5 Jan 2023 16:30:00 -0300
Subject: [PATCH 48/56] add lehman algorithm

---
 lehman.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100644 lehman.py

diff --git a/lehman.py b/lehman.py
new file mode 100644
index 0000000..07359d6
--- /dev/null
+++ b/lehman.py
@@ -0,0 +1,43 @@
+from gmpy2 import *
+
+
+def cuberoot(n):
+  a,b = iroot(n,3)
+  return a
+
+
+def lehman(n):
+    """
+    based on: https://programmingpraxis.com/2017/08/22/lehmans-factoring-algorithm/
+    """
+    if is_congruent(n, 2, 4):
+        return []
+
+    for k in range(1, cuberoot(n), 1):
+        nk4 = n*k << 2
+        ki4 = isqrt(k) << 2
+        ink4 = isqrt(nk4) + 1
+        i6,_ = iroot(n, 6)
+        ink4i6ki4 = ink4 + (i6 // (ki4))
+        for a in range(ink4, ink4i6ki4 + 1, 2):
+            b2 = (a * a) - nk4
+            if is_square(b2):
+                b = isqrt(b2)
+                p = gcd(a + b, n)
+                q = gcd(a - b, n)
+                return p, q
+    return []
+
+def tests():
+  N = [29 * 89, 3*11, 3141592651*3141592661, 7919*10861]
+  N += [94738740796943840961823530695778701408987757287583492665919730017973847138345511139064596113422435977583856843887008168202003855906708039013487349390571801141407245039011598810542232029634564848797998534872251549660454277336502838185642937637576121533945369150901808833844341421315429263207612372324026271327]
+  N += [9733382803370256893136109840971590971460094779242334919432347801491641617443615856221168611138933576118196795282443503609663168324106758595642231987246769*9733382803370256893136109840971590971460094779242334919432347801491641617443615856221168611138933576118196795282443503609663168324106758595642231987248907]
+  N=sorted(N)
+  for n in N:
+    p,q = lehman(n)
+    print(n,p,q)
+    assert p*q == n
+
+if __name__ == "__main__":
+  tests()
+

From 05cd52eb89d52ed3341920ceea786c0275994054 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Thu, 5 Jan 2023 16:49:22 -0300
Subject: [PATCH 49/56] Add kraitchik algorithm

---
 kraitchik.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)
 create mode 100644 kraitchik.py

diff --git a/kraitchik.py b/kraitchik.py
new file mode 100644
index 0000000..ff295d7
--- /dev/null
+++ b/kraitchik.py
@@ -0,0 +1,21 @@
+from gmpy2 import *
+
+def kraitchik(n):
+    x = isqrt(n)
+    while True:
+        k = 1
+        s = x * x - k * n
+        while s >= 0:
+            if is_square(s):
+                y = isqrt(s)
+                z = x + y
+                w = x - y
+                if z % n != 0 and w % n != 0:
+                    return gcd(z,n),gcd(w,n)
+            s = x * x - k * n
+            k += 1
+        x+=1
+
+print(kraitchik(2041))
+print(kraitchik(65537*65539))
+

From f8a109cf27aaba4fc6ebe0fda5bd81fd79efdb54 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Thu, 5 Jan 2023 18:21:45 -0300
Subject: [PATCH 50/56] add carmichael integer factorization algorithm

---
 carmichael.py | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)
 create mode 100644 carmichael.py

diff --git a/carmichael.py b/carmichael.py
new file mode 100644
index 0000000..32b6170
--- /dev/null
+++ b/carmichael.py
@@ -0,0 +1,26 @@
+from gmpy2 import *
+
+def factor(N):
+  """
+  Algorithm described in the wagstaf-joy of factoring book.
+  """
+  f = N1 = N - 1
+  while f & 1 == 0: f >>= 1
+  a = 2
+  while a <= N1:
+    r1 = powmod(a, f << 1, N) 
+    if r1 == 1:
+      r = powmod(a, f, N) 
+      p = gcd(r - 1, N)
+      q = gcd(r + 1, N)
+      if (q > p > 1): #and (p * q == N):
+        return factor(p) + factor(q)
+    a = next_prime(a)
+  return [N]
+
+def tests():
+  for n in [561,6601, 8911,23224518901,1901 * 3701,6601 * 8911,52193 * 3701]:
+    print(n, factor(n))
+
+if __name__ == "__main__":
+  tests()

From 2c439f59c55ace92d453ef9f8fc29f7b1cc02ef9 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Fri, 10 Mar 2023 18:20:21 -0300
Subject: [PATCH 51/56] Update trial division algorithm and apply Euler's
 criterion for prime selection in QS algorithm.

---
 qs.sage | 151 ++++++++++++++++++--------------------------------------
 1 file changed, 48 insertions(+), 103 deletions(-)

diff --git a/qs.sage b/qs.sage
index af46b8d..09d6142 100644
--- a/qs.sage
+++ b/qs.sage
@@ -2,114 +2,59 @@
 # Code taken from https://facthacks.cr.yp.to/facthacks.sage
 # FactHacks is joint work by Daniel J. Bernstein, Nadia Heninger, and Tanja Lange.
 
-from gmpy2 import gcd, isqrt, is_prime, next_prime
+from gmpy2 import gcd, isqrt, is_prime, next_prime, is_power
 import itertools
-
-
-def prebuilt_params(bits):
-    """
-    Bounds estimation
-    borrowed from msieve/mpqs.c
-    """
-    if bits <= 64:
-        return [100, 40, 1 * 65536]
-    if bits <= 128: 
-        return [450, 40, 1 * 65536]
-    if bits <= 183:
-        return [2000, 40, 1 * 65536]
-    if bits <= 200: 
-        return [3000, 50, 1 * 65536]
-    if bits <= 212: 
-        return [5400, 50, 3 * 65536]
-    if bits <= 233:
-        return [10000, 100, 3 * 65536]
-    if bits <= 249:
-        return [27000, 100, 3 * 65536]
-    if bits <= 266:
-        return [50000, 100, 3 * 65536]
-    if bits <= 283:
-        return [55000, 80, 3 * 65536]
-    if bits <= 298:
-        return [60000, 80, 9 * 65536]
-    if bits <= 315:
-        return [80000, 150, 9 * 65536]
-    if bits <= 332:
-        return [100000, 150, 9 * 65536]
-    if bits <= 348:
-        return [140000, 150, 9 * 65536]
-    if bits <= 363:
-        return [210000, 150, 13 * 65536]
-    if bits <= 379:
-        return [300000, 150, 17 * 65536]
-    if bits <= 395:
-        return [400000, 150, 21 * 65536]
-    if bits <= 415:
-        return [500000, 150, 25 * 65536] # beyond this point you're crazy 
-    if bits <= 440:
-        return [700000, 150, 33 * 65536]
-    if bits <= 465:
-        return [900000, 150, 50 * 65536]
-    if bits <= 490:
-        return [1100000, 150, 75 * 65536]
-    if bits <= 512:
-        return [1300000, 150, 100 * 65536]
-    return [1300000, 150, 100 * 65536]
-
-
-def trial_division_minus_even_powers(n,P):
-    a = []
-    r = n
-    pw = 0
-    if r & 1 == 0:
-        while r & 1 == 0:
+import sys
+
+
+def trial_division_minus_even_powers(n, P):
+    a, r, pw = [], n, 0
+    while (r & 1) == 0:
+        pw += 1
+        r >>= 1
+    if (pw & 1 == 1): a.append(2)
+    for i in range(1, len(P)):
+        pw, p = 0, P[i]
+        while (r % p) == 0:
             pw += 1
-            r >>= 1
-        if pw & 1 != 0:
-            if 2 not in a:
-                a.append(2)
-    l=len(P)
-    i=0
-    while (i < l) and (P[i] <= n):
-        if r % P[i] == 0:
-            pw = 0
-            while r % P[i] == 0:
-                pw += 1
-                r //= P[i]
-            if pw & 1 != 0: 
-                if P[i] not in a:
-                    a.append(P[i])
-        else:
-            i += 1
-    return [a,r,n]
+            r //= p
+        if (pw & 1 == 1): a.append(p)
+    return a, r, n
 
 
 def minifactor(x, P):
-    if not is_prime(x) and not is_square(x):
-        p = trial_division_minus_even_powers(x, P)
-        if p[1] == 1: # if 1 x is B-smooth
-            return p
-
-
-def QS(N,differences,y):
-    i2N = isqrt(N)
-    X = [int(a * a - N) for a in range(i2N + 1, i2N + differences)]
-    P = list(primes(2, y))
-    F = list(filter(None, map(minifactor, X, itertools.repeat(P, len(X)))))
-
-    M = matrix(GF(2), len(F), len(P), lambda i, j:P[j] in F[i][0])
-    for K in M.left_kernel().basis():
-        I = [f for f, k in zip(F, K) if k == 1]
-        x = product([int(isqrt(f[2] + N)) for f in I])
-        y = isqrt(product([f[2] for f in I]))
-        g0, g1 = gcd(N, x - y), gcd(N, x + y)
-        if N > g0 > 1 or N > g1 > 1:
-            return g0, g1
-
+    p = trial_division_minus_even_powers(x, P)
+    if p[1] == 1: return p # if 1 then x is b-smooth
+
+
+def eulercriterion(N, p):
+  return pow(N, (p - 1) >> 1, p)
+
+
+def QS(N, B1):
+    sys.stderr.write("max prime: %d\n" % B1)
+    P = [p for p in list(primes(2, B1)) if eulercriterion(N, p) == 1]
+    i2n, B2, X, offset, lp = isqrt(N), 65536, [], 0, len(P)
+    while True:
+        S = i2n + offset
+        sys.stderr.write("Relations matix: [cols: %d x rows: %d]\n" % (lp, offset + B2))
+        X += [int(a * a - N) for a in range(S + 1, S + B2)]
+        X = list(filter(lambda x:not is_power(x) and not is_prime(x), X))
+        F = list(filter(None, map(minifactor, X, itertools.repeat(P, len(X)))))
+        M = matrix(GF(2), len(F), lp, lambda i, j:P[j] in F[i][0])
+        sys.stderr.write("performing linear algebra:\n")
+        for K in M.left_kernel().basis():
+            I = [f for f, k in zip(F, K) if k == 1]
+            if len(I) > 0:
+                sys.stderr.write("On Vector basis, taking square roots x and y.\n")
+                x = product([isqrt(f[2] + N) for f in I])
+                y = isqrt(product([f[2] for f in I]))
+                g0, g1 = int(gcd(N, x - y)), int(gcd(N, x + y))
+                if N > g0 > 1: return g0, N // g0
+                if N > g1 > 1: return g1, N // g1
+        offset += B2
 
 if __name__ == "__main__":
     N = int(sys.argv[1])
-    y ,_ ,diff = prebuilt_params(N.bit_length())
-    r = QS(N, diff * 10, y * 10)
-    if r != None:
-        print("%d, %d" % r)
-
+    B1 = int(exp(sqrt(log(N) * log(log(N)))) ** (1 / sqrt(2)))  
+    print(QS(N, B1))

From c1a29543d92c98c87e9534399ff9dd1137fdae77 Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Fri, 10 Mar 2023 19:04:55 -0300
Subject: [PATCH 52/56] Improve Fermat attack algorithm by replacing isqrt with
 isqrt_rem and is_square.

---
 fermat.py | 26 +++++++++++---------------
 1 file changed, 11 insertions(+), 15 deletions(-)

diff --git a/fermat.py b/fermat.py
index 8cf8605..12026d7 100644
--- a/fermat.py
+++ b/fermat.py
@@ -1,20 +1,16 @@
 import sys
-from gmpy2 import isqrt
+from gmpy2 import isqrt, isqrt_rem, is_square
 
 def fermat(n):
-        """Fermat attack"""
-        a = isqrt(n)
-        a = b
-        b2 = (a ** 2) - n
-        count = 0
-        while (b ** 2) != b2:
-            a = a + 1
-            b2 = (a ** 2) - n
-            b = isqrt(b2)
-            count += 1
-        p = a + b
-        q = a - b
-        assert n == p * q
-        return p, q
+    a, rem = isqrt_rem(n)
+    b2 = -rem
+    c0 = (a << 1) + 1
+    c = c0
+    while not is_square(b2):
+        b2 += c
+        c += 2
+    a = (c - 1) >> 1
+    b = isqrt(b2)
+    return a - b, a + b
 
 print(fermat(int(sys.argv[1])))

From 1bddcd8a2a47bc3f0f5b7cc14696609fc932da93 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Mon, 13 Mar 2023 00:51:39 -0300
Subject: [PATCH 53/56] Update qs.sage

---
 qs.sage | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/qs.sage b/qs.sage
index 09d6142..666362f 100644
--- a/qs.sage
+++ b/qs.sage
@@ -33,14 +33,14 @@ def eulercriterion(N, p):
 
 def QS(N, B1):
     sys.stderr.write("max prime: %d\n" % B1)
-    P = [p for p in list(primes(2, B1)) if eulercriterion(N, p) == 1]
+    P = [p for p in [*(primes(2, B1))] if eulercriterion(N, p) == 1]
     i2n, B2, X, offset, lp = isqrt(N), 65536, [], 0, len(P)
     while True:
         S = i2n + offset
         sys.stderr.write("Relations matix: [cols: %d x rows: %d]\n" % (lp, offset + B2))
         X += [int(a * a - N) for a in range(S + 1, S + B2)]
-        X = list(filter(lambda x:not is_power(x) and not is_prime(x), X))
-        F = list(filter(None, map(minifactor, X, itertools.repeat(P, len(X)))))
+        X = [*(filter(lambda x:not is_power(x) and not is_prime(x), X))]
+        F = [*(filter(None, map(minifactor, X, itertools.repeat(P, len(X)))))]
         M = matrix(GF(2), len(F), lp, lambda i, j:P[j] in F[i][0])
         sys.stderr.write("performing linear algebra:\n")
         for K in M.left_kernel().basis():

From c9482669e0e43c3205a72d668c6665747326270c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dar=C3=ADo=20Clavijo?= <clavijodario@gmail.com>
Date: Mon, 13 Mar 2023 00:52:55 -0300
Subject: [PATCH 54/56] Update qs.sage

---
 qs.sage | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/qs.sage b/qs.sage
index 666362f..84f4089 100644
--- a/qs.sage
+++ b/qs.sage
@@ -7,7 +7,7 @@ import itertools
 import sys
 
 
-def trial_division_minus_even_powers(n, P):
+def trial_division_odd_pows(n, P):
     a, r, pw = [], n, 0
     while (r & 1) == 0:
         pw += 1
@@ -23,7 +23,7 @@ def trial_division_minus_even_powers(n, P):
 
 
 def minifactor(x, P):
-    p = trial_division_minus_even_powers(x, P)
+    p = trial_division_odd_pows(x, P)
     if p[1] == 1: return p # if 1 then x is b-smooth
 
 

From bb2c3bcfcc97182e6659b53b0da2c904a1ca9797 Mon Sep 17 00:00:00 2001
From: Sourcery AI <>
Date: Mon, 27 Nov 2023 20:01:40 +0000
Subject: [PATCH 55/56] 'Refactored by Sourcery'

---
 MPQS.sage          |  75 ++++++++++++++-------------
 PisanoPeriod.py    | 124 ++++++++++++++++++++-------------------------
 qs.sage            |   3 +-
 shors_classic.sage |  25 +++++----
 4 files changed, 105 insertions(+), 122 deletions(-)

diff --git a/MPQS.sage b/MPQS.sage
index db70e6b..f0507d5 100644
--- a/MPQS.sage
+++ b/MPQS.sage
@@ -69,8 +69,6 @@ def prebuilt_params(bits):
         return [900000, 150, 50 * 65536]
     if bits <= 490:
         return [1100000, 150, 75 * 65536]
-    if bits <= 512:
-        return [1300000, 150, 100 * 65536]
     return [1300000, 150, 100 * 65536]
 
 
@@ -113,10 +111,10 @@ def choose_multiplier(n, prime_list):
         if knmod8 == 5:
             scores[i] -= M_LN2;
             break;
-        if knmod8 == 3 or knmod8 == 7:
+        if knmod8 in [3, 7]:
             scores[i] -= 0.5 * M_LN2;
             break;
-        # even multipliers start with a handicap
+            # even multipliers start with a handicap
     num_multipliers = i;
 
     # for the rest of the small factor base primes 
@@ -131,7 +129,7 @@ def choose_multiplier(n, prime_list):
             knmodp = (curr_mult * modp) % prime
             # if prime i is actually in the factor base
             #   for k * n ... */
-            if (knmodp == 0 or legendre(knmodp, prime) == 1):
+            if knmodp == 0:
 
                 """ ...add its contribution. A prime p con-
                    tributes log(p) to 1 in p sieve values, plus
@@ -147,18 +145,31 @@ def choose_multiplier(n, prime_list):
                    the prime divides k*n. In that case there 
                    is only one root """
 
-                if (knmodp == 0):
-                    scores[j] -= contrib
-                else:
-                    scores[j] -= 2 * contrib
+                scores[j] -= contrib
+            elif legendre(knmodp, prime) == 1:
+
+                """ ...add its contribution. A prime p con-
+                   tributes log(p) to 1 in p sieve values, plus
+                   log(p) to 1 in p^2 sieve values, etc. The
+                   average contribution of all multiples of p 
+                   to a random sieve value is thus
+                   log(p) * (1/p + 1/p^2 + 1/p^3 + ...)
+                   = (log(p) / p) * 1 / (1 - (1/p)) 
+                   = log(p) / (p-1)
+                   This contribution occurs once for each
+                   square root used for sieving. There are two
+                   roots for each factor base prime, unless
+                   the prime divides k*n. In that case there 
+                   is only one root """
 
+                scores[j] -= contrib if (knmodp == 0) else 2 * contrib
     # use the multiplier that generates the best score 
 
     best_score = 1000.0
     best_mult = 1
 
     #print(scores)
-   
+
     for i in range(0, num_multipliers):
         score = scores[i];
         if (score < best_score):
@@ -212,15 +223,14 @@ def trial_division(n, P):
     l = len(P)
     i = 0
     for p in P:
-        if r >= p:
-            if r % p == 0:
-                pw = 0
-                while r % p == 0:
-                    pw += 1
-                    r //= p
-                a.append((int(p),int(pw)))
-        else:
+        if r < p:
             break
+        if r % p == 0:
+            pw = 0
+            while r % p == 0:
+                pw += 1
+                r //= p
+            a.append((int(p),int(pw)))
     return a,r,n
 
 
@@ -234,10 +244,7 @@ def merge_powers(ppws):
             d[p] = pw
         else:
             d[p] += pw
-    tmp2 = []
-    for p in d:
-       tmp2.append((p,d[p]))
-    return tmp2
+    return list(d.items())
 
 
 def filter_out_even_powers(ppws):
@@ -250,11 +257,7 @@ def filter_out_even_powers(ppws):
             d[p] = pw
         else:
             d[p] += pw
-    tmp2 = []
-    for p in d:
-        if d[p] & 1 != 0: # only keep odd powers
-            tmp2.append(p)
-    return tmp2
+    return [p for p, value in d.items() if value & 1 != 0]
 
 
 def trial_division_minus_even_powers(n, P):
@@ -297,13 +300,12 @@ def is_power_logprime(n, log_primes):
     Given the precomputed logs of a set of primes 
     we can test if a number n is a perfect power of that prime.
     """
-    ispow = False
     for p in log_primes:
         a = log(n) / p
         b = int(a)
         if a == b:
             iwpow = True
-    return ispow
+    return False
 
 
 def minifactor4(x, P, smooth_base):
@@ -337,7 +339,7 @@ class Poly:
         #self.counter = 0
         self.early_factors = []
 
-        if A == None and B == None and C == None:
+        if A is None and B is None and C is None:
             self.create()
         else:
             print("given poly parameters: %d %d %d" % (A,B,C))
@@ -401,15 +403,15 @@ class Poly:
         Code borrowed https://github.com/cramppet/quadratic-sieve/blob/master/QS.py
         """        
         A = self.A
-        B = self.B
         C = self.C
         n = self.n
         p = self.root_A
         if A > 1:
             g = 1
+            ainv = 1
+            B = self.B
             while g == 1:
                 p = next_prime(p)
-                ainv = 1
                 if A != 1:
                     g, inv, _ = gcdext(A, p)
                     if g != 1:
@@ -442,9 +444,8 @@ class Poly:
         """
         Return the string representation of the constructed polynomial
         """
-        m = "F(X) = %d X ^ 2 + %d X + %d with minima: %d" % (self.A, self.B, self.C, self.minima) 
-        m = m.replace("+ -","- ")
-        return m
+        m = "F(X) = %d X ^ 2 + %d X + %d with minima: %d" % (self.A, self.B, self.C, self.minima)
+        return m.replace("+ -","- ")
 
 
     def __add__(self, other):
@@ -486,9 +487,7 @@ class Poly:
         """
         Hashes unique values of the polynomial to construct an python internal id.
         """
-        h = hash("%d-%d-%d" % (self.A,self.B,self.C))
-        #print("hash: %s" % h)
-        return h
+        return hash("%d-%d-%d" % (self.A,self.B,self.C))
 
 
     def __eq__(self, other):
diff --git a/PisanoPeriod.py b/PisanoPeriod.py
index d6631a7..9316848 100644
--- a/PisanoPeriod.py
+++ b/PisanoPeriod.py
@@ -23,14 +23,10 @@ def _fib_res(self,n,p):
         """ fibonacci sequence nth item modulo p """
         if n == 0:
             return (0, 1)
-        else:
-            a, b = self._fib_res(n >> 1,p)
-            c = mod((mod(a, p) * mod(((b << 1) - a), p)), p)
-            d = mod((powmod(a, 2, p) + powmod(b, 2, p)), p)
-            if n & 1 == 0: 
-                return (c, d)
-            else:
-                return (d, mod((c + d), p))
+        a, b = self._fib_res(n >> 1,p)
+        c = mod((mod(a, p) * mod(((b << 1) - a), p)), p)
+        d = mod((powmod(a, 2, p) + powmod(b, 2, p)), p)
+        return (c, d) if n & 1 == 0 else (d, mod((c + d), p))
 
 
     def get_n_mod_d(self,n,d, use = 'mersenne'):
@@ -45,9 +41,7 @@ def get_n_mod_d(self,n,d, use = 'mersenne'):
 
     def ranged_period(self, N, start, stop, look_up_dest):
         print("ranged_period (%d,%d) start" % (start,stop))
-        tmp_look_up = {}
-        for x in range(start, stop):
-            tmp_look_up[self.get_n_mod_d(x, N)] = x
+        tmp_look_up = {self.get_n_mod_d(x, N): x for x in range(start, stop)}
         look_up_dest.update(tmp_look_up)
         #look_up_dest |= tmp_look_up
         print("ranged_period (%d,%d) end" % (start,stop))
@@ -56,24 +50,22 @@ def ranged_period(self, N, start, stop, look_up_dest):
 
     def get_period_bigint(self, N, min_accept, xdiff, verbose = False):            
         search_len = int(pow(N, (1.0 / 6) / 100))
-        
-        if search_len < min_accept:
-            search_len = min_accept
-  
+
+        search_len = max(search_len, min_accept)
         if self.verbose:
             print('Search_len: %d, log2(N): %d' % (search_len,int(log2(N))))
-        
+
         starttime = time.time()
-        diff = xdiff 
+        diff = xdiff
         p_len = int((len(str(N)) + diff) >> 1) + 1
-        begin = N - int('9'*p_len) 
+        begin = N - int('9'*p_len)
         if begin <= 0:
             begin = 1
         end = N + int('9' * p_len)
-    
+
         if self.verbose:    
             print('Search begin: %d, end: %d'%(begin, end))
-        
+
         if self.multitasked and search_len > 1000:
             C = cpu_count() * 2
             search_work = search_len // C
@@ -84,10 +76,10 @@ def get_period_bigint(self, N, min_accept, xdiff, verbose = False):
             if self.verbose:
                 print("Precompute LUT with %d tasks..." % C)
 
-            inputs = []
-            for x in range(0, search_len, search_work):
-                inputs += [(N, x, x + search_work, look_up)]
-
+            inputs = [
+                (N, x, x + search_work, look_up)
+                for x in range(0, search_len, search_work)
+            ]
             workpool = Pool(C)
 
             with workpool:
@@ -103,10 +95,10 @@ def get_period_bigint(self, N, min_accept, xdiff, verbose = False):
         if self.verbose:
             print("LUT creation ended size: %d..." % len(look_up))
             print("Searching...")
-        
 
-        while True:       
-            randi = random.randint(begin,end)            
+
+        while True:   
+            randi = random.randint(begin,end)
             res = self.get_n_mod_d(randi, N)
             if res > 0:
                 #print(res, res in look_up)
@@ -122,13 +114,15 @@ def get_period_bigint(self, N, min_accept, xdiff, verbose = False):
                         else:
                             if self.verbose:
                                 print('For N = %d\n Found res: %d, res_n: %d , T: %d\n but failed!' % (N, res, res_n, T))
-            else:
-                if randi & 1 == 0:
-                    T = randi
-                    td = int(time.time() - starttime)
-                    if self.verbose:
-                        print('First shot, For N = %d Found T:%d, randi: %d, time used %f secs.' % (N , T, randi, td))
-                    return td, T, randi
+            elif randi & 1 == 0:
+                td = int(time.time() - starttime)
+                T = randi
+                if self.verbose:
+                    print(
+                        'First shot, For N = %d Found T:%d, randi: %d, time used %f secs.'
+                        % (N, T, T, td)
+                    )
+                return td, T, T
 
   
     def _trivial_factorization_with_n_phi(self, N, T):
@@ -142,24 +136,16 @@ def _trivial_factorization_with_n_phi(self, N, T):
 
         if phi2m4d > 0:
             iphi2m4d = isqrt(phi2m4d)
-            p1.append(phi + iphi2m4d)
-            p1.append(phi - iphi2m4d)
-
+            p1.extend((phi + iphi2m4d, phi - iphi2m4d))
         if phi2p4d > 0:
             iphi2p4d = isqrt(phi2p4d)
-            p1.append(phi + iphi2p4d)
-            p1.append(phi - iphi2p4d)
-
+            p1.extend((phi + iphi2p4d, phi - iphi2p4d))
         if phi2m4d > 0:
             iphi2m4d = isqrt(phi2m4d)
-            p1.append(-phi + iphi2m4d)
-            p1.append(-phi - iphi2m4d)
-
+            p1.extend((-phi + iphi2m4d, -phi - iphi2m4d))
         if phi2p4d > 0:
             iphi2p4d = isqrt(phi2p4d)
-            p1.append(-phi + iphi2p4d)
-            p1.append(-phi - iphi2p4d)
-
+            p1.extend((-phi + iphi2p4d, -phi - iphi2p4d))
         for p in p1:
             g = gcd((p >> 1),N)
             if N > g > 1:
@@ -224,7 +210,7 @@ def test(Ns,B2=0):
         if P != None:
             phi = (P[0]-1) * (P[1]-1)
             print("phi(N): %d" % phi)
-            print("factors: %s" % str(P))
+            print(f"factors: {str(P)}")
             ff += 1
         td = time.time() - ti
         ttd = time.time() - tti
@@ -235,27 +221,27 @@ def test(Ns,B2=0):
 
 
 def test2():
-  Fib = Fibonacci()
-  N = 384237701921
-  B1s = [10**x for x in range(6,3,-1)]
-  B2s = [0,2,4,6,8]
-  n=1
-  tti = time.time()
-  for B1 in B1s:
-      for B2 in B2s:
-          ti = time.time()
-          print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,int(log2(N)),B1,B2))
-          P = Fib.factorization(N,B1,B2)
-          if P != None:
-              phi = (P[0]-1) * (P[1]-1)
-              print("phi(N): %d" % phi)
-              print("factors: %s" % str(P))
-              ff += 1
-          td = time.time() - ti
-          ttd = time.time() - tti
-          print("Runtime: %f\nFully factored:%d of %d" % (td,ff,l))
-          print("Runtime total: %f" % (ttd))
-          n += 1
+    Fib = Fibonacci()
+    N = 384237701921
+    B1s = [10**x for x in range(6,3,-1)]
+    B2s = [0,2,4,6,8]
+    n=1
+    tti = time.time()
+    for B1 in B1s:
+        for B2 in B2s:
+            ti = time.time()
+            print("Test: %d, N: %d, log2(N): %d, B1: %d, B2: %d" % (n, N,int(log2(N)),B1,B2))
+            P = Fib.factorization(N,B1,B2)
+            if P != None:
+                phi = (P[0]-1) * (P[1]-1)
+                print("phi(N): %d" % phi)
+                print(f"factors: {str(P)}")
+                ff += 1
+            td = time.time() - ti
+            ttd = time.time() - tti
+            print("Runtime: %f\nFully factored:%d of %d" % (td,ff,l))
+            print("Runtime total: %f" % (ttd))
+            n += 1
 
 
 def test3(N, B2 = 0):
@@ -291,7 +277,7 @@ def test4(l,B2=0):
         if P != None:
             phi = (P[0]-1) * (P[1]-1)
             print("phi(N): %d" % phi)
-            print("factors: %s" % str(P))
+            print(f"factors: {str(P)}")
             ff += 1
         td = time.time() - ti
         ttd = time.time() - tti
diff --git a/qs.sage b/qs.sage
index 84f4089..aa393e1 100644
--- a/qs.sage
+++ b/qs.sage
@@ -44,8 +44,7 @@ def QS(N, B1):
         M = matrix(GF(2), len(F), lp, lambda i, j:P[j] in F[i][0])
         sys.stderr.write("performing linear algebra:\n")
         for K in M.left_kernel().basis():
-            I = [f for f, k in zip(F, K) if k == 1]
-            if len(I) > 0:
+            if I := [f for f, k in zip(F, K) if k == 1]:
                 sys.stderr.write("On Vector basis, taking square roots x and y.\n")
                 x = product([isqrt(f[2] + N) for f in I])
                 y = isqrt(product([f[2] for f in I]))
diff --git a/shors_classic.sage b/shors_classic.sage
index a9fcab0..89b3164 100644
--- a/shors_classic.sage
+++ b/shors_classic.sage
@@ -22,24 +22,23 @@ def multiplicativeOrder(a,n):
    return -1
 
 def shor_factor(N,explain=False):
-  """ Shor's algorithm recreated in a classical computer
+   """ Shor's algorithm recreated in a classical computer
   Premise: given a ^ P - 1 % N  == 0, where a is random, and p = multOrd(a,N) we get a factor o N.
   Sage's implementation of multiplicative order first factors the integer (N) in question with PARI.
   In other words it can be said that the complexity of multiplicativeOrder problem 
   is equal to the complexity of the factorization problem.
   This makes this algorithm only a showcase for studying purposes and not an improvement of performance in any way."""
 
-  n = 1
-  p = q = None
-  while True:
-    a = random.randint(2,N-1)
-    b = gcd(a,N)
-    if explain:
-      print("Iter: %d" % n)
-      print("gcd(%d, %d) = %d" % (a,N,b))
-    if 1 < b < N:
-      return N//b,b
-    else:
+   n = 1
+   p = q = None
+   while True:
+      a = random.randint(2,N-1)
+      b = gcd(a,N)
+      if explain:
+        print("Iter: %d" % n)
+        print("gcd(%d, %d) = %d" % (a,N,b))
+      if 1 < b < N:
+         return N//b,b
       r = Mod(a,N).multiplicative_order() # Equivalent to quantum part of the algorithm, complexity: O(N).
       if explain:
         print("multOrd(%d, %d) = %d" % (a,N,r))
@@ -61,7 +60,7 @@ def shor_factor(N,explain=False):
           if p != None and q != None:
             assert p*q == N
             return p,q
-    n+=1
+      n+=1
 
 if __name__ == "__main__":
   """

From e5db469fe08b31bc100e8c852528d28981dfb2dd Mon Sep 17 00:00:00 2001
From: dclavijo <clavijodario@gmail.com>
Date: Tue, 18 Jun 2024 18:10:43 -0300
Subject: [PATCH 56/56] add mersenne factorization

---
 mersenne.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)
 create mode 100644 mersenne.py

diff --git a/mersenne.py b/mersenne.py
new file mode 100644
index 0000000..6cbd869
--- /dev/null
+++ b/mersenne.py
@@ -0,0 +1,22 @@
+from sympy.core.power import isqrt
+from sympy import isprime
+from functools import cache
+
+@cache
+def factorMersenne(n):
+  " Lazzy programming implementation "
+  if isprime(n):  yield n
+  p = n.bit_length()
+  for q in range(2*p + 1, isqrt(n) + 1, 2*p):
+    if [0,1,0,0,0,0,0,1][q % 8] and pow(2, p, q) == 1:
+      if n % q == 0:
+        yield from factorMersenne(q)
+        yield from factorMersenne(n//q)
+
+def g():
+  for x in range(1, 939):
+    n = (1 << x)-1
+    print(x, n,[p for p in factorMersenne(n)])
+    print(factorMersenne.cache_info())
+
+g()
