"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 
  457     def vector_add(self, i, j):
  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 
  482     def history_add(self, i, j):
  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
  507                     self.history_add(k, h)
  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))
  541     answerX_Y = []
  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)
  552         answerX_Y.append(X-Y)
  553     for k in answerX_Y:
  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))
  577     answerX_Y = []
  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:
  591             answerX_Y.append(X-Y)
  592     NN = 1
  593     for k in answerX_Y:
  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
  738         answer.append([x, ppower - x])
  739     return answer
  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