在基于二次筛法的因式分解方法中，找到的左零空间二进制 矩阵(计算模 2 的值)是关键步骤。 (这也是转置的零空间。)numpy 或 scipy 是否有工具可以快速执行此操作？
作为引用，这是我当前的代码:
# Row-reduce binary matrix
def binary_rr(m):
rows, cols = m.shape
l = 0
for k in range(min(rows, cols)):
print(k)
if l >= cols: break
# Swap with pivot if m[k,l] is 0
if m[k,l] == 0:
found_pivot = False
while not found_pivot:
if l >= cols: break
for i in range(k+1, rows):
if m[i,l]:
m[[i,k]] = m[[k,i]] # Swap rows
found_pivot = True
break
if not found_pivot: l += 1
if l >= cols: break # No more rows
# For rows below pivot, subtract row
for i in range(k+1, rows):
if m[i,l]: m[i] ^= m[k]
l += 1
return m
它几乎是高斯消元的简单实现，但由于它是用 python 编写的，因此速度非常慢。
最佳答案
qwr，我发现了一个非常快的高斯消除例程，它完成得如此之快以至于慢点是二次筛分或 SIQS 筛分步骤。高斯消元函数取自 https://raw.githubusercontent.com/skollmann/PyFactorise/master/factorise.py 处的 skollmans factorise.py。
我很快就会从头开始研究 SIQS/GNFS 实现，并希望用多线程和可能的 cython 为 python 编写一些 super 快速的东西。同时，如果你想要编译 C(Alpertons ECM 引擎)但使用 python 的东西，你可以使用:https://github.com/oppressionslayer/primalitytest/这要求您在使用 from sfactorint import p2ecm 导入 p2ecm 之前先 cd 到计算器目录并运行 make。有了它，您可以在几秒钟内分解 60 位数字。
# Requires sympy and numpy to be installed
# Adjust B and I accordingly. Set for 32 length number
# Usage:
# N=1009732533765251*1896182711927299
# factorise(N, 5000, 25000000) # Takes about 45-60 seconds on a newer computer
# N=1009732533765251*581120948477
# Linear Algebra Step finishes in 1 second, if that
# N=factorise(N, 5000, 2500000) # Takes about 5 seconds on a newer computer
# #Out[1]: 581120948477
import math
import numpy as np
from sympy import isprime
#
# siqs_ functions are the Gaussian Elimination routines right from
# skollmans factorise.py. It is the fastest Gaussian Elimination that i have
# found in python
#
def siqs_factor_from_square(n, square_indices, smooth_relations):
"""Given one of the solutions returned by siqs_solve_matrix_opt,
return the factor f determined by f = gcd(a - b, n), where
a, b are calculated from the solution such that a*a = b*b (mod n).
Return f, a factor of n (possibly a trivial one).
"""
sqrt1, sqrt2 = siqs_calc_sqrts(square_indices, smooth_relations)
assert (sqrt1 * sqrt1) % n == (sqrt2 * sqrt2) % n
return math.gcd(abs(sqrt1 - sqrt2), n)
def siqs_find_factors(n, perfect_squares, smooth_relations):
"""Perform the last step of the Self-Initialising Quadratic Field.
Given the solutions returned by siqs_solve_matrix_opt, attempt to
identify a number of (not necessarily prime) factors of n, and
return them.
"""
factors = []
rem = n
non_prime_factors = set()
prime_factors = set()
for square_indices in perfect_squares:
fact = siqs_factor_from_square(n, square_indices, smooth_relations)
if fact != 1 and fact != rem:
if isprime(fact):
if fact not in prime_factors:
print ("SIQS: Prime factor found: %d" % fact)
prime_factors.add(fact)
while rem % fact == 0:
factors.append(fact)
rem //= fact
if rem == 1:
break
if isprime(rem):
factors.append(rem)
rem = 1
break
else:
if fact not in non_prime_factors:
print ("SIQS: Non-prime factor found: %d" % fact)
non_prime_factors.add(fact)
if rem != 1 and non_prime_factors:
non_prime_factors.add(rem)
for fact in sorted(siqs_find_more_factors_gcd(non_prime_factors)):
while fact != 1 and rem % fact == 0:
print ("SIQS: Prime factor found: %d" % fact)
factors.append(fact)
rem //= fact
if rem == 1 or sfactorint_isprime(rem):
break
if rem != 1:
factors.append(rem)
return factors
def add_column_opt(M_opt, tgt, src):
"""For a matrix produced by siqs_build_matrix_opt, add the column
src to the column target (mod 2).
"""
M_opt[tgt] ^= M_opt[src]
def find_pivot_column_opt(M_opt, j):
"""For a matrix produced by siqs_build_matrix_opt, return the row of
the first non-zero entry in column j, or None if no such row exists.
"""
if M_opt[j] == 0:
return None
return lars_last_powers_of_two_trailing(M_opt[j] + 1)
def siqs_build_matrix_opt(M):
"""Convert the given matrix M of 0s and 1s into a list of numbers m
that correspond to the columns of the matrix.
The j-th number encodes the j-th column of matrix M in binary:
The i-th bit of m[i] is equal to M[i][j].
"""
m = len(M[0])
cols_binary = [""] * m
for mi in M:
for j, mij in enumerate(mi):
cols_binary[j] += "1" if mij else "0"
return [int(cols_bin[::-1], 2) for cols_bin in cols_binary], len(M), m
def siqs_solve_matrix_opt(M_opt, n, m):
"""
Perform the linear algebra step of the SIQS. Perform fast
Gaussian elimination to determine pairs of perfect squares mod n.
Use the optimisations described in [1].
[1] Koç, Çetin K., and Sarath N. Arachchige. 'A Fast Algorithm for
Gaussian Elimination over GF (2) and its Implementation on the
GAPP.' Journal of Parallel and Distributed Computing 13.1
(1991): 118-122.
"""
row_is_marked = [False] * n
pivots = [-1] * m
for j in range(m):
i = find_pivot_column_opt(M_opt, j)
if i is not None:
pivots[j] = i
row_is_marked[i] = True
for k in range(m):
if k != j and (M_opt[k] >> i) & 1: # test M[i][k] == 1
add_column_opt(M_opt, k, j)
perf_squares = []
for i in range(n):
if not row_is_marked[i]:
perfect_sq_indices = [i]
for j in range(m):
if (M_opt[j] >> i) & 1: # test M[i][j] == 1
perfect_sq_indices.append(pivots[j])
perf_squares.append(perfect_sq_indices)
return perf_squares
def sqrt_int(N):
Nsqrt = math.isqrt(N)
assert Nsqrt * Nsqrt == N
return Nsqrt
def siqs_calc_sqrts(square_indices, smooth_relations):
"""Given on of the solutions returned by siqs_solve_matrix_opt and
the corresponding smooth relations, calculate the pair [a, b], such
that a^2 = b^2 (mod n).
"""
res = [1, 1]
for idx in square_indices:
res[0] *= smooth_relations[idx][0]
res[1] *= smooth_relations[idx][1]
res[1] = sqrt_int(res[1])
return res
def quad_residue(a,n):
l=1
q=(n-1)//2
x = q**l
if x==0:
return 1
a =a%n
z=1
while x!= 0:
if x%2==0:
a=(a **2) % n
x//= 2
else:
x-=1
z=(z*a) % n
return z
def STonelli(n, p):
assert quad_residue(n, p) == 1, "not a square (mod p)"
q = p - 1
s = 0
while q % 2 == 0:
q //= 2
s += 1
if s == 1:
r = pow(n, (p + 1) // 4, p)
return r,p-r
for z in range(2, p):
#print(quad_residue(z, p))
if p - 1 == quad_residue(z, p):
break
c = pow(z, q, p)
r = pow(n, (q + 1) // 2, p)
t = pow(n, q, p)
m = s
t2 = 0
while (t - 1) % p != 0:
t2 = (t * t) % p
for i in range(1, m):
if (t2 - 1) % p == 0:
break
t2 = (t2 * t2) % p
b = pow(c, 1 << (m - i - 1), p)
r = (r * b) % p
c = (b * b) % p
t = (t * c) % p
m = i
return (r,p-r)
def build_smooth_relations(smooth_base, root_base):
smooth_relations = []
for xx in range(len(smooth_base)):
smooth_relations.append((root_base[xx], smooth_base[xx], xx))
return smooth_relations
def strailing(N):
return N>>lars_last_powers_of_two_trailing(N)
def lars_last_powers_of_two_trailing(N):
p,y=1,2
orign = N
#if orign < 17: N = N%16
N = N&15
if N == 1:
if ((orign -1) & (orign -2)) == 0: return orign.bit_length()-1
while orign&y == 0:
p+=1
y<<=1
return p
if N in [3, 7, 11, 15]: return 1
if N in [5, 13]: return 2
if N == 9: return 3
return 0
def build_matrix(factor_base, smooth_base):
factor_base = factor_base.copy()
factor_base.insert(0, 2)
sparse_matrix = []
col = 0
for xx in smooth_base:
sparse_matrix.append([])
for fx in factor_base:
count = 0
factor_found = False
while xx % fx == 0:
factor_found = True
xx=xx//fx
count+=1
if count % 2 == 0:
sparse_matrix[col].append(0)
continue
else:
if factor_found == True:
sparse_matrix[col].append(1)
else:
sparse_matrix[col].append(0)
col+=1
return np.transpose(sparse_matrix)
def get_mod_congruence(root, N, withstats=False):
r = root - N
if withstats==True:
print(f"{root} ≡ {r} mod {N}")
return r
def primes_sieve2(limit):
a = np.ones(limit, dtype=bool)
a[0] = a[1] = False
for (i, isprime) in enumerate(a):
if isprime:
yield i
for n in range(i*i, limit, i):
a[n] = False
def remove_singletons(XX):
no_singletons = []
for xx in XX:
if len(xx) != 1:
no_singletons.append(xx)
return no_singletons
def fb_sm(N, B, I):
factor_base, sieve_base, sieve_list, smooth_base, root_base = [], [], [], [], []
primes = list(primes_sieve2(B))
i,root=-1,math.isqrt(N)
for x in primes[1:]:
if quad_residue(N, x) == 1:
factor_base.append(x)
for x in range(I):
xx = get_mod_congruence((root+x)**2, N)
sieve_list.append(xx)
if xx % 2 == 0:
xx = strailing(xx+1) # using lars_last_modulus_powers_of_two(xx) bit trick
sieve_base.append(xx)
for p in factor_base:
residues = STonelli(N, p)
for r in residues:
for i in range((r-root) % p, len(sieve_list), p):
while sieve_base[i] % p == 0:
sieve_base[i] //= p
for o in range(len(sieve_list)):
# This is set to 350, which is only good for numbers
# of len < 32. Modify
# to be more dynamic for larger numbers.
if len(smooth_base) >= 350:
break
if sieve_base[o] == 1:
smooth_base.append(sieve_list[o])
root_base.append(root+o)
return factor_base, smooth_base, root_base
def isSquare(hm):
cr=math.isqrt(hm)
if cr*cr == hm:
return True
return False
def find_square(smooth_base):
for x in smooth_base:
if isSquare(x):
return (True, smooth_base.index(x))
else:
return (False, -1)
t_matrix=[]
primes=list(primes_sieve2(1000000))
def factorise(N, B=10000, I=10000000):
global primes, t_matrix
if isprime(N):
return N
for xx in primes:
if N%xx == 0:
return xx
factor_base, smooth_base, root_base = fb_sm(N,B,I)
issquare, t_matrix = find_square(smooth_base)
if issquare == True:
return math.gcd(math.isqrt(smooth_base[t_matrix])+get_mod_congruence(root_base[t_matrix], N), N)
t_matrix = build_matrix(factor_base, smooth_base)
smooth_relations = build_smooth_relations(smooth_base, root_base)
M_opt, M_n, M_m = siqs_build_matrix_opt(np.transpose(t_matrix))
perfect_squares = remove_singletons(siqs_solve_matrix_opt(M_opt, M_n, M_m))
factors = siqs_find_factors(N, perfect_squares, smooth_relations)
return factors
关于python - 在python中查找二进制矩阵的零空间，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/49287398/