Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Sign up
Appearance settings

Commit 422b4a0

Browse files
Add files via upload
1 parent 12afdc2 commit 422b4a0

24 files changed

+4058
-0
lines changed

‎README.txt‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
Thanks for looking at this project.
2+
3+
This is a project I started a few years ago, in 2019-2020.
4+
I re-worked it recently to merge some versions, make it more readable, add comments, and get all the sources I used.
5+
6+
7+
##### Introduction ######
8+
9+
This project is a Python implementation of the General number field sieve to factor integers.
10+
11+
Here are some detailed characteristics :
12+
- Kleinjung 2006 algorithm to find good polynomials
13+
- Double large prime variation used to collect relations for both algebraic and rational side
14+
- Batch smoothness test and naive smoothness test available
15+
- Union-find algorithm to detect cycles. Leveraging the graph structures (one central hub) to fasten the cycle finding when I know there is one.
16+
- Parallel sieving allowed
17+
- Parallel Polynomial search
18+
- gaussian elimination, block Lanczos, Wiedemann algorithms available for the linear algebra step (no parallelization available for now)
19+
- Lifting and Couveignes methods for computing the algebraic square root, no Montgomery algorithm available for now. (no parallelization for now)
20+
21+
##### Sources #####
22+
23+
Overall algorithm:
24+
- "Prime numbers, a computational perspective" by Richard Crandall and Carl Pomerance (really good): https://link.springer.com/book/10.1007/0-387-28979-8
25+
- "The Development of the Number Field Sieve" by many (really really good): https://link.springer.com/book/10.1007/BFb0091534
26+
- "A beginner's guide to the general number field sieve" by Michael Case: https://www.cs.umd.edu/~gasarch/TOPICS/factoring/NFSmadeeasy.pdf
27+
28+
Kleinjung polynomials search algorithm:
29+
- "On polynomial selection for the number field sieve" by THORSTEN KLEINJUNG: https://www.ams.org/journals/mcom/2006-75-256/S0025-5718年06月01日870-9/S0025-5718年06月01日870-9.pdf
30+
31+
Polynomial metric Ep_score:
32+
- "A new ranking function for polynomial selection in the number field sieve" by Nicolas David and Paul Zimmerman: https://inria.hal.science/hal-02151093v4/document
33+
34+
Double large prime
35+
- "Factoring with two large primes" by Arjen K. Lenstra and Mark S. Manasse: https://scispace.com/pdf/factoring-with-two-large-primes-1lk9719aco.pdf
36+
37+
Batch smoothness test:
38+
- "HOW TO FIND SMOOTH PARTS OF INTEGERS" by DANIEL J. BERNSTEIN: https://cr.yp.to/factorization/smoothparts-20040510.pdf
39+
40+
Gaussian elimination:
41+
- Took the gaussian elimination code from https://github.com/skollmann/PyFactorise/blob/master/factorise.py#L39
42+
43+
Block Lanczos:
44+
- "A Block Lanczos Algorithm for Finding Dependencies over GF(2)" by Peter L. Montgomery: https://scispace.com/pdf/a-block-lanczos-algorithm-for-finding-dependencies-over-gf-2-ezdu2qt0pp.pdf
45+
- "A modified block Lanczos algorithm with fewer vectors" by Emmanuel Thomé (really good): https://eprint.iacr.org/2016/329.pdf
46+
47+
Wiedemann algorithm:
48+
- "SOLVING HOMOGENEOUS LINEAR EQUATIONSOVER GF(2) VIA BLOCK WIEDEMANN ALGORITHM" by Don Coppersmith: https://www.ams.org/journals/mcom/1994-62-205/S0025-5718-1994-1192970-7/S0025-5718-1994-1192970-7.pdf
49+
50+
Square root algorithms
51+
- "Computing a square root for the number field sieve" by Jean-Marc Couveignes: https://www.math.u-bordeaux.fr/~jcouveig/publi/Cou94-2.pdf
52+
53+
##### running the algortihm #####
54+
55+
When running the main.py file, the required argument is --n the number you want to factor.
56+
57+
A suitable command looks like "python main.py --n [your_number]"
58+
59+
##### config parameters #####
60+
61+
I have defined a few parameters directly you can edit in the config file.
62+
63+
- batch_smooth_test: If True, the batch smoothness test from Bernstein will be used. Otherwise, the naive (trial division from the prime factor base) approach will be used. Batch smoothness test is generally faster
64+
- gaussian_pivot: If True, gaussian elimination will ALWAYS be used in the linear algebra step. If you want to use block Lanczos of Wiedemann algorithms, this has to be set to False.
65+
- lanczos: If True, the block Lanczos algorithm will be used. If False, the Wiedemann algorithm will be used. No matter its value, if gaussian_pivot is True then gaussian elimination will be performed.
66+
- square_root_couveignes: If True, the Couveignes algorithm is run for the square root step. It requires many primes to be inert. If False, the lifting algorithm is used
67+
- large_primes_constant: Define the constant that will be multiplied by the last prime in the factor base to obtain the bound for the single large primes. For the double large primes, the bound is this constant multiplied by the last prime in the factor base squared. Both the algebraic and rational sides have the same bounds.
68+
- block_size: If Block Lanczos or Wiedemann algorithm are used, this sets the block size. In the Wiedemann algorithm, this is only used to compute many matrix-vector products at once, and not to compute the matrix generator of the matrix sequence you obtain.
69+
- poly_search_nb_poly_coarse_eval: number of polynomials to be generated before doing precise ranking.
70+
- poly_search_nb_poly_precise_eval: number of polynomials to be kept for precise evaluation.
71+
- poly_search_prime_bound: maximum size of primes used in the polynomial generation
72+
- poly_search_nb_roots: number of roots to use in the Kleinjung polynomial generation algorithm. It is the l parameter in the original paper.
73+
- poly_search_multiplier: leading coefficient of generated polynomials is always a multiple of this multiplier.
74+
- NB_CPU_POLY_SELECTION: Sets the number of cpu used for running polynomial search. One cpu is always kept as a "leader" that collect polynomial candidates from the workers.
75+
- NB_CPU_SIEVE: Sets the number of cpu used for sieving. One cpu is always kept as a "leader" that collects relations from the sievers, and the rest of the cpus sieve, test for smoothness, and send their candidate relations to the leader.
76+
77+
##### General discussion #####
78+
79+
I tried overall to use as few exernal libraries as possible. While this makes the code way harder to write for me, and understand for you,
80+
I viewed it as a good challenge to write good code, and get to know hard and foreign concepts to me. It took a lot of work, but I got a bug free code,
81+
and got much more knowledgable on the subroutines thanks to this.
82+
83+
Here are some points I consider working on at some point:
84+
- I still have some work to do to implement the state of the art polynomial search algorithms.
85+
Right now this generates decent polynomials, but it seems that for high degree, eg 5 or 6)
86+
there are some useful techniques (rational approximation of root, translations and rotations).
87+
88+
- I have to implement the Montgomery algorithm for the algebraic square root computation.
89+
For the small numbers I have used my algorithm on, it is not a bottleneck yet.
90+
91+
- I have to optimize the Union-find algorithm for cycle detection, the version I have right now
92+
does not seem to be perfect, as I still do many operations to connect two distinct connected components.
93+
94+
- No matter how hard I tried, I am stuck on understanding the block Wiedemann algorithm. For now, the best I can do is the scalar one, with some optimizations.
95+
Namely, I use binary encoding of the blocks of vectors to compute matrix-vector product very efficiently. Then, for each scalar sequence, I compute
96+
its generator using scalar Berlekamp-Massey algorithm, and use it to update the minimal polynomial of the matrix. I consider working some way to
97+
compute the matrix generator of the matrix sequence through some algorithm (matrix Berlekamp-Massey or some other).
98+
99+
- I have to make the linear algebra step more parallel. Even though it is the fastest of the two main steps of the algorithm, no one likes to wait if it is
100+
possible to faster.
101+
102+
Here are the next steps:
103+
104+
- Some optimizations may be possible here and there, but I have already a big amount of work to get fast code. The next step on this topic is switching to C.
105+
Python is intrisically limited on the topic of speed, as it is compiled on the fly. Just like this project, I have already done part of the C project
106+
a while ago. I have to rework it to get something on the same level of quality as this Python version.

‎config/config.ini‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
[DEFAULT]
2+
batch_smooth_test=True
3+
gaussian_pivot=False
4+
lanczos=True
5+
square_root_couveignes=False
6+
large_primes_constant=100
7+
block_size=8
8+
poly_search_nb_poly_coarse_eval=1000
9+
poly_search_nb_poly_precise_eval=200
10+
poly_search_prime_bound=300
11+
poly_search_nb_roots=3
12+
poly_search_multiplier=1
13+
NB_CPU_POLY_SELECTION=4
14+
NB_CPU_SIEVE=1

‎requirements.txt‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
random
2+
math
3+
sys
4+
time
5+
datetime
6+
argparse
7+
configparser
8+
numpy

‎src/NFS.py‎

Lines changed: 218 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,218 @@
1+
# This is the file that execute the Quadratic Sieve algorithm
2+
3+
## import libraries
4+
5+
ZQDFEGSRHTDR = 35279032783773322923326661023
6+
7+
import math, time, os
8+
from datetime import datetime
9+
import log
10+
import parse_config
11+
from utils import *
12+
from polynomial_functions import *
13+
from generate_primes import *
14+
import mono_cpu_polynomial_selection
15+
import multi_cpu_polynomial_selection
16+
import mono_cpu_sieve
17+
import multi_cpu_sieve
18+
from relations import *
19+
from build_matrix import *
20+
from gaussian_elimination import *
21+
from wiedemann import *
22+
from block_lanczos import *
23+
import compute_solutions
24+
25+
## Set path to find config file
26+
27+
CONFIG_PATH = os.path.abspath("../config/config.ini")
28+
29+
## NFS functions
30+
31+
def initialize(n):
32+
33+
d = int(pow(3*math.log(n)/math.log(math.log(n)),1/3))
34+
first_B = int(math.exp(pow((8/9)*math.log(n),1/3)*pow(math.log(math.log(n)),2/3)))
35+
B = 10*first_B//70
36+
37+
return d, B
38+
39+
def initialize_2(f_x, n, m1, d, primes, leading_coeff):
40+
pairs_used,R_p,logs,tmp = [],[],[],10*d
41+
divide_leading,pow_div = [],[]
42+
for p in primes:
43+
if not leading_coeff%p:
44+
divide_leading.append(p)
45+
u = p
46+
while not leading_coeff%u: u *= p
47+
pow_div.append(u//p)
48+
logs.append(round(math.log2(p)))
49+
if p > tmp: R_p.append(find_roots_poly(f_x,p))
50+
else: R_p.append(fast_roots(f_x,p))
51+
if len(R_p[-1]) == d:
52+
if d&1: pairs_used.append([[leading_coeff*p],leading_coeff*pow(p,d),[[1,1]],p*m1,[1],[0,p],1])
53+
else: pairs_used.append([[leading_coeff*p],leading_coeff*p,[[1,1]],p*m1,[1],[0,p],1])
54+
for p in divide_leading:
55+
for i in range(len(pairs_used)): pairs_used[i].append(True)
56+
57+
B_prime = sum([len(i) for i in R_p]) + 1
58+
59+
return pairs_used, R_p, logs, divide_leading, pow_div, B_prime
60+
61+
def initialize_3(n, f_x, f_prime, const1, leading_coeff):
62+
k = 3*n.bit_length()
63+
Q = []
64+
q = const1+1
65+
if not q%2: q += 1
66+
while len(Q) < k:
67+
if is_prime(q) and leading_coeff%q:
68+
if not n%q: return q, n//q
69+
else:
70+
tmp = find_roots_poly(f_x,q)
71+
for r in tmp:
72+
if eval_mod(f_prime,r,q): Q.append([q,r])
73+
q += 2
74+
75+
return Q, k
76+
77+
def NFS(n):
78+
now = datetime.now()
79+
80+
LOG_PATH = os.path.abspath("../logs/log_"+str(now.year)+str(now.month)+str(now.day)+"_"+str(now.hour)+str(now.minute)+str(now.second)+".txt")
81+
82+
parameters = parse_config.parse_config(CONFIG_PATH)
83+
84+
flag_use_batch_smooth_test = parameters[0].lower() in ["true"]
85+
flag_gaussian_pivot = parameters[1].lower() in ["true"]
86+
flag_lanczos = parameters[2].lower() in ["true"]
87+
flag_square_root_couveignes = parameters[3].lower() in ["true"]
88+
const = int(parameters[4])
89+
BLOCK_SIZE = int(parameters[5])
90+
nb_poly_coarse_eval = int(parameters[6])
91+
nb_poly_precise_eval = int(parameters[7])
92+
prime_bound = int(parameters[8])
93+
nb_roots = int(parameters[9])
94+
multiplier = int(parameters[10])
95+
NB_CPU_POLY_SELECTION = int(parameters[11])
96+
NB_CPU_SIEVE = int(parameters[12])
97+
98+
n = int(n)
99+
100+
d, B = initialize(n)
101+
102+
primes = create_smooth_primes_base(B)
103+
for p in primes:
104+
if not n%p: return p,n//p
105+
106+
prod_primes = math.prod(primes)
107+
108+
const1, const2 = const*primes[-1], const*primes[-1]*primes[-1]
109+
110+
if NB_CPU_POLY_SELECTION == 1:
111+
f_x,m0,m1,tmp,_ = mono_cpu_polynomial_selection.poly_search(n,primes,nb_roots,prime_bound,multiplier,
112+
int(pow(n,1/(d+1))),d,nb_poly_coarse_eval,
113+
nb_poly_precise_eval,LOG_PATH)
114+
115+
else:
116+
f_x,m0,m1,tmp,_ = multi_cpu_polynomial_selection.poly_search(n,primes,nb_roots,prime_bound,multiplier,
117+
int(pow(n,1/(d+1))),d,nb_poly_coarse_eval,
118+
nb_poly_precise_eval,NB_CPU_POLY_SELECTION,LOG_PATH)
119+
120+
log.write_log(LOG_PATH, "poly search completed, parameters : m0 = "+str(m0)+" ; m1 = "+str(m1)+" ; d = "+str(d)+"\n")
121+
122+
f_x, m0, M = mono_cpu_polynomial_selection.evaluate_polynomial_quality(f_x,B,m0,m1,primes,LOG_PATH)
123+
124+
leading_coeff = f_x[0]
125+
zeros_f = get_complex_roots(f_x)
126+
zeros = [leading_coeff*i for i in zeros_f]
127+
f_prime = get_derivative(f_x)
128+
129+
g = [1,f_x[1]]
130+
for i in range(2,len(f_x)): g.append(f_x[i]*pow(leading_coeff,i-1))
131+
g_prime = get_derivative(g)
132+
133+
g_prime_sq,g_prime_eval = div_poly(poly_prod(g_prime,g_prime),g),pow(leading_coeff, d-2, n)*eval_F(m0,m1,f_prime,d-1)%n
134+
135+
delta = []
136+
for r in range(d):
137+
tmp = []
138+
for i in range(d):
139+
delt = 0
140+
for j in range(d-i):
141+
if zeros[r].real == 0 or zeros[r].imag == 0 or j == 0: delt += (abs(g[-i-j-2])*pow(leading_coeff,j)*pow(my_norm(zeros_f[r]),j))
142+
else: delt += (abs(g[-i-j-2])*pow(leading_coeff,j)*pow(my_norm(zeros_f[r]),j))+1
143+
tmp.append(delt)
144+
delta.append(tmp)
145+
146+
147+
inert_set = []
148+
for p in primes:
149+
if m1%p and leading_coeff%p and irreducibility(g,p): inert_set.append(p)
150+
151+
pairs_used, R_p, logs, divide_leading, pow_div, B_prime = initialize_2(f_x, n, m1, d, primes, leading_coeff)
152+
153+
Q, k = initialize_3(n, f_x, f_prime, B, leading_coeff)
154+
155+
log.write_log(LOG_PATH, "components of factor base : ("+str(len(primes))+","+str(B_prime)+","+str(k)+","+str(len(divide_leading))+")")
156+
log.write_log(LOG_PATH, "starting with "+str(len(pairs_used))+" free relations")
157+
158+
if NB_CPU_SIEVE <= 1:
159+
if NB_CPU_SIEVE < 1: print("NB_CPU parameter incorrectly set. Must be > 0. Sieving with 1 CPU.")
160+
pairs_used, V = mono_cpu_sieve.find_relations(f_x,leading_coeff,g,primes,R_p,Q,B_prime,divide_leading,
161+
prod_primes,pow_div,pairs_used,const1,const2,logs,m0,m1,
162+
M,d,n,flag_use_batch_smooth_test,LOG_PATH)
163+
164+
else:
165+
pairs_used, V = multi_cpu_sieve.find_relations(f_x,leading_coeff,g,primes,R_p,Q,B_prime,divide_leading,
166+
prod_primes,pow_div,pairs_used,const1,const2,logs,m0,m1,
167+
M,d,n,flag_use_batch_smooth_test,LOG_PATH,NB_CPU_SIEVE)
168+
169+
print("")
170+
log.write_log(LOG_PATH, "sieving complete, building matrix...")
171+
if flag_gaussian_pivot:
172+
matrix = build_dense_matrix(pairs_used, primes, R_p, Q, divide_leading)
173+
length = len(matrix[0])
174+
matrix, N, U = siqs_build_matrix_opt(matrix)
175+
else:
176+
matrix = build_sparse_matrix(pairs_used,primes,R_p,Q,divide_leading)
177+
matrix = transpose_sparse(matrix, V)
178+
log.write_log(LOG_PATH, "matrix built "+str(len(pairs_used))+"x"+str(V)+" reducing...")
179+
#if not flag_gaussian_pivot:
180+
# matrix,pairs_used = reduce_sparse_matrix(matrix,pairs_used)
181+
log.write_log(LOG_PATH, "matrix built "+str(len(pairs_used))+"x"+str(V)+" finding kernel...")
182+
183+
time_1 = datetime.now()
184+
185+
if flag_gaussian_pivot:
186+
null_space = siqs_solve_matrix_opt(matrix, N, U)
187+
log.write_log(LOG_PATH, str(len(null_space))+" kernel vectors found")
188+
for vec in null_space:
189+
vec = compute_solutions.convert_to_binary(vec, length)
190+
flag, res = compute_solutions.compute_factors(pairs_used,vec,n,primes,f_x,g,g_prime,g_prime_sq,g_prime_eval,
191+
m0,m1,leading_coeff,d,inert_set,zeros,delta,M,flag_square_root_couveignes,
192+
time_1,LOG_PATH)
193+
if flag: return res
194+
195+
else:
196+
if flag_lanczos:
197+
while True:
198+
null_space = block_lanczos(matrix,len(pairs_used),BLOCK_SIZE, LOG_PATH)
199+
200+
log.write_log(LOG_PATH, str(len(null_space))+" kernel vectors found")
201+
for vec in null_space:
202+
vec = compute_solutions.convert_to_binary_lanczos(vec, len(pairs_used))
203+
flag, res = compute_solutions.compute_factors(pairs_used,vec,n,primes,f_x,g,g_prime,g_prime_sq,g_prime_eval,
204+
m0,m1,leading_coeff,d,inert_set,zeros,delta,M,flag_square_root_couveignes,
205+
time_1,LOG_PATH)
206+
if flag: return res
207+
208+
else:
209+
mini_poly_estim = 1
210+
while True:
211+
null_space, mini_poly_estim = wiedemann(matrix,len(pairs_used),BLOCK_SIZE,mini_poly_estim)
212+
null_space = reduce_null_space_vectors(null_space)
213+
log.write_log(LOG_PATH, str(len(null_space))+" kernel vectors found")
214+
for vec in null_space:
215+
flag, res = compute_solutions.compute_factors(pairs_used,vec,n,primes,f_x,g,g_prime,g_prime_sq,g_prime_eval,
216+
m0,m1,leading_coeff,d,inert_set,zeros,delta,M,flag_square_root_couveignes,
217+
time_1,LOG_PATH)
218+
if flag: return res

‎src/__init__.py‎

Whitespace-only changes.

0 commit comments

Comments
(0)

AltStyle によって変換されたページ (->オリジナル) /