· 9 years ago · Feb 02, 2017, 09:16 AM
1# factmsieve.py - A Python driver for GGNFS and MSIEVE
2#
3# Copyright (c) 2010, Brian Gladman
4#
5# All rights reserved.
6#
7# Redistribution and use in source and binary forms, with
8# or without modification, are permitted provided that the
9# following conditions are met:
10#
11# Redistributions of source code must retain the above
12# copyright notice, this list of conditions and the
13# following disclaimer.
14#
15# Redistributions in binary form must reproduce the
16# above copyright notice, this list of conditions and
17# the following disclaimer in the documentation and/or
18# other materials provided with the distribution.
19#
20# The names of its contributors may not be used to
21# endorse or promote products derived from this
22# software without specific prior written permission.
23#
24# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
25# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
26# WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
27# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
28# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
29# THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
30# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
31# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
33# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
34# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
35# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
37# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
38# POSSIBILITY OF SUCH DAMAGE.
39
40# This code is a conversion of the script factmsieve.pl
41# which is Copyright 2004, Chris Monico. His contribution
42# is acknowledged, as are those of all who contributed to
43# the original Perl script.
44#
45# I also acknowledge the invaluable help I have had from
46# Jeff Gilchrist in testing amd debugging this driver.
47# Without his support during its development, it would
48# never have reached a working state.
49#
50# The support of Wingware, who kindly donated their frst
51# class Python development environment, is acknowledged.
52
53from __future__ import print_function
54from __future__ import division
55
56import os, sys, random, re, functools, string, socket, signal
57import time, subprocess, gzip, glob, math, tempfile, datetime
58import atexit, threading, collections, multiprocessing, platform
59
60VERSION = '0.86'
61
62# Set binary directory paths
63GGNFS_PATH = 'C:\ggnfs\gnfs_test'
64MSIEVE_PATH = 'C:\ggnfs\gnfs_test'
65
66# Set the number of CPU cores and threads
67NUM_CORES = 4
68THREADS_PER_CORE = 1
69
70USE_CUDA = True
71GPU_NUM = 0
72MSIEVE_POLY_TIME_LIMIT = 0
73MIN_NFS_BITS = 264
74
75# msieve polynomial search leading coefficients limits
76# search from POLY_MIN_LC + 1 to POLY_MAX_LC inclusive
77POLY_MIN_LC = 0
78POLY_MAX_LC = 4000
79
80# Set global flags to control operation
81
82CHECK_BINARIES = True
83CHECK_POLY = True
84CLEANUP = False
85DOCLASSICAL = False
86NO_DEF_NM_PARAM = False
87PROMPTS = False
88SAVEPAIRS = True
89USE_KLEINJUNG_FRANKE_PS = False
90USE_MSIEVE_POLY = True
91VERBOSE = True
92
93# End of configuration options
94
95# number of msieve poly search threads to launch
96MS_THREADS = NUM_CORES * THREADS_PER_CORE
97# number of siever threads to launch
98SV_THREADS = NUM_CORES * THREADS_PER_CORE
99# number of linear algebra threads to launch
100LA_THREADS = NUM_CORES * THREADS_PER_CORE
101
102# ggnfs and msieve executable flie names
103
104MSIEVE = 'msieve'
105
106MAKEFB = 'makefb'
107PROCRELS = 'procrels'
108
109POL51M0 = 'pol51m0b'
110POL51OPT = 'pol51opt'
111POLYSELECT = 'polyselect'
112
113PLOT = 'autogplot.sh'
114
115# default parameter files
116
117DEFAULT_PAR_FILE = 'def-par.txt'
118DEFAULT_POLSEL_PAR_FILE = 'def-nm-params.txt'
119
120# temporary files
121
122PARAMFILE = '.params'
123RELSBIN = 'rels.bin'
124
125if sys.platform.startswith('win'):
126 EXE_SUFFIX = '.exe'
127else:
128 EXE_SUFFIX = ''
129 NICE_PATH = ''
130 MSIEVE = './' + MSIEVE
131
132# static global variables
133
134PNUM = 0
135LARGEP = 3
136LARGEPRIMES = '-' + str(LARGEP) + 'p'
137nonPrefDegAdjust = 12
138polySelTimeMultiplier = 1.0
139
140# poly5 parameters
141
142pol5_p = { 'max_pst_time':0, 'search_a5step':0, 'npr':0, 'normmax':0.0,
143 'normmax1':0.0, 'normmax2':0.0, 'murphymax':0.0 }
144
145# poly select parameters
146
147pols_p = { 'degree':0, 'maxs1':0, 'maxskew':0, 'goodscore':0.0,
148 'examinefrac':0.0, 'j0':0, 'j1':0, 'estepsize':0, 'maxtime':0 }
149
150# lattice sieve parameters
151
152lats_p = { 'rlim':0, 'alim':0, 'lpbr':0, 'lpba':0, 'mfbr':0, 'mfba':0,
153 'rlambda':0.0, 'alambda':0.0, 'qintsize':0, 'lss':1,
154 'siever': 'gnfs-lasieve4I10e', 'minrels':0, 'currels':0 }
155
156# classical sieve parameters
157
158clas_p = { 'a0':0, 'a1':0, 'b0':0, 'b1':0, 'cl_a':0, 'cl_b':0 }
159
160# polynomial parameters
161
162poly_p = { 'n':0, 'degree':0, 'm':0, 'skew':0.0, 'coeff':dict() }
163
164# factorisation parameters
165
166fact_p = { 'n':0, 'dec_digits':0, 'type':'', 'knowndiv':0, 'snfs_difficulty':0,
167 'digs':0, 'qstart':0, 'qstep':0, 'q0':0, 'dd':0, 'primes':list(),
168 'comps':list(), 'divisors':list(),
169 'q_dq':collections.deque([(0,0,0)] * SV_THREADS) }
170
171# Utillity Routines
172
173# print an error message and exit
174
175def die(x, rv = -1):
176 print(x)
177 sys.exit(rv)
178
179def sig_exit(x, y):
180 die('Signal caught. Terminating...')
181
182# obtain a float or an int from a string
183
184def get_nbr(s):
185 try:
186 return float(s)
187 except ValueError:
188 try:
189 return int(s)
190 except ValueError:
191 pass
192 return None
193
194# delete a file (unlink equivalent)
195
196def delete_file(fn):
197 if os.path.exists(fn):
198 try:
199 os.unlink(fn)
200 except WindowsError:
201 pass
202
203# GREP on a list of text lines
204
205def grep_l(pat, lines):
206 r = []
207 c = re.compile(pat)
208 for l in lines:
209 if c.search(l):
210 r += [re.sub('\r|\n', ' ', l)]
211 return r
212
213# GREP on a named file
214
215def grep_f(pat, file_path):
216 r = []
217 c = re.compile(pat)
218 try:
219 with open(file_path, 'r') as in_file:
220 for l in in_file:
221 if c.search(l):
222 r += [re.sub('\r|\n', ' ', l)]
223 return r
224 except IOError:
225 die("can't find file " + file_path)
226
227# concatenate file 'app' to file 'to'
228
229def cat_f(app, to):
230 try:
231 with open(to, 'ab') as out_file:
232 try:
233 with open(app, 'rb') as in_file:
234 if VERBOSE:
235 print('appending {0:s} to {1:s}'.format(app, to))
236 buf = in_file.read(8192)
237 while buf:
238 out_file.write(buf)
239 buf = in_file.read(8192)
240 except IOError:
241 die("can't find file " + app)
242 except IOError:
243 die("can't find file " + to)
244
245# compress file 'fr' to file 'to'
246
247def gzip_f(fr, to):
248 try:
249 with open(fr, 'rb') as in_file:
250 if VERBOSE:
251 print('compressing {0:s} to {1:s}'.format(fr, to))
252 out_file = gzip.open(to, 'ab')
253 out_file.writelines(in_file)
254 out_file.close()
255 except IOError:
256 die("can't find file " + fr)
257
258# remove comment lines
259
260def chomp_comment(s):
261 return re.sub('#.*', '', s).strip()
262
263# produce date/time string for log
264
265def date_time_string() :
266 dt = datetime.datetime.today()
267 return dt.strftime('%a %b %d %H:%M:%S %Y ')
268
269# write string to log(s):
270
271def write_string_to_log(s):
272 with open(CWD + LOGNAME, 'a') as out_f:
273 print(date_time_string() + s, file = out_f)
274
275def output(s, console = True, log = True):
276 if console:
277 print(s)
278 if log:
279 write_string_to_log(s)
280
281# find processor speed
282
283def proc_speed():
284 if os.sys.platform.startswith('win'):
285 if sys.version_info[0] == 2:
286 from _winreg import OpenKey, QueryValueEx, HKEY_LOCAL_MACHINE
287 else:
288 from winreg import OpenKey, QueryValueEx, HKEY_LOCAL_MACHINE
289 h = OpenKey(HKEY_LOCAL_MACHINE,
290 'HARDWARE\\DESCRIPTION\\System\\CentralProcessor\\0')
291 mhz = float(QueryValueEx(h, '~MHz')[0])
292 else:
293 tmp = grep_f('cpu MHz\s+:\s+', '/proc/cpuinfo')
294 m = re.search('\s*cpu MHz\s+:\s+([0-9]+)', tmp[0])
295 mhz = float(m.group(1)) if m else 0.0
296 return 1e-3 * mhz
297
298# check that an executable file exists
299
300def check_binary(exe):
301 if CHECK_BINARIES:
302 pth = MSIEVE_PATH if exe == MSIEVE else GGNFS_PATH
303 pth = os.path.join(pth, exe + EXE_SUFFIX)
304 if not os.path.exists(pth):
305 print('-> Could not find the program: {0:s}.'.format(exe))
306 print('-> Did you set the paths properly in this script?')
307 print('-> They are currently set to:')
308 print('-> GGNFS_BIN_PATH = {0:s}'.format(GGNFS_PATH))
309 print('-> MSIEVE_BIN_PATH = {0:s}'.format(MSIEVE_PATH))
310 sys.exit(-1)
311
312# run an executable file
313
314def run_exe(exe, args, inp = '', in_file = None, out_file = None,
315 log = True, display = VERBOSE, wait = True):
316 al = {} if VERBOSE else {'creationflags' : 0x08000000 }
317 if sys.platform.startswith('win'):
318# priority_high = 0x00000080
319# priority_normal = 0x00000020
320# priority_idle = 0x00000040
321 al['creationflags'] = al.get('creationflags', 0) | 0x00000040
322 else:
323 al['preexec_fn'] = NICE_PATH
324
325 if in_file and os.path.exists(in_file):
326 al['stdin'] = open(in_file, 'r')
327
328 if out_file:
329 if out_file == subprocess.PIPE:
330 md = '> PIPE'
331 al['stdout'] = subprocess.PIPE
332 elif os.path.exists(out_file):
333 md = '>> ' + out_file
334 al['stdout'] = open(out_file, 'a')
335 else:
336 md = '> ' + out_file
337 al['stdout'] = open(out_file, 'w')
338
339 cs = '-> {0:s} {1:s}'.format(exe, args)
340 if in_file:
341 cs += '< {0:s}'.format(in_file)
342 if out_file:
343 cs += md
344
345 output(cs, console = display, log = log)
346
347 ex = os.path.join((GGNFS_PATH if exe != MSIEVE
348 else MSIEVE_PATH), exe + EXE_SUFFIX)
349
350 p = subprocess.Popen([ex] + args.split(' '), **al)
351
352 if not wait:
353 return p
354
355 if out_file == subprocess.PIPE:
356 if sys.version_info[0] == 3:
357 res = p.communicate(input=inp.encode())[0].decode()
358 else:
359 res = p.communicate(input=inp)[0]
360 if res:
361 res = re.split('(?:\r|\n)*', res)
362 ret = p.poll()
363 return (ret, res)
364 else:
365 return p.wait()
366
367def run_msieve(ap, extn='', parallel=False):
368 msd = os.path.abspath(os.path.join(CWD, MSIEVE_PATH))
369 os.chdir(msd)
370 rel = os.path.relpath(CWD)
371 dp = os.path.join(rel, DATNAME + extn)
372 lp = os.path.join(rel, LOGNAME + extn)
373 ip = os.path.join(rel, ININAME + extn)
374 fp = os.path.join(rel, FBNAME + extn)
375 args = ('-s {0:s} -l {1:s} -i {2:s} -nf {3:s} '
376 .format(dp, lp, ip, fp))
377 if parallel:
378 op = os.path.join(rel, OUTNAME + extn)
379 ret = run_exe(MSIEVE, args + ap, out_file=op, wait=False)
380 else:
381 ret = run_exe(MSIEVE, args + ap)
382 os.chdir(CWD)
383 return ret
384
385# generate a list of primes
386
387def prime_list(n):
388 sieve = [False, False] + [True] * (n - 1)
389 for i in range(2, int(n ** 0.5) + 1):
390 if sieve[i]:
391 m = n // i - i
392 sieve[i * i : n + 1 : i] = [False] * (m + 1)
393 return [i for i in range(n + 1) if sieve[i]]
394
395# greatest common divisor
396
397def gcd(x, y):
398 if x == 0:
399 return y
400 elif y == 0:
401 return x
402 else:
403 return gcd(y % x, x)
404
405# Miller Rabin 'probable prime' test
406#
407# returns 'False' if 'n' is definitely composite
408# returns 'True' if 'n' is prime or probably prime
409#
410# 'r' is the number of trials performed
411
412def miller_rabin(n, r = 10):
413 t = n - 1
414 s = 0
415 while not t & 1:
416 t >>= 1
417 s += 1
418 for i in range(r):
419 a = random.randint(2, n - 1)
420 x = pow(a, t, n)
421 if x != 1 and x != n - 1:
422 for j in range(s - 1):
423 x = (x * x) % n
424 if x == 1:
425 return False
426 if x == n - 1:
427 break
428 else:
429 return False
430 return True
431
432# determine if n is probably prime - return
433# 0 if n is 0, 1 or composite
434# 1 if n is probably prime
435# 2 if n is definitely prime
436
437def probable_prime_p(nn, r):
438 n = abs(nn)
439 if n <= 2:
440 return 2 if n == 2 else 0
441
442 # trial division
443 for p in prime_list(1000):
444 if not n % p:
445 return 2 if n == p else 0
446 if p * p > n:
447 return 2
448
449 # Fermat test
450 if pow(random.randint(2, n - 1), n - 1, n) != 1:
451 return 0
452
453 # Miller-Rabin test
454 return 1 if miller_rabin(n, r) else 0
455
456# count the number of lines in a file
457
458def linecount(file_path):
459 count = 0
460 if not os.path.exists(file_path):
461 die('can\'t open {0:s}'.format(file_path))
462 with open(file_path, 'r') as in_file:
463 for l in in_file:
464 count += 1
465 return count
466
467# Read the log file to see if we've found all the prime
468# divisors of N yet.
469
470def get_primes(fact_p):
471
472 with open(LOGNAME, 'r') as in_f:
473 for l in in_f:
474 m1 = re.search('r\d+=(\d+)\s+', l.rstrip())
475 if m1:
476 val = int(m1.group(1))
477 if len(val) > 1 and len(val) < len(fact_p['ndivfree']):
478 # Is this a prime divisor or composite?
479 m2 = re.search('\(pp(\d+)\)', l)
480 if m2:
481 # If this is a prime we don't already have, add it.
482 found = False
483 for p in fact_p['primes']:
484 if val == p:
485 found = True
486 if not found:
487 fact_p['primes'].append(val)
488 else:
489 fact_p['comps'].append(val)
490
491 # Now, try to figure out if we have all the prime factors:
492 x = itertools.reduce(lambda x,y : x * y, fact_p['primes'], 1)
493 if x == fact_p['ndivfree'] or probab_prime_p(fact_p['ndivfree'] // x, 10):
494 if x != fact_p['ndivfree']:
495 fact_p['primes'].append(fact_p['ndivfree'] // x)
496 for p in fact_p['primes']:
497 cs = '-> p: {0:s} (pp{1:d})'.format(val, len(val))
498 output(cs)
499 return True
500 # Here, we could try to recover other factors by division,
501 # but until we have a primality test available, this would
502 # be pointless since we couldn't really know if we're done.
503 return False
504
505# The parameter degree, if nonzero, will let this function adjust
506# parameters for SNFS factorizations with a predetermined polynomial
507# of degree that is not the optimal degree.
508
509nonPrefDegAdjust = 12
510
511def load_default_parameters(nfs_type, digits, degree,
512 fact_p, pols_p, lats_p, clas_p):
513
514 if nfs_type == 'gnfs':
515 lats_p['lss'] = 0
516 pth = os.path.join(GGNFS_PATH, DEFAULT_PAR_FILE)
517 if not os.path.exists(pth):
518 die('Could not find default parameter file {0:s}!'
519 .format(pth))
520 with open(pth, 'r') as in_f:
521 par_digits = 0
522 par_degree = 0
523 how_close = 1000
524 for l in in_f:
525 l = chomp_comment(l.strip()).strip()
526 if l:
527 tu = l.split(',')
528 t = tu[0]
529 if t == nfs_type:
530 o = 2
531 cand_digits = int(tu[1])
532 cand_degree = int(tu[o + 0])
533 s1 = abs(cand_digits - digits)
534 s2 = abs(cand_digits - nonPrefDegAdjust - digits) + nonPrefDegAdjust
535
536 # try to properly handle crossover from degree 4 to degree 5
537 if (s1 if nfs_type == 'gnfs' or not degree or degree == cand_degree
538 or par_degree == cand_degree - 1 else s2) < how_close:
539
540 how_close = (s2 if nfs_type != 'gnfs' and degree
541 and cand_degree != degree else s1)
542 par_digits = cand_digits
543 par_degree = cand_degree
544
545 pols_p['maxs1'] = int(tu[o + 1])
546 pols_p['maxskew'] = int(tu[o + 2])
547 pols_p['goodscore'] = float(tu[o + 3])
548 pols_p['examinefrac'] = float(tu[o + 4])
549 pols_p['j0'] = int(tu[o + 5])
550 pols_p['j1'] = int(tu[o + 6])
551 pols_p['estepsize'] = int(tu[o + 7])
552 pols_p['maxtime'] = int(tu[o + 8])
553
554 lats_p['rlim'] = int(tu[o + 9])
555 lats_p['alim'] = int(tu[o + 10])
556 lats_p['lpbr'] = int(tu[o + 11])
557 lats_p['lpba'] = int(tu[o + 12])
558 lats_p['mfbr'] = int(tu[o + 13])
559 lats_p['mfba'] = int(tu[o + 14])
560 lats_p['rlambda'] = float(tu[o + 15])
561 lats_p['alambda'] = float(tu[o + 16])
562 lats_p['qintsize'] = int(tu[o + 17])
563
564 clas_p['cl_a'] = int(tu[o + 18])
565 clas_p['cl_b'] = int(tu[o + 19])
566
567 fact_p['qstep'] = lats_p['qintsize']
568
569 fact_p['digs'] = digits / 0.72 if nfs_type == 'gnfs' else digits
570 # 0.72 is inspired by T.Womack's crossover 28/29 for GNFS-144 among other consideration
571
572 if fact_p['digs'] >= 160:
573 # the table parameters are easily splined; the table may be not needed at all --SB.
574 lats_p['rlim'] = lats_p['alim'] = int(0.07 * 10 ** (fact_p['digs'] / 60.0) + 0.5) * 100000
575 lats_p['lpbr'] = lats_p['lpba'] = int(21 + fact_p['digs'] / 25.0)
576 lats_p['mfbr'] = lats_p['mfba'] = 2 * lats_p['lpbr'] - (1 if fact_p['digs'] < 190 else 0)
577 lats_p['rlambda'] = lats_p['alambda'] = 2.5 if fact_p['digs'] < 200 else 2.6
578 lats_p['qintsize'] = fact_p['qstep'] = 100000
579 clas_p['cl_a'] = 4000000
580 clas_p['cl_b'] = 400
581 par_digits = digits
582
583 pols_p['degree'] = par_degree
584 print('-> Selected default factorization parameters for {0:d} digit level.'.
585 format(par_digits))
586 if nfs_type == 'gnfs':
587 r = (95, 110, 140, 158, 185, 999)
588 else:
589 if degree and par_degree != degree:
590 digits += nonPrefDegAdjust
591 r = (120, 150, 195, 240, 275, 999)
592
593 i = 1
594 for v in r:
595 if digits < v:
596 break
597 i += 1
598 else:
599 die('You are joking?')
600
601 lats_p['siever'] = 'gnfs-lasieve4I1' + str(i) + 'e'
602 output('-> Selected lattice siever: {0:s}'.format(lats_p['siever']))
603
604# These are default parameters for polynomial selection using the
605# Kleinjung/Franke tool.
606# arg0 = number of digits in N.
607
608def load_pol5_parameters(digits, pol5_p):
609
610 par_digits = 0
611 if NO_DEF_NM_PARAM or digits < 100:
612 pol5_p['search_a5step'] = 1
613 pol5_p['npr'] = int(digits / 13.0 - 4.5)
614 pol5_p['npr'] = pol5_p['npr'] if pol5_p['npr'] > 4 else 4
615 pol5_p['normmax'] = 10 ** (0.163 * digits - 1.4794)
616 pol5_p['normmax1'] = 10 ** (0.1522 * digits - 1.6969)
617 pol5_p['normmax2'] = 10 ** (0.142 * digits - 2.6429)
618 pol5_p['murphymax'] = 10 ** (-0.0569 * digits - 2.8452)
619 pol5_p['max_pst_time'] = int(0.000004 / pol5_p['murphymax'])
620 output('-> Selected polsel parameters for {0:d} digit level.'.format(digits))
621 return
622
623 pth = os.path.join(GGNFS_PATH, DEFAULT_POLSEL_PAR_FILE)
624 if not os.path.exists(pth):
625 die('Could not find default parameter file {0:s}!'
626 .format(pth))
627
628 with open(pth, 'r') as in_f:
629 how_close = 1000
630 for l in in_f:
631 l = chomp_comment(l.strip()).strip()
632 if l:
633 tu = l.split(',')
634 d = int(tu[0])
635 if abs(d - digits) < how_close:
636 o = 1
637 how_close = abs(d - digits)
638 par_digits = d
639 pol5_p['max_pst_time'] = 60 * int(tu[o + 0])
640 pol5_p['search_a5step'] = int(tu[o + 1])
641 pol5_p['npr'] = int(tu[o + 2])
642 pol5_p['normmax'] = float(tu[o + 3])
643 pol5_p['normmax1'] = float(tu[o + 4])
644 pol5_p['normmax2'] = float(tu[o + 5])
645 pol5_p['murphymax'] = float(tu[o + 6])
646
647 output('-> Selected default polsel parameters for {0:d} digit level.'
648 .format(par_digits))
649
650# The dictionary entries in this routine map each coefficient as a string
651# to the coefficient value as floating point values. The (name, value)
652# pairs on the BEGIN_POLY line are also entered into the dictionary for
653# each polynomial
654
655pname = ''
656
657def terminate_pol5(x, y):
658 die('Terminated on {0:s} by SIGINT'.format(pname))
659
660def run_pol5(fact_p, pol5_p, lats_p, clas_p):
661 global pname
662
663 check_binary(POL51M0)
664 check_binary(POL51OPT)
665 projectname = NAME + '.polsel'
666 host = socket.gethostname()
667 pname = projectname + host + '.' + str(os.getpid())
668
669 with open(pname + '.data', 'w') as out_f:
670 print('N ' + str(fact_p['n']), file = out_f)
671
672 pol5_term_flag = False
673 old_handler = signal.signal(signal.SIGINT, terminate_pol5)
674
675 load_pol5_parameters(fact_p['dec_digits'], pol5_p)
676 pol5_p['max_pst_time'] *= polySelTimeMultiplier
677 hmult = 1e3
678 load_default_parameters('gnfs', fact_p['dec_digits'], 5,
679 fact_p, pols_p, lats_p, clas_p)
680 bestpolyinf = dict()
681 bestpolyinf['Murphy_E'] = 0
682
683 start_time = time.time()
684 nerr = 0
685 h_lo = 0
686 while not pol5_term_flag and nerr < 2:
687 h_hi = h_lo + pol5_p['search_a5step']
688
689 output('-> Searching leading coefficients from {0:g} to {1:g}'
690 .format(h_lo * hmult + 1, h_hi * hmult))
691 args = ('-b {0:s} -v -v -p {1:d} -n {2:g} -a {3:g} -A {4:g}'
692 .format(pname, pol5_p['npr'], pol5_p['normmax'], h_lo, h_hi))
693
694 if not run_exe(POL51M0, args, out_file = pname + '.log'):
695
696 # lambda-comp related errors can be skipped and some polys are
697 # then found ull5-comp related are probably fatal, but let the
698 # elapsed time.time() take care of them
699
700 nerr = 0
701 suc = grep_f('success', pname + '.log')
702 changed = False
703 if suc:
704 args = ('-b {0:s} -v -v -n {1:g} -N {2:g} -e {3:g}'
705 .format(pname, pol5_p['normmax1'], pol5_p['normmax2'], pol5_p['murphymax']))
706 ret = run_exe(POL51OPT, args, out_file = pname + '.log')
707 if ret:
708 die('Abnormal return value {0:d}. Terminating...'.format(ret))
709 cat_f(pname + '.51.m', projectname + '.51.m.all')
710
711 with open(pname + '.cand', 'r') as in_f:
712
713 polyinf = dict()
714 for l in in_f:
715 l = l.rstrip()
716 if re.match('BEGIN POLY', l):
717 l = re.sub('BEGIN POLY', '', l)
718 l = re.sub('^ #', '', l)
719 tu = l.split()
720 for i in range(0, len(tu), 2):
721 polyinf[tu[i]] = float(tu[i + 1])
722 elif re.match('END POLY', l):
723 if polyinf['Murphy_E'] > bestpolyinf['Murphy_E']:
724 bestpolyinf = polyinf.copy()
725 changed = True
726 polyinf = dict()
727 else:
728 tu = l.split()
729 polyinf[re.sub('^X', 'c', tu[0])] = int(tu[1])
730
731 if changed:
732
733 for key in sorted(bestpolyinf):
734 print(key, ':', bestpolyinf[key])
735
736 with open(NAME + '.poly', 'w') as out_f:
737
738 print('name: {0:s}'.format(NAME), file = out_f)
739 print('n: {0:d}'.format(fact_p['n']), file = out_f)
740
741 for key in reversed(sorted(bestpolyinf)):
742 if re.match('c\d+', key) or re.match('Y\d+', key):
743 print('{0:s}: {1:d}'
744 .format(key, bestpolyinf[key]), file = out_f)
745 elif re.match('skewness', key):
746 print('skew: {0:<10.2f}'
747 .format(bestpolyinf[key]), file = out_f)
748 elif key == 'M':
749 print(key, bestpolyinf[key])
750 print('# {0:s} {1:d}'
751 .format(key, bestpolyinf[key]), file = out_f)
752 else:
753 print('# {0:s} {1:g}'
754 .format(key, bestpolyinf[key]), file = out_f)
755
756 print('type: {0:s}'.format('gnfs'), file = out_f)
757 print('rlim: {0:d}'.format(lats_p['rlim']), file = out_f)
758 print('alim: {0:d}'.format(lats_p['alim']), file = out_f)
759 print('lpbr: {0:d}'.format(lats_p['lpbr']), file = out_f)
760 print('lpba: {0:d}'.format(lats_p['lpba']), file = out_f)
761 print('mfbr: {0:d}'.format(lats_p['mfbr']), file = out_f)
762 print('mfba: {0:d}'.format(lats_p['mfba']), file = out_f)
763 print('rlambda: {0:g}'.format(lats_p['rlambda']), file = out_f)
764 print('alambda: {0:g}'.format(lats_p['alambda']), file = out_f)
765 print('qintsize: {0:d}'.format(lats_p['qintsize']), file = out_f)
766
767 cat_f(pname +'.cand', projectname + '.cand.all')
768 delete_file(pname + '.cand')
769
770 print('-> =====================================================')
771 output('-> Best score so far: {0:g} (good_score={1:g})'
772 .format(bestpolyinf['Murphy_E'], pol5_p['murphymax']))
773 print('-> =====================================================')
774 delete_file(pname + '.log')
775 delete_file(pname + '.51.m')
776
777 if time.time() > start_time + pol5_p['max_pst_time']:
778 pol5_term_flag = True
779 h_lo = h_hi
780
781 delete_file(pname + '.data')
782 signal.signal(signal.SIGINT, old_handler)
783
784# We will start with a higher leading coefficient divisor. When it
785# appears that we are searching in an interesting range, it will
786# be backed down so that the resulting range can be searched with
787# a finer resolution. This means that from time to time, the same
788# poly will be found several times as we hone in on a region.
789
790def run_poly_select(fact_p, pols_p, lats_p, clas_p):
791
792 check_binary(POLYSELECT)
793 lcd_choices = (2, 4, 4, 12, 12, 24, 24, 48, 48, 144, 144, 720, 5040)
794 lcd_level = 4 + (fact_p['dec_digits'] - 70) // 10
795 if lcd_level < 0:
796 lcd_level = 0
797 if lcd_level > 12:
798 lcd_level = 12
799
800 output('-> Starting search with leading coefficient divisor {0:d}'
801 .format(lcd_choices[lcd_level]))
802
803 load_default_parameters('gnfs', fact_p['dec_digits'], 0,
804 fact_p, pols_p, lats_p, clas_p)
805
806 pols_p['maxtime'] *= polySelTimeMultiplier
807 best_score = 0.0
808 start_time = time.time()
809 done = False
810 best_lc = 1
811 last_lc = 0
812 multiplier = 0.75
813 e_lo = 1
814
815 while not done:
816 e_hi = e_lo + pols_p['estepsize']
817
818 lcd = lcd_choices[lcd_level]
819 with open(NAME + '.polsel', 'w') as out_f:
820
821 print('name: {0:s}'.format(NAME), file = out_f)
822 print('n: {0:d}'.format(fact_p['n']), file = out_f)
823 print('deg: {0:d}'.format(pols_p['degree']), file = out_f)
824 print('bf: {0:s}.best.poly'.format(NAME), file = out_f)
825 print('maxs1: {0:d}'.format(pols_p['maxs1']), file = out_f)
826 print('maxskew: {0:d}'.format(pols_p['maxskew']), file = out_f)
827 print('enum: {0:d}'.format(lcd), file = out_f)
828 print('e0: {0:d}'.format(e_lo), file = out_f)
829 print('e1: {0:d}'.format(e_hi), file = out_f)
830 print('cutoff: {0:g}'.format(0.75 * pols_p['goodscore']), file = out_f)
831 print('examinefrac: {0:g}'.format(pols_p['examinefrac']), file = out_f)
832 print('j0: {0:d}'.format(pols_p['j0']), file = out_f)
833 print('j1: {0:d}'.format(pols_p['j1']), file = out_f)
834
835 args = '-if {0:s}.polsel'.format(NAME)
836 ret = run_exe(POLYSELECT, args)
837 if ret:
838 die('Return value res. Terminating...'.format(ret))
839
840 # Find the score of the best polynomial: E(F1,F2) =
841 inp = True
842 try:
843 tmp = grep_f('E\(F1,F2\) =', NAME + '.best.poly')
844 except IOError:
845 inp = False
846
847 if inp:
848 score = tmp[0]
849 score = float(re.sub('.*E\(F1,F2\) =', '', tmp[0]))
850 tmp = grep_f('^c' + str(pols_p['degree']), NAME + '.best.poly')
851 m = re.search('^c' + str(pols_p['degree']) + ':\s([\+\-]*\d*)', tmp[0])
852 lc = int(m.group(1))
853 else:
854 score = 0.0
855 if (score > multiplier * pols_p['goodscore'] and lc != best_lc
856 and lc != last_lc ):
857 multipler = 1.1 * multiplier if multiplier < 0.9 / 1.1 else multiplier
858 last_lc = lc
859 new_lcd_level = lcd_level - 1 if lcd > 0 else 0
860 e_lo = e_lo * lcd_choices[lcd_level] // lcd_choices[new_lcd_level] - pols_p['estepsize']
861
862 if lcd_level != new_lcd_level:
863 print('-> Leading coefficient divisor dropped from {0:d} to {1:d}.'
864 .format(lcd_choices[lcd_level], lcd_choices[new_lcd_level]))
865 lcd_level = new_lcd_level
866
867 if score > best_score and lc != best_lc:
868 best_score = score
869 best_lc = lc
870 last_best_time = time.time()
871 os.rename(NAME + '.best.poly', NAME + '.thebest.poly')
872
873 # We should now fill in the missing parameters with the default
874 # loaded in from table. Do this now, so that the user has the
875 # option to kill the script and still have a viable poly file.
876
877 with open(NAME + '.thebest.poly', 'r') as in_f:
878 with open(NAME + '.poly', 'w') as out_f:
879
880 for l in in_f:
881 l = l.rstrip()
882 if l.find('rlim:') != -1:
883 l = 'rlim: {0:d}'.format(lats_p['rlim'])
884 elif l.find('alim:') != -1:
885 l = 'alim: {0:d}'.format(lats_p['alim'])
886 elif l.find('lpbr:') != -1:
887 l = 'lpbr: {0:d}'.format(lats_p['lpbr'])
888 elif l.find('lpba:') != -1:
889 l = 'lpba: {0:d}'.format(lats_p['lpba'])
890 elif l.find('mfbr:') != -1:
891 l = 'mfbr: {0:d}'.format(lats_p['mfbr'])
892 elif l.find('mfba:') != -1:
893 l = 'mfba: {0:d}'.format(lats_p['mfba'])
894 elif l.find('rlambda:') != -1:
895 l = 'rlambda: {0:g}'.format(lats_p['rlambda'])
896 elif l.find('alambda:') != -1:
897 l = 'alambda: {0:g}'.format(lats_p['alambda'])
898 elif l.find('qintsize:') != -1:
899 l = 'qintsize: {0:d}'.format(lats_p['qintsize'])
900 if l:
901 print(l, file = out_f)
902 print('type: gnfs', file = out_f)
903
904 delete_file(NAME + '.thebest.poly')
905 delete_file(NAME + '.best.poly')
906
907 print('-> =====================================================')
908 output('-> Best score so far: {0:g} (good_score={1:g}) '
909 .format(best_score, pols_p['goodscore']))
910 print('-> =====================================================')
911
912 # We will allow another 5 minutes just in case there happens to
913 # be a really good poly nearby (or the 'good_score' value was too low)
914 if best_score > 1.4 * pols_p['goodscore']:
915 done = (time.time() > last_best_time + 300)
916
917 done |= (time.time() > start_time + pols_p['maxtime'])
918 e_lo = e_hi
919
920 output('-> Using poly with score={0:f}'.format(best_score))
921 delete_file(NAME + '.polsel')
922
923def msieve_mpqs(fact_p):
924
925 with open(ININAME, 'w') as out_f:
926 out_f.write('{0:d}'.format(fact_p['n']))
927 ret = run_msieve('-v')
928 if ret:
929 die('Msieve Error: return value {0:d}... terminating...'.format(ret))
930 return True
931
932procs = []
933
934def terminate_threads(sig=signal.SIGINT, frame=None):
935 global procs
936 try:
937 for p in procs:
938 if p.poll() == None:
939 p.terminate()
940 for p in procs:
941 if p.poll() == None:
942 p.wait()
943 except:
944 pass
945
946def fb_to_poly():
947 with open(NAME + '.poly', 'w') as out_f:
948 with open(NAME + '.fb', 'r') as in_f:
949 for l in in_f:
950
951 m = re.match('N\s+(\d+)', l)
952 if m:
953 print('n: {0:s}'.format(m.group(1)), file = out_f)
954
955 m = re.match('SKEW\s+(\d*\.\d*)', l)
956 if m:
957 skew = float(m.group(1))
958
959 m = re.match('R(\d+)\s+([+-]?\d+)', l)
960 if m:
961 print('Y{0:s}: {1:s}'.format(m.group(1),
962 m.group(2)), file = out_f)
963
964 m = re.match('A(\d+)\s+([+-]?\d+)', l)
965 if m:
966 print('c{0:s}: {1:s}'.format(m.group(1),
967 m.group(2)), file = out_f)
968 try:
969 print('skew: {0:<10.2f}'.format(skew), file = out_f)
970 print('type: gnfs', file = out_f)
971 except:
972 die('{0:s}.fb is not in the correct format'.format(NAME))
973
974def find_best(fn, bestscore, poly):
975 try:
976 in_f = open(fn, 'r')
977 except IOError:
978 die("can't read " + fn)
979
980 # first input line (followed by polynomial lines) format:
981 # norm 1.498232e-10 alpha -5.818790 e 9.344e-10 rroots 3
982 # 1 2 3 4 5 6 7 8
983 for l in in_f:
984 word = l.split()
985 if word[1] == 'norm':
986 better = False
987 if float(word[6]) > bestscore:
988 bestscore = float(word[6])
989 output("Best score so far: {0:s}".format(l.strip()))
990 better = True
991 poly = []
992 if better:
993 poly += [l]
994 in_f.close()
995 return (bestscore, poly)
996
997def run_msieve_poly(fact_p, minv, maxv):
998 global procs
999
1000 poly = []
1001 bestscore = 0
1002 if USE_CUDA or MS_THREADS == 1:
1003 with open(ININAME, 'w') as out_f:
1004 out_f.write('{0:d}'.format(fact_p['n']))
1005 if USE_CUDA:
1006 ap = '-g {0:d} -v -np'.format(GPU_NUM)
1007 bp = ' Is CUDA enabled?'
1008 else:
1009 ap = '-v -np'
1010 bp = ''
1011 if MSIEVE_POLY_TIME_LIMIT:
1012 ap = '-d {0:d} '.format(MSIEVE_POLY_TIME_LIMIT) + ap
1013 else:
1014 ap += ' {0:d},{1:d} -t {2:d}'.format(minv if minv else 1, maxv, NUM_CORES)
1015 ret = run_msieve(ap)
1016 if ret:
1017 die('Msieve Error: return value {0:d}.{1:s} Terminating...'.format(ret, bp))
1018 bestscore, poly = find_best(NAME + '.msieve.dat.p', bestscore, poly)
1019
1020 else:
1021 old_handler = signal.signal(signal.SIGINT, terminate_threads)
1022 for j in range(MS_THREADS):
1023 extn = '.T' + str(j)
1024 with open(ININAME + extn, 'w') as out_f:
1025 out_f.write('{0:d}'.format(fact_p['n']))
1026
1027 procs = [None] * MS_THREADS
1028 ret = 0
1029 base = minv
1030 if os.path.exists(NAME + '.lcf'):
1031 try:
1032 with open(NAME + '.lcf', 'r') as in_f:
1033 t = int(in_f.readline().strip())
1034 if minv < t < maxv:
1035 base = t
1036 elif t >= maxv:
1037 die("previous run conflict (delete {0:s}.lcf to run again)".format(NAME))
1038 except IOError:
1039 die("Can't read " + NAME + '.lcf')
1040 term = False
1041 while base < maxv and not ret:
1042 for j in range(MS_THREADS):
1043 extn = '.T' + str(j)
1044 retc = 0 if not procs[j] else procs[j].poll()
1045 if retc != None:
1046 ret |= retc
1047 if os.path.exists(NAME + '.lcf'):
1048 try:
1049 with open(NAME + '.lcf', 'r') as in_f:
1050 base = int(in_f.readline().strip())
1051 except IOError:
1052 die("Can't read " + NAME + '.lcf')
1053 try:
1054 with open(NAME + '.lcf', 'w') as out_f:
1055 print(base + 100, file=out_f)
1056 except IOError:
1057 die("Can't write to " + NAME + '.lcf')
1058 procs[j] = run_msieve('-np {0:d},{1:d} -v '.format(base + 1, base + 100), extn, parallel=True)
1059 base += 100
1060 try:
1061 time.sleep(1)
1062 except:
1063 term = True
1064 break
1065 if ret:
1066 die('an error occurred during polynomial search')
1067 if term:
1068 terminate_threads()
1069 die('msieve terminated')
1070
1071 for j in range(MS_THREADS):
1072 nm = NAME + '.fb.T' + str(j)
1073 if os.path.exists(nm):
1074 delete_file(nm)
1075 output('deleted ' + nm)
1076 nm = NAME + '.dat.T' + str(j) + '.p'
1077 bestscore, poly = find_best(NAME + '.dat.T' + str(j) + '.p', bestscore, poly)
1078
1079 if bestscore == 0:
1080 output('could not find any polynomials')
1081 die('could not find any polynomials')
1082
1083 with open(NAME + '.poly', 'w') as out_f:
1084 with open(NAME + '.n', 'r') as in_f:
1085 for l in in_f:
1086 out_f.writelines([l])
1087 out_f.writelines(poly + ['type: gnfs'])
1088 out_f.close()
1089
1090 return True
1091
1092def change_parameters(lats_p):
1093
1094 check_binary(MAKEFB)
1095 check_binary(PROCRELS)
1096 cs = ('-> Parameter change detected, dumping relations... ')
1097 output(cs)
1098
1099 args = '-fb {0:s}.fb -prel {1:s} -dump'.format(NAME, RELSBIN)
1100 ret = run_exe(PROCRELS, args)
1101 if ret:
1102 die('Return value {0:d}. Terminating...'.format(ret))
1103
1104 for f in glob.iglob(RELSBIN + '*'):
1105 delete_file(f)
1106 for f in glob.iglob('cols*'):
1107 delete_file(f)
1108 for f in glob.iglob('deps*'):
1109 delete_file(f)
1110 delete_file('factor.easy')
1111 for f in glob.iglob('lpindex*'):
1112 delete_file(f)
1113 for f in glob.iglob(NAME + '.*.afb.*'):
1114 delete_file(f)
1115
1116 output('-> Making new factor base files...')
1117 args = ('-rl {0[rlim]:d} -al {0[alim]:d} -lpbr {0[lpbr]:d} -lpba {0[lpba]:d}'
1118 '{1:s} -of {2:s}.fb -if {2:s}.poly'.format(lats_p, LARGEPRIMES, NAME))
1119 ret = run_exe(MAKEFB, args)
1120 if ret:
1121 die('Return value {0:d}. Terminating...'.format(ret))
1122
1123 output('-> Reprocessing siever output...')
1124 i = 0
1125 while os.path.exists('spairs.dump.i'):
1126 args = (' -fb {0:s}.fb -prel {1:s} -newrel spairs.dump.i -nolpcount'
1127 .format(NAME, RELSBIN))
1128 ret = run_exe(PROCRELS, args)
1129 i = i + 1
1130
1131 # Update the paramfile, used by this script to detect change
1132 # in parameter settings.
1133 with open(PARAMFILE, 'w') as out_f:
1134 print('rlim: {0:d}'.format(lats_p['rlim']), file = out_f)
1135 print('alim: {0:d}'.format(lats_p['alim']), file = out_f)
1136 print('lbpr: {0:s}'.format(LBPR), file = out_f)
1137 print('lbpa: {0:s}'.format(LBPA), file = out_f)
1138
1139def check_for_parameter_change(lats_p):
1140
1141 if not os.path.exists(PARAMFILE):
1142
1143 output('-> Creating param file to detect parameter changes...')
1144 with open(PARAMFILE, 'w') as out_f:
1145 print('rlim: {0:d}'.format(lats_p['rlim']), file = out_f)
1146 print('alim: {0:d}'.format(lats_p['alim']), file = out_f)
1147 print('lpbr: {0:d}'.format(lats_p['lpbr']), file = out_f)
1148 print('lpba: {0:d}'.format(lats_p['lpba']), file = out_f)
1149 return
1150
1151 # Okay - it exists. We should check it to see if the
1152 # parameters are the same as the current ones.
1153 with open(PARAMFILE, 'r') as in_f:
1154 old_params = in_f.readlines()
1155
1156 tmp = grep_l('rlim:', old_params)
1157 old_rlim = int(re.sub('.*rlim:', '', tmp[0]))
1158 tmp = grep_l('alim:', old_params)
1159 old_alim = int(re.sub('.*alim:', '', tmp[0]))
1160 tmp = grep_l('lpbr:', old_params)
1161 old_lpbr = int(re.sub('.*lpbr:', '', tmp[0]))
1162 tmp = grep_l('lpba:', old_params)
1163 old_lpba = int(re.sub('.*lpba:', '', tmp[0]))
1164 if (old_rlim != lats_p['rlim'] or old_alim != lats_p['alim'] or
1165 old_lpbr != lats_p['lpbr'] or old_lpba != lats_p['lpba']):
1166 change_parameters(lats_p)
1167 else:
1168 output('-> No parameter change detected, resuming... ')
1169
1170def plot_lp():
1171
1172 if CHECK_BINARIES:
1173 if not os.path.exists(PLOT):
1174 return
1175
1176 if not os.path.exists(LOGNAME):
1177 return
1178
1179 with open(LOGNAME, 'r') as in_f:
1180 with open('.lprels', 'w') as out_f:
1181 print('0, 0', file = out_f)
1182 with open('.rels', 'w') as out_rels:
1183 print('0, 0', file = out_rels)
1184
1185 for l in in_f:
1186 l = l.rstrip().upper()
1187 m = re.search('LARGEPRIMES')
1188 if m:
1189 l = re.sub('\[.*]\s*','', l) # strip date
1190 l = re.sub('(LARGEPRIMES: |RELATIONS: )', '', l) # strip the labels
1191 tu = l.strip().split(',')
1192 print('{0:s}, {1:d}'.format(tu[1], int(tu[0]) - int(tu[1])), FILE = outf)
1193
1194 m = re.search('FINALFF')
1195 if m:
1196 l = re.sub('\[.*]\s*','', l) # strip date
1197 l = re.sub('(LARGEPRIMES: |RELATIONS: )', '', l) # strip the labels
1198 tu = l.strip().split(',')
1199 print('{0:s}, {1:s}'.format(tu[0], tu[1]), file = out_rels)
1200
1201 os.putenv('XAXIS', 'Total relations')
1202 os.putenv('YAXIS', '')
1203 args = 'xprimes.png \'ExcessLargePrimes\' .lprels'
1204 ret = run_exe(PLOT, args)
1205 if ret:
1206 die('Return value {0:d}. Terminating...'.format(ret))
1207
1208 os.putenv('XAXIS', 'Total relations')
1209 os.putenv('YAXIS', 'Full relation-sets')
1210 args = 'relations.png \'TotalFF\' .rels'
1211 ret = run_exe(PLOT, args)
1212 if ret:
1213 die('Return value {0:d}. Terminating...'.format(ret))
1214
1215 delete_file('.lprels')
1216 delete_file('.rels')
1217
1218class is_missing(Exception):
1219 def __init__(self, x):
1220 self.value = x
1221 def __str__(self):
1222 return self.value
1223
1224def check_parameters(fact_p, poly_p, lats_p):
1225
1226 check_binary(MSIEVE)
1227 check_binary(lats_p['siever'])
1228 if CHECK_BINARIES:
1229 delete_file('.lprels')
1230 delete_file('.rels')
1231 with open('.rels', 'w') as out_f:
1232 print('{0:d}'.format(fact_p['n']), file = out_f)
1233 delete_file('.lprels')
1234 delete_file('.rels')
1235
1236 if not fact_p['n']:
1237 die('Error: \'n\' not supplied!')
1238 if not (poly_p['m'] or 'Y1' in poly_p['coeff']):
1239 die('Error: \'m\' not supplied!')
1240 if not 'c' + str(poly_p['degree']) in poly_p['coeff']:
1241 die('Error: polynomial not supplied!')
1242 if not poly_p['skew']:
1243 poly_p['skew'] = (abs( poly_p['coeff']['c0'] / poly_p['coeff']['c'
1244 + str(poly_p['degree'])]) ** (1.0 / poly_p['degree']))
1245 print('-> Using calculated skew {0:<10.2f}'.format(poly_p['skew']))
1246
1247 try :
1248 if not lats_p['rlim']: raise is_missing('rlim')
1249 if not lats_p['alim']: raise is_missing('alim')
1250 if not lats_p['lpbr']: raise is_missing('lpbr')
1251 if not lats_p['lpba']: raise is_missing('lpba')
1252 if not lats_p['mfbr']: raise is_missing('mfbr')
1253 if not lats_p['mfba']: raise is_missing('mfba')
1254 if not lats_p['rlambda']: raise is_missing('rlambda')
1255 if not lats_p['alambda']: raise is_missing('alambda')
1256 if not fact_p['qstep']: raise is_missing('qstep')
1257 except is_missing as x:
1258 die('Error: \'' + str(x) + '\' not supplied!')
1259
1260def get_parm_int(data, match):
1261 tmp = grep_l('^' + match + ':', data)
1262 if tmp:
1263 m = re.search('^' + match + ':\s*(-?\d+)', tmp[0])
1264 if m:
1265 return int(m.group(1))
1266 return None
1267
1268# Read default parameters first. Then we will just override
1269# by any user-supplied parameters.
1270
1271def read_parameters(fact_p, poly_p, lats_p):
1272
1273 with open(NAME + '.poly', 'r') as in_f:
1274 file_lines = in_f.readlines()
1275
1276 fact_p['n'] = get_parm_int(file_lines, 'n')
1277 fact_p['dec_digits'] = len(str(fact_p['n']))
1278
1279 # Find the polynomial degree.
1280 coefvals = dict()
1281 coeff_lines = grep_l('^c\d+:', file_lines)
1282 d = 0
1283 for l in coeff_lines:
1284 # Grab the coefficient index
1285 # First char of line 'c' followed by a digit string.
1286 m = re.search('^c(\d+):\s*(-?\d+)', l)
1287 if m:
1288 key = int(m.group(1))
1289 val = int(m.group(2))
1290 if key > d:
1291 d = key
1292 if key in coefvals:
1293 print('-> Warning: redefining c{0:d}'.format(key))
1294 coefvals[key] = val
1295
1296 poly_p['degree'] = get_parm_int(file_lines, 'deg')
1297 if not poly_p['degree']:
1298 poly_p['degree'] = d
1299 if poly_p['degree'] != d:
1300 print('-> Error: poly file specifies degree {0:d} but highest poly'
1301 .format(poly_p['degree']))
1302 print('-> coefficient given is c{0:d}!'.format(d))
1303 sys.exit(-1)
1304
1305 for i in range(d):
1306 if i not in coefvals:
1307 coefvals[i] = 0
1308
1309 commonfac = 0
1310 first = True
1311 for key in reversed(sorted(coefvals)):
1312 if first:
1313 commonfac = coefvals[key]
1314 first = False
1315 else:
1316 commonfac = gcd(commonfac, coefvals[key])
1317
1318 if CHECK_POLY and commonfac > 1:
1319 print('-> Error: poly coefficients have a common factor commonfac.'
1320 ' Please divide it out.')
1321 sys.exit(-1)
1322
1323 denom = get_parm_int(file_lines, 'Y1')
1324 numer = get_parm_int(file_lines, 'Y0')
1325 poly_p['m'] = get_parm_int(file_lines, 'm')
1326
1327 if denom and numer:
1328 if denom < 0:
1329 denom = -denom
1330 else:
1331 numer = -numer
1332 # print '-> Common root is numer / denom\n'
1333 # paranoia if CHECK_POLY is set
1334
1335 if CHECK_POLY:
1336 cf = gcd(numer, denom)
1337 if cf != 1:
1338 print('-> Error: {0:d} and {1:d} have a common factor {2:d}.'
1339 ' Please divide it out.'.format(denom, numer, cf))
1340 sys.exit()
1341
1342 if CHECK_POLY and poly_p['m'] and (denom * poly_p['m'] - numer) % fact_p['n']:
1343 print('-> Error: {0:d} * {1:d} + {2:d} != 0 mod {0:d}!'
1344 .format(denom, poly_p['m'], numer, fact_p['n']))
1345 sys.exit(-1)
1346
1347 polyval = 0
1348 if denom and numer:
1349 for i in range(poly_p['degree'] + 1):
1350 polyval += coefvals[i] * numer ** i * denom ** (poly_p['degree'] - i)
1351 elif poly_p['m']:
1352 # print('-> Common root is {0:d}'.format(poly_p['m']))
1353 for i in range(poly_p['degree'], -1, -1):
1354 polyval = poly_p['m'] * polyval + coefvals[i]
1355
1356 if CHECK_POLY and polyval == 0:
1357 print('-> Warning: evaluated polynomial value polyval is negative or zero.')
1358 print('-> This is at least a little strange.')
1359
1360 if CHECK_POLY and polyval % fact_p['n'] != 0:
1361 print('-> Error: evaluated polynomial value polyval is not a multiple of n!')
1362 sys.exit(-1)
1363
1364 # print '-> Evaluated value is polyval\n'
1365
1366 tmp = grep_l('type:', file_lines)
1367 m = re.search('.*type:\s+(gnfs|snfs)', tmp[0])
1368 if not m:
1369 print('-> Error: poly file should contain one of the following lines:')
1370 print('-> type: snfs')
1371 print('-> type: gnfs')
1372 print('-> Please add the appropriate line and re-run.')
1373 sys.exit(-1)
1374 elif m.group(1) == 'gnfs':
1375 fact_p['type'] = 'gnfs'
1376 load_default_parameters('gnfs', fact_p['dec_digits'], poly_p['degree'],
1377 fact_p, pols_p, lats_p, clas_p)
1378 else:
1379 fact_p['type'] = 'snfs'
1380 # We need the difficulty level of the number, which may be
1381 # noticably larger than the number of digits.
1382 fact_p['snfs_difficulty'] = len(str(polyval)) # + math.log(int(str(polyval)[0:1])) - 1
1383 print('-> SNFS_DIFFICULTY is about {0:d}.'.format(fact_p['snfs_difficulty']))
1384 load_default_parameters('snfs', fact_p['snfs_difficulty'], poly_p['degree'],
1385 fact_p, pols_p, lats_p, clas_p)
1386
1387 # Now look for user-supplied parameters.
1388 fact_p['q0'] = 0
1389 with open(NAME + '.poly', 'r') as in_f:
1390 for l in in_f:
1391 l = l.rstrip()
1392 l = chomp_comment(l)
1393 l = l.strip()
1394 if len(l) and l.find(':') != - 1:
1395 token, val = l.split(':')
1396 if len(token) > 0 and len(val) > 0:
1397 if token == 'n':
1398 fact_p['n'] = int(val)
1399 elif token == 'm':
1400 poly_p['m'] = int(val)
1401 elif token == 'rlim':
1402 lats_p['rlim'] = int(val)
1403 elif token == 'alim':
1404 lats_p['alim'] = int(val)
1405 elif token == 'lpbr':
1406 lats_p['lpbr'] = int(val)
1407 elif token == 'lpba':
1408 lats_p['lpba'] = int(val)
1409 elif token == 'mfbr':
1410 lats_p['mfbr'] = int(val)
1411 elif token == 'mfba':
1412 lats_p['mfba'] = int(val)
1413 elif token == 'rlambda':
1414 lats_p['rlambda'] = float(val)
1415 elif token == 'alambda':
1416 lats_p['alambda'] = float(val)
1417 elif token == 'knowndiv':
1418 fact_p['knowndiv'] = int(val)
1419 elif token == 'skew':
1420 poly_p['skew'] = float(val)
1421 elif token == 'q0':
1422 fact_p['q0'] = int(val)
1423 elif token == 'qintsize':
1424 lats_p['qintsize'] = int(val)
1425 elif token == 'lss':
1426 lats_p['lss'] = int(val)
1427 else:
1428 m = re.search('(c|Y).', token)
1429 if m :
1430 poly_p['coeff'][token] = int(val)
1431
1432 fact_p['qstep'] = lats_p['qintsize']
1433 if fact_p['knowndiv']:
1434 if fact_p['n'] % fact_p['knowndiv']:
1435 die('-> Error: knowndiv {0:d} does not divide {0:d}!'
1436 .format(fact_p['knowndiv'], fact_p['n']))
1437 fact_p['ndivfree'] = fact_p['n'] // fact_p['knowndiv']
1438 fact_p['n'] = fact_p['ndivfree']
1439 else:
1440 fact_p['ndivfree'] = fact_p['n']
1441
1442 if probable_prime_p(fact_p['ndivfree'], 10):
1443 die('-> Error: {0:d} is probably prime!'.format(fact_p['n']))
1444 if fact_p['q0'] == 0:
1445 fact_p['q0'] = (lats_p['rlim'] if lats_p['lss'] else lats_p['alim']) // 2
1446 fact_p['qstart'] = fact_p['q0']
1447 check_for_parameter_change(lats_p)
1448
1449# sieving intervals are triples (q_start, q_pos, q_end) that
1450# are kept in a Python deque list in fact_p['q_dq']. In each
1451# sieving step, the q interval (q0 .. q0 + qstep) is divided
1452# into sub-intervals for sieving by the available threads on
1453# this machine. When this step is completed, the intervals
1454# are updated for the next step.
1455
1456# When termainated by a user interrupt the siever will write
1457# a file named 'LASTSPQ<N>', with <N> = 100 * PNUM + TNUM if
1458# multi-threaded or PNUM if single threaded (thread number =
1459# TNUM, processor number = PNUM). This gives the q position
1460# when the siever was stopped. These files from each thread
1461# are used to create a file, JOBNAME.resume, that is used to
1462# restart the sieving. It has the fileds:
1463#
1464# Q0: the Q0 value when terminated
1465# QSTEP: the QSTEP value when terminated
1466# QQ<N>: the sieving interval (ql, qp, qh)
1467# for thread<N> when it was stopped
1468
1469# The following three routines manipulate sieve intervals
1470# and are collected together here to assist in further
1471# improvements. The algorithms are primitive right now.
1472
1473def init_q_dq(fact_p, n_threads, qbase):
1474 inc = int(fact_p['qstep'] // n_threads)
1475 for j in range(n_threads):
1476 if j >= len(fact_p['q_dq']):
1477 fact_p['q_dq'].append(0)
1478 if not fact_p['q_dq'][j]:
1479 fact_p['q_dq'][j] = ((qbase, qbase, qbase + inc))
1480 qbase += inc
1481
1482def update_q_dq(fact_p, n_threads, n_clients):
1483 for j in range(n_threads):
1484 ql, qp, qh = fact_p['q_dq'].popleft()
1485 ql += n_clients * fact_p['qstep']
1486 qh += n_clients * fact_p['qstep']
1487 fact_p['q_dq'].append((ql, ql, qh))
1488
1489def divide_q_dq(fact_p, n_threads):
1490 while len(fact_p['q_dq']) < n_threads:
1491 ql, qp, qh = fact_p['q_dq'].popleft()
1492 t = int((qh - qp) // 2)
1493 fact_p['q_dq'].append((ql, qp, qp + t))
1494 fact_p['q_dq'].append((ql, qp + t, qh))
1495
1496def write_resume_file(n_threads, fact_p):
1497 with open(JOBNAME + '.resume', 'w') as out_f:
1498 if 'q0' in fact_p:
1499 print('Q0: {0:d}'.format(fact_p['q0']), file = out_f)
1500 if 'qstep' in fact_p:
1501 print('QSTEP: {0:d}'.format(fact_p['qstep']), file = out_f)
1502 for j in range(n_threads):
1503 ql, qp, qh = fact_p['q_dq'][j]
1504 print('QQ{0:d}: {1:d}, {2:d}, {3:d}'
1505 .format(j, ql, qp, qh), file = out_f)
1506
1507def read_resume_file():
1508 rn = JOBNAME + '.resume'
1509 fact_p['q_dq'].clear()
1510 with open(rn, 'r') as in_f:
1511 tmp = in_f.readlines()
1512 for l in tmp:
1513 m = re.search('Q0:\s+(\d+)', l)
1514 if m:
1515 q0 = int(m.group(1))
1516 m = re.search('QSTEP:\s+(\d+)', l)
1517 if m:
1518 qstep = int(m.group(1))
1519 m = re.search('QQ(\d+):\s*(\d+),\s*(\d+),\s*(\d+)', l)
1520 if m:
1521 while int(m.group(1)) >= len(fact_p['q_dq']):
1522 fact_p['q_dq'].append(0)
1523 fact_p['q_dq'][int(m.group(1))] = (int(m.group(2)),
1524 int(m.group(3)), int(m.group(4)))
1525 return (q0, qstep)
1526
1527def setup(fact_p, poly_p, lats_p, client, n_clients, n_threads):
1528
1529 # Should we resume from an earlier run? #
1530 resume = os.path.exists(JOBNAME + '.resume')
1531
1532 if resume:
1533 if PROMPTS:
1534 ip = ('-> It appears that an earlier attempt was interrupted. Resume? (y/n) ')
1535 while True:
1536 if sys.version_info[0] < 3:
1537 r = raw_input(ip)
1538 else:
1539 r = input(ip)
1540
1541 if r == 'y' or r == 'Y':
1542 resume = True
1543 break
1544 elif r == 'n' or r == 'N':
1545 resume = False
1546 break
1547
1548 if fact_p['type'] == 'gnfs':
1549 if lats_p['lpbr'] == 25:
1550 lats_p['minrels'] = int(38000.0 * (fact_p['dec_digits'] - 47))
1551 elif lats_p['lpbr'] == 26:
1552 lats_p['minrels'] = int(91000.0 * (fact_p['dec_digits'] - 55))
1553 elif lats_p['lpbr'] == 27:
1554 lats_p['minrels'] = int(150000.0 * (fact_p['dec_digits'] - 61))
1555 elif lats_p['lpbr'] == 28:
1556 lats_p['minrels'] = int(440000.0 * (fact_p['dec_digits'] - 89))
1557 else:
1558 lats_p['minrels'] = int(10.0 ** (fact_p['dec_digits'] / 41.0 + 4.0))
1559 else:
1560 lats_p['minrels'] = int(10.0 ** (fact_p['digs'] / 70.0 + 4.6))
1561# lats_p['minrels'] = int(0.2 * 1.442695 *((2 ** lats_p['lpba']) /lats_p['lpba'] + (2 ** lats_p['lpbr']) / lats_p['lpbr']))
1562 output('-> Estimated minimum relations needed: {0:g}'.format(1.0 * lats_p['minrels']))
1563
1564 # Setup the parameters for sieving ranges. #
1565 if not resume:
1566 if client_id == 1:
1567 if os.path.exists(NAME + '.n') and not os.path.exists(NAME + '.poly'):
1568 delete_file(LOGNAME)
1569 output('-> cleaning up before a restart')
1570
1571 # Clean up any junk leftover from an earlier attempt.
1572 for f in glob.iglob('cols*'):
1573 delete_file(f)
1574 for f in glob.iglob('deps*'):
1575 delete_file(f)
1576 delete_file('factor.easy')
1577 for f in glob.iglob('lpindex*'):
1578 delete_file(f)
1579 for f in glob.iglob(RELSBIN + '*'):
1580 delete_file(f)
1581 delete_file('spairs.out')
1582 delete_file('spairs.out.gz')
1583 delete_file(NAME + '.fb')
1584 for f in glob.iglob('*.afb.0'):
1585 delete_file(f)
1586 for f in glob.iglob('*.afb.1'):
1587 delete_file(f)
1588 delete_file('.params')
1589
1590 # Was a discriminant divisor supplied?
1591 if fact_p['dd']:
1592 with open(LOGNAME, 'a') as out_f:
1593 print('{0:d}'.format(fact_p['dd']), file = out_f)
1594
1595 # Create msieve file
1596 with open(ININAME, 'w') as out_f:
1597 print('{0:d}'.format(fact_p['n']), file = out_f)
1598
1599 with open(DATNAME, 'w') as out_f:
1600 print('N {0:d}'.format(fact_p['n']), file = out_f)
1601
1602 with open(FBNAME, 'w') as out_f:
1603 print('N {0:d}'.format(fact_p['n']), file = out_f)
1604 print('SKEW {0:<10.2f}'.format(poly_p['skew']), file = out_f)
1605
1606 if not 'Y1' in poly_p['coeff']:
1607 poly_p['coeff']['Y1'] = 1
1608 poly_p['coeff']['Y0'] = -poly_p['m']
1609
1610 for key in reversed(sorted(poly_p['coeff'])):
1611 if key[0] == 'c':
1612 print('A{0:d} {1:d}'.format(int(key[1:]), poly_p['coeff'][key]), file = out_f)
1613
1614 print('R1 {0:d}'.format(poly_p['coeff']['Y1']), file = out_f)
1615 print('R0 {0:d}'.format(poly_p['coeff']['Y0']), file = out_f)
1616 print('FAMAX {0:d}'.format(lats_p['alim']), file = out_f)
1617 print('FRMAX {0:d}'.format(lats_p['rlim']), file = out_f)
1618 print('SALPMAX {0:d}'.format(2 ** lats_p['lpba']), file = out_f)
1619 print('SRLPMAX {0:d}'.format(2 ** lats_p['lpbr']), file = out_f)
1620 fact_p['q0'] = fact_p['qstart']
1621 fact_p['q_dq'].clear()
1622
1623 else:
1624
1625 fact_p['q0'] = 0
1626 try:
1627 t = read_resume_file()
1628 fact_p['q0'], fact_p['qstep'] = t
1629 output('-> resuming a block for q from {0:d} to {1:d}'
1630 .format(fact_p['q0'], fact_p['q0'] + fact_p['qstep']))
1631 except IOError:
1632 print('-> Could not determine a starting q value!')
1633 print('-> Please enter a starting point for the special q: ')
1634 if sys.version_info[0] < 3:
1635 tmp = raw_input()
1636 else:
1637 tmp = input()
1638 fact_p['q0'] = int(tmp)
1639 if not fact_p['q0']:
1640 fact_p['q0'] = fact_p['qstart']
1641 fact_p['q_dq'].clear()
1642
1643 # add sieve intervals if the sieve interval queue is not long enough
1644 if 0 < len(fact_p['q_dq']) < n_threads:
1645 divide_q_dq(fact_p, n_threads)
1646 else:
1647 t_q = int(fact_p['q0'] + ((client - 1) % n_clients) * fact_p['qstep'])
1648 init_q_dq(fact_p, n_threads, t_q)
1649
1650def job_out0(sieve_lim, lats_p, out_f):
1651
1652 print('rlim: {0:d}'.format(sieve_lim[1]), file = out_f)
1653 print('alim: {0:d}'.format(sieve_lim[0]), file = out_f)
1654 print('lpbr: {0:d}'.format(lats_p['lpbr']), file = out_f)
1655 print('lpba: {0:d}'.format(lats_p['lpba']), file = out_f)
1656 print('mfbr: {0:d}'.format(lats_p['mfbr']), file = out_f)
1657 print('mfba: {0:d}'.format(lats_p['mfba']), file = out_f)
1658 print('rlambda: {0:g}'.format(lats_p['rlambda']), file = out_f)
1659 print('alambda: {0:g}'.format(lats_p['alambda']), file = out_f)
1660
1661def job_out1(q0, q1, out_f):
1662 print('q0: {0:d}'.format(q0), file = out_f)
1663 print('qintsize: {0:d}'.format(q1 - q0), file = out_f)
1664 print('#q1:{0:d}'.format(q1), file = out_f)
1665
1666def job_out2(clas_p, out_f):
1667 print('a0: {0:s}'.format(clas_p['a0']), file = out_f)
1668 print('a1: {0:s}'.format(clas_p['a1']), file = out_f)
1669 print('b0: {0:s}'.format(clas_p['b0']), file = out_f)
1670 print('b1: {0:s}'.format(clas_p['b1']), file = out_f)
1671
1672def make_sieve_jobfile(FNAME, fact_p, poly_p, lats_p, clas_p = None):
1673
1674 q0 = fact_p['q0']
1675 sieve_type = 1 if q0 > 0 or fact_p['qstep'] > 0 else 0
1676
1677 sieve_lim = [lats_p['alim'], lats_p['rlim'], None]
1678 if sieve_type == 1:
1679 if lats_p['lss'] and lats_p['rlim'] > q0:
1680 sieve_lim[1] = q0 - 1
1681 sieve_lim[2] = '.afb.1'
1682 elif not lats_p['lss'] and lats_p['alim'] > q0:
1683 sieve_lim[0] = q0 - 1
1684 sieve_lim[2] = '.afb.0'
1685
1686 for j in range(SV_THREADS):
1687 ql, qp, qh = fact_p['q_dq'][j]
1688 t_fname = FNAME + '.T' + str(j)
1689 delete_file(t_fname)
1690
1691 output('-> making sieve job for q = {1:d} in {0:d} .. {2:d} as file {3:s}'
1692 .format(ql, qp, qh, t_fname))
1693
1694 with open(t_fname, 'w') as out_f:
1695
1696 print('n: {0:d}'.format(fact_p['n']), file = out_f)
1697 if poly_p['m']:
1698 print('m: {0:d}'.format(poly_p['m']), file = out_f)
1699
1700 # The polynomial coefficients:
1701 for key in reversed(sorted(poly_p['coeff'])) :
1702 print('{0:s} {1:d}'.format(key + ':', poly_p['coeff'][key]), file = out_f)
1703
1704 print('skew: {0:<10.2f}'.format(poly_p['skew']), file = out_f)
1705
1706 job_out0(sieve_lim, lats_p, out_f)
1707 if sieve_type == 1:
1708 job_out1(qp, qh, out_f)
1709 else:
1710 job_out2(clas_p, out_f)
1711
1712 if sieve_lim[2]:
1713 delete_file(t_fname + sieve_lim[2])
1714
1715 return sieve_lim[2]
1716
1717# arg0 = A value, to sieve [-A,A]
1718# arg1 = B0
1719# arg2 = B1, to sieve for b values in [B0, B1].
1720
1721def run_classical_sieve(a, b0, b1):
1722 if not DOCLASSICAL:
1723 return
1724
1725 max_b = 0
1726 lastline = 0
1727
1728 # First, scan the line file and find the largest b-value that was sieved.
1729 line_file = DATNAME + '.line'
1730
1731 if os.path.exists(line_file):
1732 with open(line_file, 'r') as in_f:
1733 tmp = in_f.readline()
1734 tmp = in_f.readline()
1735 bb = int(tmp)
1736 max_b = bb if bb > max_b else max_b
1737
1738 max_b = b0 if max_b < b0 else max_b
1739 if max_b < b1:
1740
1741 with open(FBNAME, 'a') as out_f:
1742 print('SLINE {0:d}'.format(a), file = out_f)
1743
1744 print('-> Line file scanned: resuming classical sieve from b = {0:d}.'
1745 .format(max_b))
1746 ret = run_msieve('-t {0:d} -ns {1:d},{2:d}'.format(LA_THREADS, max_b, b1))
1747 if ret:
1748 die('Interrupted. Terminating...')
1749
1750def log_sieve_time(val):
1751 write_string_to_log('LatSieveTime: {0:g}'.format(val))
1752
1753def read_spq(fact_p):
1754 for j in range(SV_THREADS):
1755 ql, qp, qh = fact_p['q_dq'][j]
1756 try:
1757 with open('.last_spq' + str(100 * PNUM + j), 'r') as in_f:
1758 try:
1759 t = int(in_f.readline().strip())
1760 except ValueError:
1761 pass
1762 else:
1763 if t > qp:
1764 fact_p['q_dq'][j] = (ql, t, qh)
1765 except IOError:
1766 pass
1767
1768def monitor_sieve_threads(procs):
1769 running, ret = False, 0
1770 for p in procs:
1771 retc = p.poll()
1772 if retc == None:
1773 running = True
1774 else:
1775 ret |= retc
1776 return (running, ret)
1777
1778def run_siever(client_id, n_clients, n_threads, fact_p, lats_p):
1779 global procs
1780 old_handler = signal.signal(signal.SIGINT, terminate_threads)
1781 procs = [None] * SV_THREADS
1782 siever_option = '-r' if lats_p['lss'] else '-a'
1783 siever_side = 'rational' if lats_p['lss'] else 'algebraic'
1784 output('-> entering sieving loop')
1785
1786 while not os.path.exists(COLSNAME):
1787
1788 sieve_lim = make_sieve_jobfile(JOBNAME, fact_p, poly_p, lats_p)
1789
1790 if fact_p['q0'] >= 2 ** lats_p['lpba']:
1791 print('-> {0:s} : Severe error!'.format(NAME))
1792 print('-> Current special q = {0:s} has exceeded max. large alg. prime = {1:d)!'
1793 .format(fact_p['q0'], 2 ** lats_p['lpba']))
1794 print('-> You can try increasing LPBA, re-launch this script and cross your fingers.')
1795 print('-> But be aware that if you\'re seeing this, your factorization is taking')
1796 print('-> much longer than it would have with better parameters.')
1797 sys.exit(-1)
1798
1799 write_resume_file(n_threads, fact_p)
1800 output('-> Lattice sieving {0:s} q from {1:d} to {2:d}.'
1801 .format(siever_side, fact_p['q0'], fact_p['q0'] + fact_p['qstep']))
1802 start_time = time.time()
1803
1804 for j in range(n_threads):
1805 sn = SIEVER_OUTPUTNAME + '.T' + str(j)
1806 delete_file(sn)
1807 args = ('-k -o {0:s} -v -n {1:d} {2:s} {3:s}'
1808 .format(sn, 100 * PNUM + j, siever_option, JOBNAME + '.T' + str(j)))
1809 procs[j] = run_exe(lats_p['siever'], args, wait = False)
1810
1811 ret, running, term = 0, True, False
1812 while running and not term:
1813 try:
1814 time.sleep(1)
1815 except IOError:
1816 term = True
1817 break
1818 else:
1819 read_spq(fact_p)
1820 running, ret = monitor_sieve_threads(procs)
1821 if ret or term:
1822 terminate_threads()
1823
1824 for j in range(n_threads):
1825 if sieve_lim:
1826 delete_file(JOBNAME + '.T' + str(j) + sieve_lim)
1827 cat_f(SIEVER_OUTPUTNAME + '.T' + str(j), SIEVER_OUTPUTNAME)
1828 delete_file(SIEVER_OUTPUTNAME + '.T' + str(j))
1829
1830 if ret and ret != -1073741819 or term:
1831 if isinstance(ret, int):
1832 output('-> Return value {0:d}. Updating job file and terminating...'
1833 .format(ret))
1834 cat_f(SIEVER_OUTPUTNAME, SIEVER_ADDNAME)
1835 delete_file(SIEVER_OUTPUTNAME)
1836
1837 write_resume_file(n_threads, fact_p)
1838 # Record the time to the logfile.
1839 log_sieve_time(time.time() - start_time)
1840 die('Terminating...')
1841 else:
1842 fact_p['q0'] += n_clients * fact_p['qstep']
1843 update_q_dq(fact_p, n_threads, n_clients)
1844 write_resume_file(0, fact_p)
1845
1846 if not os.path.exists(SIEVER_OUTPUTNAME):
1847 die('Some error ocurred and no relations were found! Examine the log file.')
1848
1849 if client_id > 1:
1850 cat_f(SIEVER_OUTPUTNAME, SIEVER_ADDNAME)
1851 delete_file(SIEVER_OUTPUTNAME)
1852 else:
1853 # Are there relations coming from somewhere else which should be added in?
1854 if os.path.exists(SIEVER_ADDNAME):
1855 cat_f(SIEVER_ADDNAME, SIEVER_OUTPUTNAME)
1856 delete_file(SIEVER_ADDNAME)
1857
1858 for i in range(n_clients):
1859 san = 'spairs.add.' + str(i + 1)
1860 if os.path.exists(san):
1861 cat_f(san, SIEVER_OUTPUTNAME)
1862 delete_file(san)
1863
1864 cat_f(SIEVER_OUTPUTNAME, DATNAME)
1865 if SAVEPAIRS:
1866 gzip_f(SIEVER_OUTPUTNAME, 'spairs.save.gz')
1867 delete_file(SIEVER_OUTPUTNAME)
1868
1869 if os.path.exists('minrels.txt'):
1870 with open('minrels.txt', 'r') as in_f:
1871 tmp = in_f.readline()
1872 m = re.search('(\d+)', tmp)
1873 if m:
1874 lats_p['minrels'] = int(m.group(1))
1875
1876 lats_p['currels'] = curr_rels = linecount(DATNAME)
1877 pc = (100.0 * curr_rels) / lats_p['minrels']
1878 output('Found {0:d} relations, {1:2.1f}% of the estimated minimum ({2:d}).'
1879 .format(curr_rels, pc, lats_p['minrels']))
1880 if curr_rels > lats_p['minrels']:
1881 ret = run_msieve('-t {0:d} -nc1'.format(LA_THREADS))
1882 if ret:
1883 die('Return value {0:d}. Terminating...'.format(ret))
1884
1885 for j in range(n_threads):
1886 delete_file(JOBNAME + '.T' + str(j))
1887 log_sieve_time(time.time() - start_time)
1888
1889 fact_p['last_spq'] = dict()
1890 delete_file(JOBNAME + '.resume')
1891 signal.signal(signal.SIGINT, old_handler)
1892
1893def summary_name(name, fact_p):
1894 if fact_p['type'] == 'gnfs':
1895 s = 'g' + str(fact_p['dec_digits'])
1896 else :
1897 s = 's' + str(fact_p['snfs_difficulty'])
1898 t = os.path.split(name)
1899 return os.path.join(t[0], s + '-' + t[1] + '.txt')
1900
1901poly_time = 0.0
1902
1903def output_summary(name, fact_p, pols_p, poly_p, lats_p):
1904 global poly_time
1905
1906 # Set the summary file name
1907 sum_name = summary_name(name, fact_p)
1908
1909 # Figure the time scale for this machine.
1910 output('-> Computing time scale for this machine...')
1911 (ret, res) = run_exe(PROCRELS, '-speedtest', out_file = subprocess.PIPE)
1912 if res:
1913 tmp = grep_l('timeunit:', res)
1914 timescale = float(re.sub('timeunit:\s*', '', tmp[0]))
1915 else:
1916 timescale = 0.0
1917
1918 # And gather up some stats.
1919 sieve_time = 0.0
1920 relproc_time = 0.0
1921 matrix_time = 0
1922 sqrt_time = 0
1923 min_q = max_q = 0
1924 prunedmat = rels = ''
1925 rprimes = aprimes = 0
1926 version = 'Msieve-1.40'
1927
1928 with open(LOGNAME, 'r') as in_f:
1929
1930 for l in in_f:
1931 l = l.rstrip()
1932 l = re.sub('\s+$', '', l)
1933 tu = l.split()
1934
1935 m = re.search('for q from (\d+) to (\d+) as file', l)
1936 if m:
1937 t = int(m.group(1))
1938 min_q = t if t < min_q or not min_q else min_q
1939 t = int(m.group(2))
1940 max_q = t if t > max_q or not max_q else max_q
1941
1942 m = re.search('LatSieveTime:\s*(\d+)', l)
1943 if m:
1944 sieve_time += float(m.group(1))
1945
1946 m = re.search('(Msieve.*)$', l)
1947 if m:
1948 version = m.group(1)
1949
1950 m = re.search('RelProcTime: (\S+)', l)
1951 if m:
1952 relproc_time += int(m.group(1))
1953
1954 m = re.search('BLanczosTime: (\S+)', l)
1955 if m:
1956 matrix_time += int(m.group(1))
1957
1958 m = re.search('sqrtTime: (\S+)', l)
1959 if m:
1960 sqrt_time += int(m.group(1))
1961
1962 m = re.search('rational ideals', l)
1963 if m:
1964 rprimes = str(tu[6]) + ' ' + str(tu[7]) + ' ' + str(tu[5])
1965
1966 m = re.search('algebraic ideals', l)
1967 if m:
1968 aprimes = str(tu[6]) + ' ' + str(tu[7]) + ' ' + str(tu[5])
1969
1970 m = re.search('unique relations', l)
1971 if m:
1972 rels = str(tu[-3]) + ' ' + str(tu[-1])
1973
1974 m = re.search('matrix is', l)
1975 if m:
1976 prunedmat = str(tu[7]) + ' x ' + str(tu[9])
1977
1978 m = re.search('p[0-9]+\s+factor', l)
1979 if m:
1980 tmp = int(re.sub('.*factor: ', '', l))
1981 if 1 < len(str(tmp)) < fact_p['dec_digits']:
1982 fact_p['divisors'].append(tmp)
1983
1984 with open('ggnfs.log', 'a') as out_f:
1985 print('Number: {0:s}'.format(NAME), file = out_f)
1986 print('N = {0:d}'.format(fact_p['n']), file = out_f)
1987 for dd in sorted(set(fact_p['divisors'])):
1988 print('factor: {0:d}'.format(dd), file = out_f)
1989
1990 # Convert times from seconds to hours.
1991 poly_time /= 3600.0
1992 sieve_time /= 3600.0
1993 relproc_time /= 3600.0
1994 matrix_time /= 3600.0
1995 sqrt_time /= 3600.0
1996 total_time = poly_time + sieve_time + relproc_time + matrix_time + sqrt_time
1997
1998 with open(sum_name, 'w') as out_f:
1999 print('Number: {0:s}'.format(NAME), file = out_f)
2000 print('N = {0:d} ({1:d} digits)'.format(fact_p['n'], fact_p['dec_digits']), file = out_f)
2001
2002 if fact_p['type'] == 'snfs':
2003 print('SNFS difficulty: {0:d} digits.'.format(fact_p['snfs_difficulty']), file = out_f)
2004
2005 print('Divisors found:', file = out_f)
2006 if fact_p['knowndiv']:
2007 print(' knowndiv: {0:d}'.format(fact_p['knowndiv']), file = out_f)
2008
2009 r = 1
2010 for dd in sorted(set(fact_p['divisors'])):
2011 print('r{0:d}={1:d} (pp{2:d})'.format(r, dd, len(str(dd))), file = out_f)
2012 r += 1
2013
2014 print('Version: {0:s}'.format(version), file = out_f)
2015
2016 if fact_p['type'] == 'mpqs':
2017 print('Completed using mpqs mode', file=out_f)
2018 else:
2019 print('Total time: {0:1.2f} hours.'.format(total_time), file = out_f)
2020 print('Scaled time: {0:1.2f} units (timescale= {1:1.3f}).'
2021 .format(total_time * timescale, timescale))
2022
2023 try:
2024 with open(NAME + '.poly', 'r') as in_f:
2025 print('Factorization parameters were as follows:', file = out_f)
2026 for l in in_f:
2027 l = l.rstrip()
2028 print(l, file = out_f)
2029 except IOError:
2030 pass
2031 siever_side = 'rational' if lats_p['lss'] else 'algebraic'
2032 print('Factor base limits: {0:d}/{1:d}'.format(lats_p['rlim'], lats_p['alim']), file = out_f)
2033 print('Large primes per side: {0:d}'.format(LARGEP), file = out_f)
2034 print('Large prime bits: {0:d}/{1:d}'.format(lats_p['lpbr'], lats_p['lpba']), file = out_f)
2035 print('Sieved {0:s} special-q in [{1:d}, {2:d})'.format(siever_side, min_q, max_q), file = out_f)
2036 print('Total raw relations: {0:d}'.format(lats_p['currels']), file = out_f)
2037 print('Relations: {0:s}'.format(rels), file = out_f)
2038 print('Pruned matrix : {0:s}'.format(prunedmat), file = out_f)
2039 if poly_time:
2040 print('Polynomial selection time: {0:1.2f} hours.'.format(poly_time), file = out_f)
2041 print('Total sieving time: {0:1.2f} hours.'.format(sieve_time), file = out_f)
2042 print('Total relation processing time: {0:1.2f} hours.'.format(relproc_time), file = out_f)
2043 print('Matrix solve time: {0:1.2f} hours.'.format(matrix_time), file = out_f)
2044 print('time per square root: {0:1.2f} hours.'.format(sqrt_time), file = out_f)
2045
2046 if fact_p['type'] == 'snfs':
2047
2048 fact_p['digs'] = fact_p['snfs_difficulty']
2049 DFL = ('{0[type]:s},{0[digs]:d},{1[degree]:d},{3:d},{4:d},{5:g},{6:g},{7:d},'
2050 '{8:d},{9:d},{10:d},{2[rlim]:d},{2[alim]:d},{2[lpbr]:d},{2[lpba]:d},'
2051 '{2[mfbr]:d},{2[mfba]:d},{2[rlambda]:g},{2[alambda]:g},{2[qintsize]:d}'
2052 .format(fact_p, poly_p, lats_p, 0,0,0,0,0,0,0,0))
2053
2054 elif fact_p['type'] == 'gnfs':
2055
2056 fact_p['digs'] = fact_p['dec_digits'] - 1
2057 DFL = ('{0[type]:s},{0[digs]:d},{1[degree]:d},{2[maxs1]:d},{2[maxskew]:d},'
2058 '{2[goodscore]:g},{2[examinefrac]:g},{2[j0]:d},{2[j1]:d},{2[estepsize]:d},'
2059 '{2[maxtime]:d},{3[rlim]:d},{3[alim]:d},{3[lpbr]:d},{3[lpba]:d},'
2060 '{3[mfbr]:d},{3[mfba]:d},{3[rlambda]:g},{3[alambda]:g},{3[qintsize]:d}'
2061 .format(fact_p, poly_p, pols_p, lats_p))
2062
2063 print('Prototype def-par.txt line would be: {0:s}'.format(DFL), file = out_f)
2064 print('total time: {0:1.2f} hours.'.format(total_time), file = out_f)
2065
2066 print(platform.processor(), file = out_f)
2067 print('processors: {0:d}, speed: {1:.2f}GHz'
2068 .format(multiprocessing.cpu_count(), proc_speed()), file = out_f)
2069 t = platform.platform()
2070 if len(t):
2071 print(platform.platform(), file = out_f)
2072 t = platform.python_version_tuple()
2073 print('Running Python {0[0]:s}.{0[1]:s}'.format(t), file = out_f)
2074 output('-> Factorization summary written to {0:s}'.format(sum_name))
2075
2076# ###########################################
2077# ########## Begin execution here ###########
2078# ###########################################
2079
2080print('-> ________________________________________________________________')
2081print('-> | Running factmsieve.py, a Python driver for MSIEVE with GGNFS |')
2082print('-> | sieving support. It is Copyright, 2010, Brian Gladman and is |')
2083print('-> | a conversion of factmsieve.pl that is Copyright, 2004, Chris |')
2084print('-> | Monico. Version {0:s} (Python 2.6 or later) 20th June 2011. |'.format(VERSION))
2085print('-> |______________________________________________________________|')
2086
2087if len(sys.argv) != 2 and len(sys.argv) != 4:
2088 print('USAGE: {0:s} <number file | poly file | msieve poly file> [ id num]'
2089 .format(sys.argv[0]))
2090 print(' where <polynomial file> is a file with the poly info to use')
2091 print(' or <number file> is a file containing the number to factor.')
2092 print(' Optional: id/num specifies that this is client <id> of <num>')
2093 print(' clients total (clients are numbered 1,2,3,...).')
2094 print(' (Multi-client mode is still very experimental - it should only')
2095 print(' be used for testing, and is only intended as a hack for running')
2096 print(' conveniently on a very small number of machines - perhaps 5 - 10).')
2097 sys.exit(-1)
2098
2099GGNFS_PATH = os.path.abspath(GGNFS_PATH)
2100MSIEVE_PATH = os.path.abspath(MSIEVE_PATH)
2101CWD = os.getcwd() + '/'
2102NAME = sys.argv[1]
2103NAME = re.sub('\.(poly|fb|n)', '', NAME)
2104
2105# NOTE: These names are global
2106
2107LOGNAME = NAME + '.log'
2108JOBNAME = NAME + '.job'
2109ININAME = NAME + '.ini'
2110DATNAME = NAME + '.dat'
2111FBNAME = NAME + '.fb'
2112OUTNAME = NAME + '.msp'
2113DEPFNAME = DATNAME + '.dep'
2114COLSNAME = DATNAME + '.cyc'
2115SIEVER_OUTPUTNAME = 'spairs.out'
2116SIEVER_ADDNAME = 'spairs.add'
2117
2118signal.signal(signal.SIGINT, sig_exit)
2119
2120if len(sys.argv) == 4:
2121 client_id = int(sys.argv[2])
2122 num_clients = int(sys.argv[3])
2123 if client_id < 1 or client_id > num_clients:
2124 die('-> Error: client id should be between 1 and the number of clients {0:d}'
2125 .format(num_clients))
2126 PNUM += client_id
2127else:
2128 num_clients = 1
2129 client_id = 1
2130
2131output('-> factmsieve.py (v{0:s})'.format(VERSION), console=False)
2132t = platform.python_version_tuple()
2133print('-> Running Python {0[0]:s}.{0[1]:s}'.format(t))
2134output('-> This is client {0:d} of {1:d}'.format(client_id, num_clients))
2135output('-> Running on {0:d} Core{1:s} with {2:d} hyper-thread{3:s} per Core'
2136 .format(NUM_CORES, ('' if NUM_CORES == 1 else 's'),
2137 THREADS_PER_CORE, ('' if THREADS_PER_CORE == 1 else 's')))
2138output('-> Working with NAME = {0:s}'.format(NAME))
2139
2140if client_id > 1:
2141 JOBNAME += '.' + str(client_id)
2142 SIEVER_OUTPUTNAME += '.' + str(client_id)
2143 SIEVER_ADDNAME += '.' + str(client_id)
2144
2145check_binary(MSIEVE)
2146
2147# Is there a poly file already, or do we need to create one?
2148if not os.path.exists(NAME + '.poly'):
2149 if os.path.exists(NAME + '.fb'):
2150 output('-> Converting msieve polynomial (*.fb) to ggnfs (*.poly) format.')
2151 fb_to_poly()
2152 elif not os.path.exists(NAME + '.n'):
2153 die('-> Could not find the file {0:s}.n'.format(NAME))
2154 else:
2155 delete_file(PARAMFILE)
2156 with open(NAME + '.n', 'r') as in_f:
2157 numberinf = in_f.readlines()
2158 t = grep_l('^n:', numberinf)
2159 if len(t):
2160 m = re.search('n:\s*(\d+)', t[0])
2161 if len(t) and m and len(m.group(1)):
2162 fact_p['n'] = int(m.group(1))
2163 fact_p['dec_digits'] = len(m.group(1))
2164 print('-> Found n = {0:d}.'.format(fact_p['n']))
2165 else:
2166 output('-> Could not find a number in the file {0:s}.n'.format(NAME))
2167 die('-> Did you forget the \'n:\' tag?')
2168 if probable_prime_p(fact_p['n'], 10):
2169 die ('-> Error: n is probably prime!')
2170
2171 if fact_p['n'] >= 2 ** MIN_NFS_BITS:
2172 poly_start = time.time()
2173 print('-> Polynomial file {0:s}.poly does not exist!'.format(NAME))
2174 if num_clients > 1:
2175 die('-> Script does not support polynomial search across multiple clients!')
2176 output('-> Running polynomial selection ...')
2177 if fact_p['dec_digits'] < 98:
2178 USE_KLEINJUNG_FRANKE_PS = 0
2179 if USE_KLEINJUNG_FRANKE_PS:
2180 run_pol5(fact_p, pol5_p, lats_p, clas_p)
2181 elif USE_MSIEVE_POLY:
2182 if not run_msieve_poly(fact_p, POLY_MIN_LC, POLY_MAX_LC):
2183 sys.exit(0)
2184 else:
2185 run_poly_select(fact_p, pols_p, lats_p, clas_p)
2186 if not os.path.exists(NAME + '.poly'):
2187 die('Polynomial selection failed.')
2188 poly_time = time.time() - poly_start
2189 else:
2190 output('-> Running mpqs')
2191 fact_p['type'] = 'mpqs'
2192 if msieve_mpqs(fact_p):
2193 output_summary(NAME, fact_p, pols_p, poly_p, lats_p)
2194 sys.exit(0)
2195
2196# Read and verify parameters for the factorization.
2197if not os.path.exists(NAME + '.poly'):
2198 die('Cannot find {0:s}.poly'.format(NAME))
2199
2200read_parameters(fact_p, poly_p, lats_p)
2201check_parameters(fact_p, poly_p, lats_p)
2202if not os.path.exists(DEPFNAME):
2203 output('-> Running setup ...')
2204 setup(fact_p, poly_p, lats_p, client_id, num_clients, SV_THREADS)
2205
2206# Do some classical sieving, if needed/applicable.
2207if client_id > 1:
2208 DOCLASSICAL = 0
2209if DOCLASSICAL:
2210 output('-> Running classical siever ...')
2211 run_classical_sieve(clas_p['cl_a'], 1, clas_p['cl_b'])
2212
2213# Finally, sieve until we have reached the desired min # of FF's.
2214output('-> Running lattice siever ...')
2215run_siever(client_id, num_clients, SV_THREADS, fact_p, lats_p)
2216
2217if client_id > 1:
2218 print('Client {0:d} terminating...'.format(client_id))
2219 sys.exit(0)
2220
2221# Obviously, the matrix solving step.
2222if not os.path.exists(DEPFNAME):
2223 output('-> Running matrix solving step ...')
2224 ret = run_msieve('-t {0:d} -nc2'.format(LA_THREADS))
2225 if ret:
2226 die('Return value {0:d}. Terminating...'.format(ret))
2227 if not os.path.exists(DEPFNAME):
2228 die('Some error occurred and matsolve did not record dependencies.')
2229else:
2230 print('-> File \'deps\' already exists. Proceeding to sqrt step.')
2231
2232# Do as many square root jobs as needed to get the final factorization
2233if not os.path.exists(summary_name(NAME, fact_p)):
2234 output('-> Running square root step ...')
2235 ret = run_msieve('-t {0:d} -nc3'.format(LA_THREADS))
2236
2237if CLEANUP:
2238 for f in glob.iglob('cols*'):
2239 delete_file(f)
2240 for f in glob.iglob('deps*'):
2241 delete_file(f)
2242 delete_file('factor.easy')
2243 for f in glob.iglob('rels*'):
2244 delete_file(f)
2245 delete_file('spairs.out')
2246 delete_file('spairs.out.gz')
2247 delete_file(NAME + '.fb')
2248 for f in glob.iglob('*.afb.0'):
2249 delete_file(f)
2250 delete_file('tmpdata.000')
2251 delete_file(PARAMFILE)
2252
2253output_summary(NAME, fact_p, pols_p, poly_p, lats_p)
2254for i in range(5):
2255 sys.stdout.write('\a')
2256 sys.stdout.flush()
2257 time.sleep(0.5)