## "Fossies" - the Fresh Open Source Software Archive

### Member "NZMATH-1.2.0/nzmath/factor/mpqs.py" (19 Nov 2012, 26965 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "mpqs.py" see the Fossies "Dox" file reference documentation.

```    1 import math
2 import time
3 import logging
4 import nzmath.arith1 as arith1
5 import nzmath.gcd as gcd
6 import nzmath.prime as prime
7 import nzmath.factor.find as find
8
9 # TODO: add psyco to config entry
10 ##         try:
11 ##             import psyco
12 ##             psyco.full()
13 ##         except ImportError:
14 ##             pass
15
16 _log = logging.getLogger('nzmath.factor.mpqs')
17
18 class QS(object):
19     def __init__(self, n, sieverange, factorbase):
20         self.number = n
21         self.sqrt_n = int(math.sqrt(n))
22         for i in [2,3,5,7,11,17,19]:
23             if n % i == 0:
24                 raise ValueError("This number is divided by %d" % i)
25
26         self.digit = arith1.log(self.number, 10) + 1
27         self.Srange = sieverange
28         self.FBN = factorbase
29
30         self.move_range = range(self.sqrt_n-self.Srange, self.sqrt_n+self.Srange+1)
31
32         i = 0
33         k = 0
34         factor_base = [-1]
35         FB_log = [0]
36         while True:
37             ii = primes_table[i]
38             if arith1.legendre(self.number, ii) == 1:
39                 factor_base.append(ii)
40                 FB_log.append(primes_log_table[i])
41                 k += 1
42                 i += 1
43                 if k == self.FBN:
44                     break
45             else:
46                 i += 1
47
48         self.FB = factor_base
49         self.FB_log = FB_log
50         self.maxFB = factor_base[-1]
51         N_sqrt_list = []
52         for i in self.FB:
53             if i != 2 and i != -1:
54                 e = int(math.log(2*self.Srange, i))
55                 N_sqrt_modp = sqroot_power(self.number, i, e)
56                 N_sqrt_list.append(N_sqrt_modp)
57         self.solution = N_sqrt_list  #This is square roots of N on Z/pZ, p in factor base.
58
59         poly_table = []
60         log_poly = []
61         minus_val = []
62         for j in self.move_range:
63             jj = (j**2)-self.number
64             if jj < 0:
65                 jj = -jj
66                 minus_val.append(j-self.sqrt_n+self.Srange)
67             elif jj == 0:
68                 jj = 1
69             lj = int((math.log(jj)*30)*0.97)  # 0.97 is an erroe
70             poly_table.append(jj)
71             log_poly.append(lj)
72         self.poly_table = poly_table  # This is Q(x) value , x in [-M+sqrt_n,M+sqrt_n].
73         self.log_poly = log_poly      # This is log(Q(x)) value.
74         self.minus_check = minus_val # This is "x" that Q(x) is minus value.
75
76     def run_sieve(self):
77         T = time.time()
78         M = self.Srange
79         start_location = []
80         logp = [0]*(2*M+1)
81         j = 2
82         for i in self.solution:
83             k = 0
84             start_p = []
85             while k < len(i):
86                 q = int((self.sqrt_n)/(self.FB[j]**(k+1)))
87                 s_1 = q*(self.FB[j]**(k+1))+i[k][0]
88                 s_2 = q*(self.FB[j]**(k+1))+i[k][1]
89                 while True:
90                     if s_1 < self.sqrt_n-M:
91                         s_1 = s_1+(self.FB[j]**(k+1))
92                         break
93                     else:
94                         s_1 = s_1-(self.FB[j]**(k+1))
95                 while True:
96                     if s_2 < self.sqrt_n-M:
97                         s_2 = s_2+(self.FB[j]**(k+1))
98                         break
99                     else:
100                         s_2 = s_2-(self.FB[j]**(k+1))
101                 start_p.append([s_1-self.sqrt_n+M,s_2-self.sqrt_n+M])
102
103                 k += 1
104             start_location.append(start_p)
105             j += 1
106         self.start_location = start_location
107
108         if self.poly_table[0] & 1 == 0:
109             i = 0
110             while i <= 2*M:
111                 j = 1
112                 while True:
113                     if self.poly_table[i] % (2**(j+1)) == 0:
114                         j+=1
115                     else:
116                         break
117                 logp[i] += self.FB_log[1]*j
118                 i += 2
119         else:
120             i = 1
121             while i <= 2*M:
122                 j = 1
123                 while True:
124                     if self.poly_table[i] % (2**(j+1)) == 0:
125                         j += 1
126                     else:
127                         break
128                 logp[i] += self.FB_log[1]*j
129                 i += 2
130         L = 2
131         for j in self.start_location:
132             k = 0
133             while k < len(j):
134                 s_1 = j[k][0]
135                 s_2 = j[k][1]
136                 h_1 = 0
137                 h_2 = 0
138                 while s_1+h_1 <= 2*M:
139                     logp[s_1+h_1] += self.FB_log[L]
140                     h_1 += self.FB[L]**(k+1)
141                 while s_2+h_2 <= 2*M:
142                     logp[s_2+h_2] += self.FB_log[L]
143                     h_2 += self.FB[L]**(k+1)
144                 k += 1
145             L += 1
146
147         self.logp = logp
148         smooth = []
149         for t in range(2*M+1):
150             if logp[t] >= self.log_poly[t]:
151                 poly_val = self.poly_table[t]
152                 index_vector = []
153                 for p in self.FB:
154                     if p == -1:
155                         if t in self.minus_check:
156                             index_vector.append(1)
157                         else:
158                             index_vector.append(0)
159                     else:
160                         r = 0
161                         while poly_val % (p**(r+1)) == 0:
162                             r += 1
163                         v = r & 1
164                         index_vector.append(v)
165                 smooth.append([index_vector, (poly_val, t+self.sqrt_n-M)])
166         _log.info(" Sieve time is %f sec" % (time.time()-T))
167         _log.info(" Found smooth numbers are %d / %d" % (len(smooth), len(self.FB)))
168         self.smooth = smooth
169         return smooth
170
171
172 class MPQS(object):
173     def __init__(self, n, sieverange=0, factorbase=0, multiplier=0):
174         self.number = n
175         _log.info("The number is %d MPQS starting" % n)
176
177         if prime.primeq(self.number):
178             raise ValueError("This number is Prime Number")
179         for i in [2,3,5,7,11,13]:
180             if n % i == 0:
181                 raise ValueError("This number is divided by %d" % i)
182
183         self.sievingtime = 0
184         self.coefficienttime = 0
185         self.d_list = []
186         self.a_list = []
187         self.b_list = []
188
189         #Decide prameters for each digits
190         self.digit = arith1.log(self.number, 10) + 1
191
192         if sieverange != 0:
193             self.Srange = sieverange
194             if factorbase != 0:
195                 self.FBN = factorbase
196             elif self.digit < 9:
197                 self.FBN = parameters_for_mpqs[0][1]
198             else:
199                 self.FBN = parameters_for_mpqs[self.digit-9][1]
200         elif factorbase != 0:
201             self.FBN = factorbase
202             if self.digit < 9:
203                 self.Srange = parameters_for_mpqs[0][0]
204             else:
205                 self.Srange = parameters_for_mpqs[self.digit-9][0]
206         elif self.digit < 9:
207             self.Srange = parameters_for_mpqs[0][0]
208             self.FBN = parameters_for_mpqs[0][1]
209         elif self.digit > 53:
210             self.Srange = parameters_for_mpqs[44][0]
211             self.FBN = parameters_for_mpqs[44][1]
212         else:
213             self.Srange = parameters_for_mpqs[self.digit-9][0]
214             self.FBN = parameters_for_mpqs[self.digit-9][1]
215
216         self.move_range = range(-self.Srange, self.Srange+1)
217
218         # Decide k such that k*n = 1 (mod4) and k*n has many factor base
219         if multiplier == 0:
220             self.sqrt_state = []
221             for i in [3,5,7,11,13]:
222                 s = arith1.legendre(self.number, i)
223                 self.sqrt_state.append(s)
224
225             if self.number % 8 == 1 and self.sqrt_state == [1,1,1,1,1]:
226                 k=1
227             else:
228                 index8 = (self.number & 7) >> 1
229                 j = 0
230                 while self.sqrt_state != prime_8[index8][j][1]:
231                     j += 1
232                 k = prime_8[index8][j][0]
233         else:
234             if n & 3 == 1:
235                 k = 1
236             else:
237                 if multiplier == 1:
238                     raise ValueError("This number is 3 mod 4 ")
239                 else:
240                     k = multiplier
241         self.number = k*self.number
242         self.multiplier = k
243
244         _log.info("%d - digits Number" % self.digit)
245         _log.info("Multiplier is %d" % self.multiplier)
246
247         # Table of (log p) , p in FB
248         i = 0
249         k = 0
250         factor_base = [-1]
251         FB_log = [0]
252         while k < self.FBN:
253             ii = primes_table[i]
254             if arith1.legendre(self.number,ii) == 1:
255                 factor_base.append(ii)
256                 FB_log.append(primes_log_table[i])
257                 k += 1
258             i += 1
259
260         self.FB = factor_base
261         self.FB_log = FB_log
262         self.maxFB = factor_base[-1]
263
264         # Solve x^2 = n (mod p^e)
265         N_sqrt_list = []
266         for i in self.FB:
267             if i != 2 and i != -1:
268                 e = int(math.log(2*self.Srange, i))
269                 N_sqrt_modp = sqroot_power(self.number, i, e)
270                 N_sqrt_list.append(N_sqrt_modp)
271         self.Nsqrt = N_sqrt_list
272
273     def make_poly(self):
274         """
275         Make coefficients of f(x)= ax^2+b*x+c
276         """
277         T = time.time()
278         if self.d_list == []:
279             d = int(math.sqrt((math.sqrt(self.number)/(math.sqrt(2)*self.Srange))))
280             if d& 1 == 0:
281                 if (d+1)& 3 == 1: #case d=0 mod4
282                     d += 3
283                 else:
284                     d += 1       #case d=2 mod4
285             elif d& 3 == 1:       #case d=1 mod4
286                 d += 2
287                                  #case d=3 mod4
288         else:
289             d = self.d_list[-1]
290
291         while d in self.d_list or not prime.primeq(d) or \
292               arith1.legendre(self.number,d) != 1 or d in self.FB:
293             d += 4
294         a = d**2
295         h_0 = pow(self.number, (d-3) >> 2, d)
296         h_1 = (h_0*self.number) % d
297         h_2 = ((arith1.inverse(2,d)*h_0*(self.number - h_1**2))//d) % d
298         b = (h_1 + h_2*d) % a
299         if b& 1 == 0:
300             b = b - a
301
302         self.d_list.append(d)
303         self.a_list.append(a)
304         self.b_list.append(b)
305
306         # Get solution of  F(x) = 0 (mod p^i)
307         solution = []
308         i = 0
309         for s in self.Nsqrt:
310             k = 0
311             p_solution = []
312             ppow = 1
313             while k < len(s):
314                 ppow *= self.FB[i+2]
315                 a_inverse = arith1.inverse(2*self.a_list[-1], ppow)
316                 x_1 = ((-b + s[k][0])*a_inverse) % ppow
317                 x_2 = ((-b + s[k][1])*a_inverse) % ppow
318                 p_solution.append([x_1, x_2])
319                 k += 1
320             i += 1
321             solution.append(p_solution)
322         self.solution = solution
323         self.coefficienttime += time.time() - T
324
325     def run_sieve(self):
326         self.make_poly()
327         T = time.time()
328         M = self.Srange
329         a = self.a_list[-1]            #
330         b = self.b_list[-1]            # These are coefficients of F(x)
331         c = (b**2-self.number)//(4*a)   #
332         d = self.d_list[-1]            #
333
334         self.poly_table = []  # This is F(x) value , x in [-M,M].
335         self.log_poly = []    # This is log(F(x)) value.
336         self.minus_check = [] # This is "x" that F(x) is minus value.
337         for j in self.move_range:
338             jj = (a * j + b) * j + c
339             if jj < 0:
340                 jj = -jj
341                 self.minus_check.append(j + M)
342             elif jj == 0:
343                 jj = 1
344             lj = int((math.log(jj)*30)*0.95)  # 0.95 is an erroe
345             self.poly_table.append(jj)
346             self.log_poly.append(lj)
347
348         y = arith1.inverse(2*d, self.number)
349         start_location = []
350         logp = [0]*(2*M+1)
351         j = 2
352         for i in self.solution:
353             start_p = []
354             ppow = 1
355             for k in range(len(i)):
356                 ppow *= self.FB[j]
357                 q = -M // ppow
358                 s_1 = (q + 1) * ppow + i[k][0]
359                 s_2 = (q + 1) * ppow + i[k][1]
360                 while s_1 + M >= ppow:
361                     s_1 = s_1 - ppow
362                 while s_2 + M >= ppow:
363                     s_2 = s_2 - ppow
364                 start_p.append([s_1+M, s_2+M])
365             start_location.append(start_p)
366             j += 1
367         self.start_location = start_location
368
369         i = self.poly_table[0] & 1
370         while i <= 2*M:
371             j = 1
372             while self.poly_table[i] % (2**(j+1)) == 0:
373                 j += 1
374             logp[i] += self.FB_log[1]*j
375             i += 2
376         L = 2
377         for plocation in self.start_location:
378             for k in range(len(plocation)):
379                 s_1 = plocation[k][0]
380                 s_2 = plocation[k][1]
381                 ppow = self.FB[L]**(k+1)
382                 while s_1 <= 2*M:
383                     logp[s_1] += self.FB_log[L]
384                     s_1 += ppow
385                 while s_2 <= 2*M:
386                     logp[s_2] += self.FB_log[L]
387                     s_2 += ppow
388             L += 1
389
390         self.logp = logp
391         smooth = []
392         for t in range(2*M+1):
393             if logp[t] >= self.log_poly[t]:
394                 poly_val = self.poly_table[t]
395                 index_vector = []
396                 H = (y*(2*a*(t-self.Srange)+b))%self.number
397                 for p in self.FB:
398                     if p == -1:
399                         if t in self.minus_check:
400                             index_vector.append(1)
401                         else:
402                             index_vector.append(0)
403                     else:
404                         r = 0
405                         while poly_val % (p**(r+1)) == 0:
406                             r += 1
407                         v = r & 1
408                         index_vector.append(v)
409                 smooth.append([index_vector, (poly_val, H)])
410         self.sievingtime += time.time() - T
411         return smooth
412
413     def get_vector(self):
414         T = time.time()
415         P = len(self.FB)
416         if P < 100:
417             V = -5
418         else:
419             V = 0
420         smooth = []
421         i = 0
422         while P*1 > V:
423             n = self.run_sieve()
424             V += len(n)
425             smooth += n
426             i += 1
427             if i % 50 == 0:
428                 if i == 50:
429                     _log.info("*/ Per 50 times changing poly report /* ")
430                 _log.info("Time of deciding coefficient = %f sec" % self.coefficienttime)
431                 _log.info("Sieving Time = %f sec" % self.sievingtime)
432                 _log.info("Find smooth numbers are %d / %d" % (V, P))
433         if P < 100:
434             V += 5
435         self.smooth = smooth
436         _log.info("*/ Total %d times changing poly report /*" % i)
437         _log.info("Time of deciding coefficient = %f sec" % self.coefficienttime)
438         _log.info("Sieving Time = %f sec" % self.sievingtime)
439         _log.info("Found smooth numbers are %d / %d" % (V, P))
440         _log.info("Total time of getting enough smooth numbers = %f sec" % (time.time()-T))
441         return smooth
442
443
444 class Elimination(object):
445     def __init__(self, smooth):
446         self.vector = []
447         self.history = []
448         i = 0
449         for vec in smooth:
450             self.vector.append(vec[0])
451             self.history.append({i:1})
452             i += 1
453         self.FB_number = len(self.vector[0])
454         self.row_size = len(self.vector)
455         self.historytime = 0
456
458         V_i = self.vector[i]
459         V_j = self.vector[j]
460         k = 0
461         while k < len(V_i):
462             if V_i[k] == 1:
463                 if V_j[k] == 1:
464                     V_j[k] = 0
465                 else:
466                     V_j[k] = 1
467             k += 1
468
469     def transpose(self):
470         Transe_vector = []
471         i = 0
472         while i < self.FB_number:
473             j = 0
474             vector = []
475             while j < self.row_size:
476                 vector.append(self.vector[j][i])
477                 j += 1
478             Transe_vector.append(vector)
479             i += 1
480         self.Transe = Transe_vector
481
483         T = time.time()
484         H_i = self.history[i].keys()
485         H_j = self.history[j].keys()
486         for k in H_i:
487             if k in H_j:
488                 del self.history[j][k]
489             else:
490                 self.history[j][k] = 1
491         self.historytime += time.time() - T
492
493     def gaussian(self):
494         T = time.time()
495         pivot = []
496         FBnum = self.FB_number
497         Smooth = len(self.vector)
498         for j in range(self.FB_number):
499             for k in range(Smooth):
500                 if k in pivot or not self.vector[k][j]:
501                     continue
502                 pivot.append(k)
503                 V_k = self.vector[k]
504                 for h in range(Smooth):
505                     if h in pivot or not self.vector[h][j]:
506                         continue
508                     V_h = self.vector[h]
509                     for q in range(j, FBnum):
510                         if V_k[q]:
511                             V_h[q] = not V_h[q]
512                 break
513         self.pivot = pivot
514
515         zero_vector = []
516         for check in range(Smooth):
517             if check not in pivot:
518                 g = 0
519                 while g < FBnum:
520                     if self.vector[check][g] == 1:
521                         break
522                     g += 1
523                 if g == FBnum:
524                     zero_vector.append(check)
525         _log.info("Time of Gaussian Elimination = %f sec" % (time.time() - T))
526         return zero_vector
527
528 def qs(n, s, f):
529     """
530     This is main function of QS
531     Arguments are (composite_number, sieve_range, factorbase_size)
532     You must input these 3 arguments.
533     """
534     a = time.time()
535     Q = QS(n, s, f)
536     _log.info("Sieve range is [ %d , %d ] , Factorbase size = %d , Max Factorbase %d" % (Q.move_range[0], Q.move_range[-1], len(Q.FB), Q.maxFB))
537     Q.run_sieve()
538     V = Elimination(Q.smooth)
539     A = V.gaussian()
540     _log.info("Found %d linearly dependent relations" % len(A))
542     N_factors = []
543     for i in A:
544         B = V.history[i].keys()
545         X = 1
546         Y = 1
547         for j in B:
548             X *= Q.smooth[j][1][0]
549             Y *= Q.smooth[j][1][1]
550             Y = Y % Q.number
551         X = sqrt_modn(X, Q.number)
554         if k != 0:
555             factor = gcd.gcd(k, Q.number)
556             if factor not in N_factors and factor != 1 and \
557                factor != Q.number and prime.primeq(factor) == 1:
558                 N_factors.append(factor)
559     N_factors.sort()
560     _log.info("Total time = %f sec" % (time.time()-a))
561     _log.info(str(N_factors))
562
563 def mpqs(n, s=0, f=0, m=0):
564     """
565     This is main function of MPQS.
566     Arguments are (composite_number, sieve_range, factorbase_size, multiplier)
567     You must input composite_number at least.
568     """
569     T = time.time()
570     M = MPQS(n, s, f, m)
571     _log.info("Sieve range is [ %d , %d ] , Factorbase size = %d , Max Factorbase %d" % (M.move_range[0], M.move_range[-1], len(M.FB), M.maxFB))
572     M.get_vector()
573     N = M.number // M.multiplier
574     V = Elimination(M.smooth)
575     A = V.gaussian()
576     _log.info("Found %d linerly dependent relations" % len(A))
578     N_prime_factors = []
579     N_factors = []
580     output = []
581     for i in A:
582         B = V.history[i].keys()
583         X = 1
584         Y = 1
585         for j in B:
586             X *= M.smooth[j][1][0]
587             Y *= M.smooth[j][1][1]
588             Y = Y % M.number
589         X = sqrt_modn(X, M.number)
590         if X != Y:
592     NN = 1
594         factor = gcd.gcd(k, N)
595         if factor not in N_factors and factor != 1 and factor != N \
596                and factor not in N_prime_factors:
597             if prime.primeq(factor):
598                 NN = NN*factor
599                 N_prime_factors.append(factor)
600             else:
601                 N_factors.append(factor)
602
603     _log.info("Total time = %f sec" % (time.time() - T))
604
605     if NN == N:
606         _log.debug("Factored completely!")
607         N_prime_factors.sort()
608         for p in N_prime_factors:
609             N = N // p
610             i = arith1.vp(N, p, 1)[0]
611             output.append((p, i))
612         return output
613     elif NN != 1:
614         f = N // NN
615         if prime.primeq(f):
616             N_prime_factors.append(f)
617             _log.debug("Factored completely !")
618             N_prime_factors.sort()
619             for p in N_prime_factors:
620                 N = N // p
621                 i = arith1.vp(N, p, 1)[0]
622                 output.append((p, i))
623             return output
624
625     for F in N_factors:
626         for FF in N_factors:
627             if F != FF:
628                 Q = gcd.gcd(F, FF)
629                 if prime.primeq(Q) and Q not in N_prime_factors:
630                     N_prime_factors.append(Q)
631                     NN = NN*Q
632
633     N_prime_factors.sort()
634     for P in N_prime_factors:
635         i, N = arith1.vp(N, P)
636         output.append((P, i))
637
638     if  N == 1:
639         _log.debug("Factored completely!! ")
640         return output
641
642     for F in N_factors:
643         g = gcd.gcd(N, F)
644         if prime.primeq(g):
645             N_prime_factors.append(g)
646             N = N // g
647             i = arith1.vp(N, g, 1)[0]
648             output.append((g, i))
649     if N == 1:
650         _log.debug("Factored completely !! ")
651         return output
652     elif prime.primeq(N):
653         output.append((N, 1))
654         _log.debug("Factored completely!!! ")
655         return output
656     else:
657         N_factors.sort()
658         _log.error("Sorry, not factored completely")
659         return output, N_factors
660
661
662
663 ########################################################
664 #                                                      #
665 # Following functions are subfunction for main program #
666 #                                                      #
667 ########################################################
668
669 def eratosthenes(n):
670     if n < 2:
671         raise ValueError("Input value must be bigger than 1")
672     else:
673         primes = [2]
674         sieves = [1] * (((n+1) >> 1) - 1)
675         k = 3
676         i = 0
677         sieves_len = len(sieves)
678         while k <= math.sqrt(n):
679             if sieves[i] == 1:
680                 j = 1
681                 while i+j*k <= sieves_len-1:
682                     sieves[i+j*k] = 0
683                     j = j+1
684                 k, i = k+2, i+1
685             else:
686                 k, i = k+2, i+1
687         y = 3
688         for x in sieves:
689             if x == 1:
690                 primes.append(y)
691             y = y+2
692         return primes
693
694 def prime_mod8(n):
695     """
696     This is table for choosing multiplier which makes N to have
697     factorbase(2,3,5,7,11,13)
698     """
699     primes = eratosthenes(n)
700     PrimeList = {1:[], 3:[], 5:[], 7:[]}
701     LegendreList = {1:[], 3:[], 5:[], 7:[]}
702     sp = [2, 3, 5, 7, 11, 13]
703     for p in primes:
704         if p not in sp:
705             leg = [arith1.legendre(p, q) for q in sp[1:]]
706             if leg not in PrimeList[p & 7]:
707                 LegendreList[p & 7].append(leg)
708                 PrimeList[p & 7].append([p, leg])
709     return [PrimeList[1], PrimeList[3], PrimeList[5], PrimeList[7]]
710
711 def eratosthenes_log(n):
712     primes = eratosthenes(n)
713     primes_log = []
714     for i in primes:
715         l = int(math.log(i)*30) #30 is scale
716         primes_log.append(l)
717     return primes_log
718
719 def sqrt_modn(n, modulo):
720     import nzmath.factor.methods as methods
721     factorOfN = methods.trialDivision(n)
722     prod = 1
723     for p, e in factorOfN:
724         prod = (prod * pow(p, e >> 1, modulo)) % modulo
725     return prod
726
727 def sqroot_power(a, p, n):
728     """
729     return squareroot of a mod p^k for k = 2,3,...,n
730     """
731     x = arith1.modsqrt(a, p)
732     answer = [[x, p - x]]
733     ppower = p
734     inverse = arith1.inverse(x << 1, p)
735     for i in range(n - 1):
736         x += (a - x ** 2) // ppower * inverse % p * ppower
737         ppower *= p
740
741 #################
742 # Initial items #
743 #################
744 primes_table = eratosthenes(10**5)
745 primes_log_table = eratosthenes_log(10**5)
746 prime_8 = prime_mod8(8090)
747 parameters_for_mpqs = [[100,20], #9 -digits
748                       [100,21], #10
749                       [100,22], #11
750                       [100,24], #12
751                       [100,26], #13
752                       [100,29], #14
753                       [100,32], #15
754                       [200,35], #16
755                       [300,40], #17
756                       [300,60], #18
757                       [300,80], #19
758                       [300,100], #20
759                       [300,100], #21
760                       [300,120], #22
761                       [300,140], #23
762                       [600,160], #24
763                       [900,180], #25
764                       [1200,200], #26
765                       [1000,220], #27
766                       [2000,240], #28
767                       [2000,260], #29
768                       [2000,325], #30
769                       [2000,355], #31
770                       [2000,375], #32
771                       [3000,400], #33
772                       [2000,425], #34
773                       [2000,550], #35
774                       [3000,650], #36
775                       [5000,750], #37
776                       [4000,850], #38
777                       [4000,950], #39
778                       [5000,1000], #40
779                       [14000,1150], #41
780                       [15000,1300], #42
781                       [15000,1600], #43
782                       [15000,1900], #44
783                       [15000,2200], #45
784                       [20000,2500], #46
785                       [25000,2500], #47
786                       [27500,2700], #48
787                       [30000,2800], #49
788                       [35000,2900], #50
789                       [40000,3000], #51
790                       [50000,3200], #52
791                       [50000,3500]] #53
792
793
794 ###
795 ### only find a factor
796 ###
797
798 def mpqsfind(n, s=0, f=0, m=0, verbose=False):
799     """
800     This is main function of MPQS.
801     Arguments are (composite_number, sieve_range, factorbase_size, multiplier)
802     You must input composite_number at least.
803     """
804     # verbosity
805     if verbose:
806         _log.setLevel(logging.DEBUG)
807         _log.debug("verbose")
808     else:
809         _log.setLevel(logging.NOTSET)
810
811     starttime = time.time()
812     M = MPQS(n, s, f, m)
813     _log.info("Sieve range is [%d, %d]" % (M.move_range[0], M.move_range[-1]))
814     _log.info("Factorbase size = %d, Max Factorbase %d" % (len(M.FB), M.maxFB))
815     M.get_vector()
816     N = M.number // M.multiplier
817     V = Elimination(M.smooth)
818     A = V.gaussian()
819     _log.info("Found %d linearly dependent relations" % len(A))
820     differences = []
821     for i in A:
822         B = V.history[i].keys()
823         X = 1
824         Y = 1
825         for j in B:
826             X *= M.smooth[j][1][0]
827             Y *= M.smooth[j][1][1]
828             Y = Y % M.number
829         X = arith1.floorsqrt(X) % M.number
830         if X != Y:
831             differences.append(X-Y)
832
833     for diff in differences:
834         divisor = gcd.gcd(diff, N)
835         if 1 < divisor < N:
836             _log.info("Total time = %f sec" % (time.time() - starttime))
837             return divisor
```