"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/prime.py" (19 Nov 2012, 34853 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 "prime.py" see the Fossies "Dox" file reference documentation.

    1 """
    2 A module for generating primes and testing primality.
    3 """
    4 
    5 import logging
    6 import warnings
    7 import nzmath.arith1 as arith1
    8 import nzmath.gcd as gcd
    9 import nzmath.bigrandom as bigrandom
   10 import nzmath.bigrange as bigrange
   11 import nzmath.poly.array as array_poly
   12 from nzmath.config import GRH
   13 from nzmath.plugins import MATHMODULE as math
   14 
   15 _log = logging.getLogger('nzmath.prime')
   16 _log.setLevel(logging.DEBUG)
   17 
   18 
   19 PRIMES_LE_31 = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31)
   20 PRIMONIAL_31 = arith1.product(PRIMES_LE_31)
   21 
   22 
   23 def trialDivision(n, bound=0):
   24     """
   25     Trial division primality test for an odd natural number.
   26     Optional second argument is a search bound of primes.
   27     If the bound is given and less than the sqaure root of n
   28     and True is returned, it only means there is no prime factor
   29     less than the bound.
   30     """
   31 
   32     if bound:
   33         m = min(bound, arith1.floorsqrt(n))
   34     else:
   35         m = arith1.floorsqrt(n)
   36     for p in bigrange.range(3, m+1, 2):
   37         if not (n % p):
   38             return False
   39     return True
   40 
   41 
   42 def spsp(n, base, s=None, t=None):
   43     """
   44     Strong Pseudo-Prime test.  Optional third and fourth argument
   45     s and t are the numbers such that n-1 = 2**s * t and t is odd.
   46     """
   47     if s is None or t is None:
   48         s, t = arith1.vp(n-1, 2)
   49     z = pow(base, t, n)
   50     if z != 1 and z != n-1:
   51         j = 0
   52         while j < s:
   53             j += 1
   54             z = pow(z, 2, n)
   55             if z == n-1:
   56                 break
   57         else:
   58             return False
   59     return True
   60 
   61 
   62 def miller(n):
   63     """
   64     Miller's primality test.
   65 
   66     This test is valid under GRH.
   67     """
   68     s, t = arith1.vp(n - 1, 2)
   69     # The O-constant 2 by E.Bach
   70     ## ln(n) = log_2(n)/log_2(e) = ln(2) * log_2(n)
   71     ## 2 * ln(2) ** 2 <= 0.960906027836404
   72     bound = min(n - 1, 960906027836404 * (arith1.log(n) + 1) ** 2 // 10 ** 15 + 1)
   73     _log.info("bound: %d" % bound)
   74     for b in range(2, bound):
   75         if not spsp(n, b, s, t):
   76             return False
   77     return True
   78 
   79 def millerRabin(n, times=20):
   80     """
   81     Miller-Rabin pseudo-primality test.  Optional second argument
   82     times (default to 20) is the number of repetition.  The error
   83     probability is at most 4**(-times).
   84     """
   85     s, t = arith1.vp(n - 1, 2)
   86     for _ in range(times):
   87         b = bigrandom.randrange(2, n-1)
   88         if not spsp(n, b, s, t):
   89             return False
   90     return True
   91 
   92 
   93 def bigprimeq(z):
   94     """
   95     Giving up rigorous proof of primality, return True for a probable
   96     prime.
   97     """
   98     if int(z) != z:
   99         raise ValueError("non-integer for primeq()")
  100     elif z <= 1:
  101         return False
  102     elif gcd.gcd(z, PRIMONIAL_31) > 1:
  103         return (z in PRIMES_LE_31)
  104     return millerRabin(z)
  105 
  106 
  107 def Lucas_chain(n, f, g, x_0, x_1):
  108     """
  109     Given an integer n, two functions f and g, and initial value (x_0, x_1),
  110     compute (x_n, x_{n+1}), where the sequence {x_i} is defined as:
  111       x_{2i} = f(x_i)
  112       x_{2i+1} = g(x_i, x_{i+1})
  113     """
  114     binary = arith1.expand(n, 2)
  115     u = x_0
  116     v = x_1
  117     while binary:
  118         if 1 == binary.pop():
  119             u, v = g(u, v), f(v)
  120         else:
  121             u, v = f(u), g(u, v)
  122     return u, v
  123 
  124 
  125 def _lucas_test_sequence(n, a, b):
  126     """
  127     Return x_0, x_1, x_m, x_{m+1} of Lucas sequence of parameter a, b,
  128     where m = (n - (a**2 - 4*b / n)) >> 1.
  129     """
  130     d = a**2 - 4*b
  131     if (d >= 0 and arith1.floorsqrt(d) ** 2 == d) \
  132     or not(gcd.coprime(n, 2*a*b*d)):
  133         raise ValueError("Choose another parameters.")
  134 
  135     x_0 = 2
  136     inv_b = arith1.inverse(b, n)
  137     x_1 = ((a ** 2) * inv_b - 2) % n
  138 
  139     # Chain functions
  140     def even_step(u):
  141         """
  142         'double' u.
  143         """
  144         return (u**2 - x_0) % n
  145 
  146     def odd_step(u, v):
  147         """
  148         'add' u and v.
  149         """
  150         return (u*v - x_1) % n
  151 
  152     m = (n - arith1.legendre(d, n)) >> 1
  153     x_m, x_mplus1 = Lucas_chain(m, even_step, odd_step, x_0, x_1)
  154 
  155     return x_0, x_1, x_m, x_mplus1
  156 
  157 
  158 def lpsp(n, a, b):
  159     """
  160     Lucas test.
  161     Return True if n is a Lucas pseudoprime of parameters a, b,
  162     i.e. with respect to x**2-a*x+b.
  163     """
  164     x_0, x_1, x_m, x_mplus1 = _lucas_test_sequence(n, a, b)
  165 
  166     return (x_1 * x_m - x_0 * x_mplus1) % n == 0
  167 
  168 
  169 def fpsp(n, a, b):
  170     """
  171     Frobenius test.
  172     Return True if n is a Frobenius pseudoprime of parameters a, b,
  173     i.e. with respect to x**2-a*x+b.
  174     """
  175     x_0, x_1, x_m, x_mplus1 = _lucas_test_sequence(n, a, b)
  176 
  177     if (x_1 * x_m - x_0 * x_mplus1) % n == 0:
  178         euler_pow = pow(b, (n-1) >> 1, n)
  179         return (euler_pow * x_m) % n == 2
  180     else:
  181         return False
  182 
  183 
  184 def by_primitive_root(n, divisors):
  185     """
  186     Return True iff n is prime.
  187 
  188     The method proves the primality of n by existence of a primitive
  189     root.  It requires a sequence of all prime divisors of n - 1.
  190 
  191     Lehmer,D.H., Tests for primality by the converse of Fermat's
  192     theorem, Bull.Amer.Math.Soc, Vol.33, pp.327-340, 1927.
  193     """
  194     m_order = n - 1
  195     primes = tuple(p for p in divisors if p > 2)
  196     for g in bigrange.range(2, n):
  197         jacobi = arith1.legendre(g, n)
  198         if jacobi == 0:
  199             return False
  200         elif jacobi == 1:
  201             continue
  202         if pow(g, m_order, n) != 1:
  203             return False
  204         if all(pow(g, m_order // p, n) != 1 for p in primes):
  205             return True
  206     return True # only when n=2, flow reaches this line
  207 
  208 
  209 def full_euler(n, divisors):
  210     """
  211     Return True iff n is prime.
  212 
  213     The method proves the primality of n by the equality
  214     phi(n) = n - 1, where phi denotes the Euler totient.
  215     It requires a sequence of all prime divisors of n - 1.
  216 
  217     Brillhart,J. & Selfridge,J.L., Some Factorizations of $2^n\pm 1$
  218     and Related Results, Math.Comp., Vol.21, 87-96, 1967
  219     """
  220     m_order = n - 1
  221     primes = set(divisors)
  222     for g in bigrange.range(2, n):
  223         if pow(g, m_order, n) != 1:
  224             return False
  225         if 2 in primes:
  226             jacobi = arith1.legendre(g, n)
  227             if jacobi == 0:
  228                 return False
  229             elif jacobi == -1:
  230                 primes.remove(2)
  231         satisfied = [p for p in primes if p > 2 and pow(g, m_order // p, n) != 1]
  232         if satisfied:
  233             primes.difference_update(satisfied)
  234         if not primes:
  235             return True
  236     return True # only when n=2, flow reaches this line
  237 
  238 
  239 def prime(s):
  240     """
  241     prime(n) returns the n-th prime number.
  242     """
  243     if s != int(s):
  244         raise ValueError("non-integer for prime()")
  245     elif s <= 0:
  246         raise ValueError("non-positive-integer for prime()")
  247     i = 1
  248     for p in generator():
  249         if i == s:
  250             return p
  251         i += 1
  252     # The following line should not be reached:
  253     raise ValueError("Too big number %d for prime(i)." % s)
  254 
  255 
  256 def generator():
  257     """
  258     Generate primes from 2 to infinity.
  259     """
  260     yield 2
  261     yield 3
  262     yield 5
  263     coprimeTo30 = (7, 11, 13, 17, 19, 23, 29, 31)
  264     times30 = 0
  265     while True:
  266         for i in coprimeTo30:
  267             if primeq(i + times30):
  268                 yield i + times30
  269         times30 += 30
  270 
  271 
  272 def generator_eratosthenes(n):
  273     """
  274     Generate primes up to n (inclusive) using Eratosthenes sieve.
  275     """
  276     if n < 2:
  277         raise StopIteration
  278 
  279     yield 2
  280     if n <= 2:
  281         return
  282 
  283     yield 3
  284 
  285     # make list for sieve
  286     sieve_len_max = (n+1) >> 1
  287     sieve = [True, False, True]
  288     sieve_len = 3
  289     k = 5
  290     i = 2
  291     while sieve_len < sieve_len_max:
  292         if sieve[i]:
  293             yield k
  294             sieve_len *= k
  295             if sieve_len_max < sieve_len:
  296                 sieve_len //= k
  297                 # adjust sieve list length
  298                 sieve *= sieve_len_max // sieve_len
  299                 sieve += sieve[:(sieve_len_max - len(sieve))]
  300                 sieve_len = sieve_len_max
  301             else:
  302                 sieve = sieve * k
  303             for j in range(i, sieve_len, k):
  304                 sieve[j] = False
  305         k += 2
  306         i += 1
  307 
  308     # sieve
  309     limit = arith1.floorsqrt(n)
  310     while k <= limit:
  311         if sieve[i]:
  312             yield k
  313             j = (k ** 2 - 1) >> 1
  314             while j < sieve_len_max:
  315                 sieve[j] = False
  316                 j += k
  317         k += 2
  318         i += 1
  319 
  320     # output result
  321     limit = (n - 1) >> 1
  322     while i <= limit:
  323         if sieve[i]:
  324             yield 2 * i + 1
  325         i += 1
  326 
  327 
  328 def nextPrime(n):
  329     """
  330     Return the smallest prime bigger than the given integer.
  331     """
  332     if n <= 1:
  333         return 2
  334     if n == 2:
  335         return 3
  336     n += (1 + (n & 1)) # make n be odd.
  337     while not primeq(n):
  338         n += 2
  339     return n
  340 
  341 
  342 def randPrime(n):
  343     """
  344     Return a random n-digits prime
  345     """
  346     if n <= 0:
  347         raise ValueError("input number must be natural number")
  348 
  349     if n == 1:
  350         return bigrandom.map_choice(lambda i: (2, 3, 5, 7)[i], 4)
  351 
  352     p = bigrandom.randrange(10**(n-1)+1, 10**n, 2)
  353     while not primeq(p):
  354         p += 2
  355     if p < 10**n:
  356         return p
  357 
  358     # in very rare case or n is too small case,
  359     # search continues from the lower bound.
  360     p = 10**(n-1) + 1
  361     while not primeq(p):
  362         p += 2
  363     return p
  364 
  365 
  366 def smallSpsp(n, s=None, t=None):
  367     """
  368     4 spsp tests are sufficient to determine whether an integer less
  369     than 10**12 is prime or not.  Optional third and fourth argument
  370     s and t are the numbers such that n - 1 = 2**s * t and t is odd.
  371     """
  372     if s is None or t is None:
  373         s, t = arith1.vp(n - 1, 2)
  374     for p in (2, 13, 23, 1662803):
  375         if not spsp(n, p, s, t):
  376             return False
  377     return True
  378 
  379 
  380 def primeq(n):
  381     """
  382     A convenient function for primatilty test. It uses one of
  383     trialDivision, smallSpsp or apr depending on the size of n.
  384     """
  385     if int(n) != n:
  386         raise ValueError("non-integer for primeq()")
  387     if n <= 1:
  388         return False
  389 
  390     if gcd.gcd(n, PRIMONIAL_31) > 1:
  391         return (n in PRIMES_LE_31)
  392     if n < 2000000:
  393         return trialDivision(n)
  394     if not smallSpsp(n):
  395         return False
  396     if n < 10 ** 12:
  397         return True
  398     if not GRH:
  399         return apr(n)
  400     else:
  401         return miller(n)
  402 
  403 
  404 def primonial(p):
  405     """
  406     Return 2*3*...*p for given prime p.
  407     """
  408     return arith1.product(generator_eratosthenes(p))
  409 
  410 
  411 # defs for APR algorithm
  412 
  413 def _isprime(n, pdivisors=None):
  414     """
  415     Return True iff n is prime.
  416 
  417     The check is done without APR.
  418     Assume that n is very small (less than 10**12) or
  419     prime factorization of n - 1 is known (prime divisors are passed to
  420     the optional argument pdivisors as a sequence).
  421     """
  422     if gcd.gcd(n, PRIMONIAL_31) > 1:
  423         return (n in PRIMES_LE_31)
  424     elif n < 10 ** 12:
  425         # 1369 == 37**2
  426         # 1662803 is the only prime base in smallSpsp which has not checked
  427         return n < 1369 or n == 1662803 or smallSpsp(n)
  428     else:
  429         return full_euler(n, pdivisors)
  430 
  431 def properDivisors(n):
  432     """
  433     Return proper divisors of n (divisors of n excluding 1 and n).
  434 
  435     It is only useful for a product of small primes.
  436     One can use FactoredInteger.proper_divisors() as well.
  437 
  438     DEPRECATION: this function will be removed in version 1.3.
  439     """
  440     warnings.warn(DeprecationWarning("properDivisor is deprecated"))
  441 
  442     if n in (1, 2, 3, 5, 7, 11, 13, 17, 19, 23):
  443         return []
  444     else:
  445         l = [1]
  446         for p, e in _factor(n):
  447             for j in range(1, e + 1):
  448                 l += [k*pow(p, j) for k in l if k % p != 0]
  449         l.remove(1)
  450         l.remove(n)
  451         l.sort()
  452         return l
  453 
  454 def _factor(n, bound=0):
  455     """
  456     Trial division factorization for a natural number.
  457     Optional second argument is a search bound of primes.
  458     If the bound is given and less than the sqaure root of n,
  459     result is not proved to be a prime factorization.
  460     """
  461     factors = []
  462     if not (n & 1):
  463         v2, n = arith1.vp(n, 2)
  464         factors.append((2, v2))
  465     m = _calc_bound(n, bound)
  466     p = 3
  467     while p <= m:
  468         if not (n % p):
  469             v, n = arith1.vp(n, p)
  470             factors.append((p, v))
  471             m = _calc_bound(n, bound)
  472         p += 2
  473     if n > 1:
  474         factors.append((n, 1))
  475     return factors
  476 
  477 def _calc_bound(n, bound=0):
  478     if bound:
  479         m = min((bound, arith1.floorsqrt(n)))
  480     else:
  481         m = arith1.floorsqrt(n)
  482     return m
  483 
  484 def primitive_root(p):
  485     """
  486     Return a primitive root of p.
  487     """
  488     pd = FactoredInteger(p - 1).proper_divisors()
  489     for i in bigrange.range(2, p):
  490         for d in pd:
  491             if pow(i, (p - 1)//d, p) == 1:
  492                 break
  493         else:
  494             return i
  495 
  496 
  497 class Zeta(object):
  498     """
  499     Represent linear combinations of roots of unity.
  500     """
  501     def __init__(self, size, pos=None, val=1):
  502         self.size = size
  503         self.z = [0] * self.size
  504         if pos is not None:
  505             self.z[pos % self.size] = val
  506 
  507     def __add__(self, other):
  508         if self.size == other.size:
  509             m = self.size
  510             zr_a = Zeta(m)
  511             for i in range(m):
  512                 zr_a.z[i] = self.z[i] + other.z[i]
  513             return zr_a
  514         else:
  515             m = gcd.lcm(self.size, other.size)
  516             return self.promote(m) + other.promote(m)
  517 
  518     def __mul__(self, other):
  519         if not isinstance(other, Zeta):
  520             zr_m = Zeta(self.size)
  521             zr_m.z = [x*other for x in self.z]
  522             return zr_m
  523         elif self.size == other.size:
  524             zr_m = Zeta(self.size)
  525             other = +other
  526             for k in range(other.size):
  527                 if not other.z[k]:
  528                     continue
  529                 elif other.z[k] == 1:
  530                     zr_m = zr_m + (self<<k)
  531                 else:
  532                     zr_m = zr_m + (self<<k)*other.z[k]
  533             return zr_m
  534         else:
  535             m = gcd.lcm(self.size, other.size)
  536             return self.promote(m)*other.promote(m)
  537 
  538     __rmul__ = __mul__
  539 
  540     def __lshift__(self, offset):
  541         """
  542         The name is shift but the meaning of function is rotation.
  543         """
  544         new = Zeta(self.size)
  545         new.z = self.z[-offset:] + self.z[:-offset]
  546         return new
  547 
  548     def __pow__(self, e, mod=0):
  549         if mod:
  550             return self._pow_mod(e, mod)
  551 
  552         r = Zeta(self.size, 0)
  553         if e == 0:
  554             return r
  555         z = +self
  556         while True:
  557             if e & 1:
  558                 r = z*r
  559                 if e == 1:
  560                     return r
  561             e //= 2
  562             z = z._square()
  563 
  564     def _pow_mod(self, e, mod):
  565         r = Zeta(self.size, 0)
  566         if e == 0:
  567             return r
  568         z = self % mod
  569         while True:
  570             if e & 1:
  571                 r = z * r % mod
  572                 if e == 1:
  573                     return r
  574             e //= 2
  575             z = z._square_mod(mod)
  576 
  577     def _square(self):
  578         return self * self
  579 
  580     def _square_mod(self, mod):
  581         return self * self % mod
  582 
  583     def __pos__(self):
  584         m = self.size
  585         z_p = Zeta(m)
  586         if m & 1 == 0:
  587             mp = m >> 1
  588             for i in range(mp):
  589                 if self.z[i] > self.z[i+mp]:
  590                     z_p.z[i] = self.z[i] - self.z[i+mp]
  591                 else:
  592                     z_p.z[i+mp] = self.z[i+mp] - self.z[i]
  593         else:
  594             p = 3
  595             while m % p:
  596                 p += 2
  597             mp = m // p
  598             for i in range(mp):
  599                 mini = self.z[i]
  600                 for j in range(mp + i, m, mp):
  601                     if mini > self.z[j]:
  602                         mini = self.z[j]
  603                 for j in range(i, m, mp):
  604                     z_p.z[j] = self.z[j] - mini
  605         return z_p
  606 
  607     def __mod__(self, mod):
  608         """
  609         Return a new Zeta instance modulo 'mod'.
  610 
  611         All entries are reduced to the least absolute residues.
  612         """
  613         new = Zeta(self.size)
  614         half = (mod - 1) >> 1 # -1 to make 1 (mod 2) be 1, not -1.
  615         new.z = [(x + half) % mod - half for x in self.z]
  616         return new
  617 
  618     def __setitem__(self, position, value):
  619         self.z[position % self.size] = value
  620 
  621     def __getitem__(self, position):
  622         return self.z[position % self.size]
  623 
  624     def promote(self, size):
  625         if size == self.size:
  626             return +self
  627         new = Zeta(size)
  628         r = size // self.size
  629         for i in range(self.size):
  630             new.z[i*r] = self.z[i]
  631         return new
  632 
  633     def __len__(self):
  634         return self.size
  635 
  636     def __eq__(self, other):
  637         for i in range(self.size):
  638             if self.z[i] != other.z[i]:
  639                 return False
  640         return True
  641 
  642     def __hash__(self):
  643         return sum([hash(self.z[i]) for i in range(1, self.size)])
  644 
  645     def weight(self):
  646         """
  647         Return the number of nonzero entries.
  648         """
  649         return len(filter(None, self.z))
  650 
  651     def mass(self):
  652         """
  653         Return the sum of all entries.
  654         """
  655         return sum(self.z)
  656 
  657 
  658 class FactoredInteger(object):
  659     """
  660     Integers with factorization information.
  661     """
  662     def __init__(self, integer, factors=None):
  663         """
  664         FactoredInteger(integer [, factors])
  665 
  666         If factors is given, it is a dict of type {prime:exponent}
  667         and the product of prime**exponent is equal to the integer.
  668         Otherwise, factorization is carried out in initialization.
  669         """
  670         self.integer = int(integer)
  671         if factors is None:
  672             self.factors = dict(_factor(self.integer))
  673         else:
  674             self.factors = dict(factors)
  675 
  676     @classmethod
  677     def from_partial_factorization(cls, integer, partial):
  678         """
  679         Construct a new FactoredInteger object from partial
  680         factorization information given as dict of type
  681         {prime:exponent}.
  682         """
  683         partial_factor = 1
  684         for p, e in partial.iteritems():
  685             partial_factor *= p**e
  686         assert not integer % partial_factor, "wrong factorization"
  687         return cls(integer // partial_factor) * cls(partial_factor, partial)
  688 
  689     def __iter__(self):
  690         """
  691         Default iterator
  692         """
  693         return self.factors.iteritems()
  694 
  695     def __mul__(self, other):
  696         if isinstance(other, FactoredInteger):
  697             integer = self.integer * other.integer
  698             new_factors = self.factors.copy()
  699             for p in other.factors:
  700                 new_factors[p] = new_factors.get(p, 0) + other.factors[p]
  701             return self.__class__(integer, new_factors)
  702         else:
  703             return self * FactoredInteger(other)
  704 
  705     __rmul__ = __mul__
  706 
  707     def __pow__(self, other):
  708         new_integer = self.integer**other
  709         new_factors = {}
  710         for p in self.factors:
  711             new_factors[p] = self.factors[p] * other
  712         return self.__class__(new_integer, new_factors)
  713 
  714     def __pos__(self):
  715         return self.copy()
  716 
  717     def __str__(self):
  718         return str(self.integer)
  719 
  720     def __eq__(self, other):
  721         try:
  722             return self.integer == other.integer
  723         except AttributeError:
  724             return self.integer == int(other)
  725 
  726     def __hash__(self):
  727         return hash(self.integer)
  728 
  729     def __ne__(self, other):
  730         return not self.__eq__(other)
  731 
  732     def __le__(self, other):
  733         try:
  734             return self.integer <= other.integer
  735         except AttributeError:
  736             return self.integer <= int(other)
  737 
  738     def __lt__(self, other):
  739         try:
  740             return self.integer < other.integer
  741         except AttributeError:
  742             return self.integer < int(other)
  743 
  744     def __gt__(self, other):
  745         try:
  746             return self.integer > other.integer
  747         except AttributeError:
  748             return self.integer > int(other)
  749 
  750     def __ge__(self, other):
  751         try:
  752             return self.integer >= other.integer
  753         except AttributeError:
  754             return self.integer >= int(other)
  755 
  756     def __long__(self):
  757         return int(self.integer)
  758 
  759     __int__ = __long__
  760 
  761     def __mod__(self, other):
  762         # maybe you want self.is_divisible_by(other)
  763         try:
  764             if other.integer in self.factors:
  765                 return 0
  766             return self.integer % other.integer
  767         except AttributeError:
  768             if int(other) in self.factors:
  769                 return 0
  770             return self.integer % int(other)
  771 
  772     def __rfloordiv__(self, other):
  773         # assume other is an integer.
  774         return other // self.integer
  775 
  776     def copy(self):
  777         return self.__class__(self.integer, self.factors.copy())
  778 
  779     def is_divisible_by(self, other):
  780         """
  781         Return True if other divides self.
  782         """
  783         if int(other) in self.factors:
  784             # other is prime and divides
  785             return True
  786         return not self.integer % int(other)
  787 
  788     def exact_division(self, other):
  789         """
  790         Divide by a factor.
  791         """
  792         divisor = int(other)
  793         quotient = self.copy()
  794         if divisor in quotient.factors:
  795             if quotient.factors[divisor] == 1:
  796                 del quotient.factors[divisor]
  797             else:
  798                 quotient.factors[divisor] -= 1
  799         elif not isinstance(other, FactoredInteger):
  800             dividing = divisor
  801             for p, e in self.factors.iteritems():
  802                 while not dividing % p:
  803                     dividing //= p
  804                     if quotient.factors[p] == 1:
  805                         del quotient.factors[p]
  806                         assert dividing % p, dividing
  807                     else:
  808                         quotient.factors[p] -= 1
  809                 if dividing == 1:
  810                     break
  811             assert dividing == 1
  812         else:
  813             for p, e in other.factors.iteritems():
  814                 assert p in quotient.factors and quotient.factors[p] >= e
  815                 if quotient.factors[p] == e:
  816                     del quotient.factors[p]
  817                 else:
  818                     quotient.factors[p] -= e
  819         quotient.integer //= divisor
  820         return quotient
  821 
  822     # maybe this is what you want, isn't it?
  823     __floordiv__ = exact_division
  824 
  825     def divisors(self):
  826         """
  827         Return all divisors.
  828         """
  829         divs = [FactoredInteger(1)]
  830         for p, e in self.factors.iteritems():
  831             q = FactoredInteger(1)
  832             pcoprimes = list(divs)
  833             for j in range(1, e + 1):
  834                 q *= FactoredInteger(p, {p:1})
  835                 divs += [k * q for k in pcoprimes]
  836         return divs
  837 
  838     def proper_divisors(self):
  839         """
  840         Return the proper divisors (divisors of n excluding 1 and n).
  841         """
  842         return self.divisors()[1:-1]
  843 
  844     def prime_divisors(self):
  845         """
  846         Return the list of primes that divides the number.
  847         """
  848         return self.factors.keys()
  849 
  850 
  851 class TestPrime(object):
  852     primes = PRIMES_LE_31
  853     primecache = set(primes)
  854 
  855     def __init__(self, t=12):
  856         if isinstance(t, (int, long)):
  857             self.t = FactoredInteger(t)
  858         else:
  859             assert isinstance(t, FactoredInteger)
  860             self.t = t
  861         powerof2 = self.t.factors[2] + 2
  862         self.et = FactoredInteger(2 ** powerof2, {2:powerof2})
  863         for d in self.t.divisors():
  864             # d is an instance of FactoredInteger
  865             p = d.integer + 1
  866             if p & 1 and (p in self.primecache or _isprime(p, d.factors)):
  867                 self.et = self.et * FactoredInteger(p, {p:1})
  868                 if p in self.t.factors:
  869                     e = self.t.factors[p]
  870                     self.et = self.et * FactoredInteger(p**e, {p:e})
  871                 self.primecache.add(p)
  872 
  873     def next(self):
  874         eu = []
  875         for p in self.primes:
  876             if p in self.t.factors:
  877                 eu.append((p - 1) * p**(self.t.factors[p] - 1))
  878             else:
  879                 eu.append(p - 1)
  880                 break
  881         p = self.primes[eu.index(min(eu))]
  882         return self.__class__(self.t * FactoredInteger(p, {p:1}))
  883 
  884 
  885 class Status:
  886     """
  887     status collector for apr.
  888     """
  889     def __init__(self):
  890         self.d = {}
  891 
  892     def yet(self, key):
  893         """
  894         set key's status be 'yet'.
  895         """
  896         self.d[key] = 0
  897 
  898     def done(self, key):
  899         """
  900         set key's status be 'done'.
  901         """
  902         self.d[key] = 1
  903 
  904     def yet_keys(self):
  905         """
  906         Return keys whose stati are 'yet'.
  907         """
  908         return [k for k in self.d if not self.d[k]]
  909 
  910     def isDone(self, key):
  911         """
  912         Return whether key's status is 'done' or not.
  913         """
  914         return self.d[key]
  915 
  916     def subodd(self, p, q, n, J):
  917         """
  918         Return the sub result for odd key 'p'.
  919         If it is True, the status of 'p' is flipped to 'done'.
  920         """
  921         s = J.get(1, p, q)
  922         Jpq = J.get(1, p, q)
  923         m = s.size
  924         for x in range(2, m):
  925             if x % p == 0:
  926                 continue
  927             sx = Zeta(m)
  928             i = x
  929             j = 1
  930             while i > 0:
  931                 sx[j] = Jpq[i]
  932                 i = (i + x) % m
  933                 j += 1
  934             sx[0] = Jpq[0]
  935             sx = pow(sx, x, n)
  936             s = s*sx % n
  937         s = pow(s, n//m, n)
  938         r = n % m
  939         t = 1
  940         for x in range(1, m):
  941             if x % p == 0:
  942                 continue
  943             c = (r*x) // m
  944             if c:
  945                 tx = Zeta(m)
  946                 i = x
  947                 j = 1
  948                 while i > 0:
  949                     tx[j] = Jpq[i]
  950                     i = (i + x) % m
  951                     j += 1
  952                 tx[0] = Jpq[0]
  953                 tx = pow(tx, c, n)
  954                 t = t*tx % n
  955         s = +(t*s % n)
  956         if s.weight() == 1 and s.mass() == 1:
  957             for i in range(1, m):
  958                 if gcd.gcd(m, s.z.index(1)) == 1:
  959                     self.done(p)
  960                 return True
  961         return False
  962 
  963     def sub8(self, q, k, n, J):
  964         s = J.get(3, q)
  965         J3 = J.get(3, q)
  966         m = len(s)
  967         sx_z = {1:s}
  968         x = 3
  969         step = 2
  970         while m > x:
  971             z_4b = Zeta(m)
  972             i = x
  973             j = 1
  974             while i != 0:
  975                 z_4b[j] = J3[i]
  976                 i = (i + x) % m
  977                 j += 1
  978             z_4b[0] = J3[0]
  979             sx_z[x] = z_4b
  980             s = pow(sx_z[x], x, n) * s
  981             step = 8 - step
  982             x += step
  983 
  984         s = pow(s, n//m, n)
  985 
  986         r = n % m
  987         step = 2
  988         x = 3
  989         while m > x:
  990             c = r*x
  991             if c > m:
  992                 s = pow(sx_z[x], c//m, n) * s
  993             step = 8 - step
  994             x += step
  995         r = r & 7
  996         if r == 5 or r == 7:
  997             s = J.get(2, q).promote(m) * s
  998         s = +(s % n)
  999 
 1000         if s.weight() == 1 and s.mass() == 1:
 1001             if gcd.gcd(m, s.z.index(1)) == 1 and pow(q, (n-1) >> 1, n) == n-1:
 1002                 self.done(2)
 1003             return True
 1004         elif s.weight() == 1 and s.mass() == n-1:
 1005             if gcd.gcd(m, s.z.index(n-1)) == 1 and pow(q, (n-1) >> 1, n) == n-1:
 1006                 self.done(2)
 1007             return True
 1008         return False
 1009 
 1010     def sub4(self, q, n, J):
 1011         j2 = J.get(1, 2, q)**2
 1012         s = q*j2 % n
 1013         s = pow(s, n >> 2, n)
 1014         if n & 3 == 3:
 1015             s = s*j2 % n
 1016         s = +(s % n)
 1017         if s.weight() == 1 and s.mass() == 1:
 1018             i = s.z.index(1)
 1019             if (i == 1 or i == 3) and pow(q, (n-1) >> 1, n) == n-1:
 1020                 self.done(2)
 1021             return True
 1022         return False
 1023 
 1024     def sub2(self, q, n):
 1025         s = pow(n-q, (n-1) >> 1, n)
 1026         if s == n-1:
 1027             if n & 3 == 1:
 1028                 self.done(2)
 1029         elif s != 1:
 1030             return False
 1031         return True
 1032 
 1033     def subrest(self, p, n, et, J, ub=200):
 1034         if p == 2:
 1035             q = 5
 1036             while q < 2*ub + 5:
 1037                 q += 2
 1038                 if not _isprime(q) or et % q == 0:
 1039                     continue
 1040                 if n % q == 0:
 1041                     _log.info("%s divides %s.\n" % (q, n))
 1042                     return False
 1043                 k = arith1.vp(q-1, 2)[0]
 1044                 if k == 1:
 1045                     if n & 3 == 1 and not self.sub2(q, n):
 1046                         return False
 1047                 elif k == 2:
 1048                     if not self.sub4(q, n, J):
 1049                         return False
 1050                 else:
 1051                     if not self.sub8(q, k, n, J):
 1052                         return False
 1053                 if self.isDone(p):
 1054                     return True
 1055             else:
 1056                 raise ImplementLimit("limit")
 1057         else:
 1058             step = p*2
 1059             q = 1
 1060             while q < step*ub + 1:
 1061                 q += step
 1062                 if not _isprime(q) or et % q == 0:
 1063                     continue
 1064                 if n % q == 0:
 1065                     _log.info("%s divides %s.\n" % (q, n))
 1066                     return False
 1067                 if not self.subodd(p, q, n, J):
 1068                     return False
 1069                 if self.isDone(p):
 1070                     return True
 1071             else:
 1072                 raise ImplementLimit("limit")
 1073 
 1074 
 1075 class JacobiSum:
 1076     """
 1077     A repository of the Jacobi sums.
 1078     """
 1079     def __init__(self):
 1080         self.shelve = {}
 1081 
 1082     def get(self, group, p, q=None):
 1083         if q:
 1084             assert group == 1
 1085             if (group, p, q) not in self.shelve:
 1086                 self.make(q)
 1087             return self.shelve[group, p, q]
 1088         else:
 1089             assert group == 2 or group == 3
 1090             if (group, p) not in self.shelve:
 1091                 self.make(p)
 1092             return self.shelve[group, p]
 1093 
 1094     def make(self, q):
 1095         fx = self.makefx(q)
 1096         qpred = q - 1
 1097         qt = _factor(qpred)
 1098         qt2 = [k for (p, k) in qt if p == 2][0]
 1099         k, pk = qt2, 2**qt2
 1100         if k >= 2:
 1101             J2q = Zeta(pk, 1 + fx[1])
 1102             for j in range(2, qpred):
 1103                 J2q[j + fx[j]] = J2q[j + fx[j]] + 1
 1104             self.shelve[1, 2, q] = +J2q
 1105             if k >= 3:
 1106                 J2 = Zeta(8, 3 + fx[1])
 1107                 J3 = Zeta(pk, 2 + fx[1])
 1108                 for j in range(2, qpred):
 1109                     J2[j*3 + fx[j]] = J2[j*3 + fx[j]] + 1
 1110                     J3[j*2 + fx[j]] = J3[j*2 + fx[j]] + 1
 1111                 self.shelve[3, q] = +(self.shelve[1, 2, q]*J3)
 1112                 self.shelve[2, q] = +(J2*J2)
 1113         else:
 1114             self.shelve[1, 2, q] = 1
 1115         for (p, k) in qt:
 1116             if p == 2:
 1117                 continue
 1118             pk = p**k
 1119             Jpq = Zeta(pk, 1 + fx[1])
 1120             for j in range(2, qpred):
 1121                 Jpq[j + fx[j]] = Jpq[j + fx[j]] + 1
 1122             self.shelve[1, p, q] = +Jpq
 1123         del fx
 1124 
 1125     @staticmethod
 1126     def makefx(q):
 1127         """
 1128         Return a dict called 'fx'.
 1129         The dict stores the information that fx[i] == j iff
 1130           g**i + g**j = 1 mod q
 1131         for g, a primitive root of the prime q.
 1132         """
 1133         g = primitive_root(q)
 1134         qpred = q - 1
 1135         qd2 = qpred >> 1
 1136         g_mf = [0, g]
 1137         for _ in range(2, qpred):
 1138             g_mf.append((g_mf[-1]*g) % q)
 1139         fx = {}
 1140         for i in range(1, qd2):
 1141             if i in fx:
 1142                 continue
 1143             # search j s.t. g**j + g**i = 1 mod q
 1144             j = g_mf.index(q + 1 - g_mf[i])
 1145             fx[i] = j
 1146             fx[j] = i
 1147             fx[qpred - i] = (j - i + qd2) % qpred
 1148             fx[fx[qpred - i]] = qpred - i
 1149             fx[qpred - j] = (i - j + qd2) % qpred
 1150             fx[fx[qpred - j]] = qpred - j
 1151         del g_mf
 1152         return fx
 1153 
 1154 
 1155 def apr(n):
 1156     """
 1157     apr is the main function for Adleman-Pomerance-Rumery primality test.
 1158     Assuming n has no prime factors less than 32.
 1159     Assuming n is spsp for several bases.
 1160     """
 1161     L = Status()
 1162 
 1163     rb = arith1.floorsqrt(n) + 1
 1164     el = TestPrime()
 1165     while el.et <= rb:
 1166         el = el.next()
 1167 
 1168     plist = el.t.factors.keys()
 1169     plist.remove(2)
 1170     L.yet(2)
 1171     for p in plist:
 1172         if pow(n, p-1, p*p) != 1:
 1173             L.done(p)
 1174         else:
 1175             L.yet(p)
 1176     qlist = el.et.factors.keys()
 1177     qlist.remove(2)
 1178     J = JacobiSum()
 1179     for q in qlist:
 1180         for p in plist:
 1181             if (q-1) % p != 0:
 1182                 continue
 1183             if not L.subodd(p, q, n, J):
 1184                 return False
 1185         k = arith1.vp(q-1, 2)[0]
 1186         if k == 1:
 1187             if not L.sub2(q, n):
 1188                 return False
 1189         elif k == 2:
 1190             if not L.sub4(q, n, J):
 1191                 return False
 1192         else:
 1193             if not L.sub8(q, k, n, J):
 1194                 return False
 1195     for p in L.yet_keys():
 1196         if not L.subrest(p, n, el.et, J):
 1197             return False
 1198     r = int(n)
 1199     for _ in bigrange.range(1, el.t.integer):
 1200         r = (r*n) % el.et.integer
 1201         if n % r == 0 and r != 1 and r != n:
 1202             _log.info("%s divides %s.\n" %(r, n))
 1203             return False
 1204     return True
 1205 
 1206 
 1207 class ImplementLimit (Exception):
 1208     """
 1209     Exception throwed when the execution hits the implementation limit.
 1210     """
 1211     pass
 1212 
 1213 
 1214 # AKS
 1215 def aks(n):
 1216     """
 1217     AKS ( Cyclotomic Congruence Test ) primality test for a natural number.
 1218     
 1219     Input: a natural number n ( n > 1 ).
 1220     Output: n is prime => return True
 1221             n is not prime => return False.
 1222     """
 1223     import nzmath.multiplicative as multiplicative
 1224 
 1225     if arith1.powerDetection(n)[1] != 1: #Power Detection
 1226         return False
 1227 
 1228     lg = math.log(n, 2)
 1229     k = int(lg * lg)
 1230 
 1231     start = 3
 1232     while 1:
 1233         d = gcd.gcd(start, n)
 1234         if 1 < d < n:
 1235             return False
 1236         x = n % start
 1237         N = x
 1238         for i in range(1, k + 1):
 1239             if N == 1:
 1240                 break
 1241             N = (N * x) % start
 1242         if i == k:
 1243             r = start
 1244             break
 1245         start += 1
 1246     d = gcd.gcd(r, n)
 1247     if 1 < d < n:
 1248         return False
 1249     if n <= r:
 1250         return True
 1251 
 1252     e = multiplicative.euler(r) #Cyclotomic Conguence
 1253     e = math.sqrt(e)
 1254     e = int(e * lg)
 1255     for b in range(1, e + 1):
 1256         f = array_poly.ArrayPolyMod([b, 1], n)
 1257         total = array_poly.ArrayPolyMod([1], n)
 1258         count = n
 1259         while count > 0:
 1260             if count & 1:
 1261                 total = total * f
 1262                 total = _aks_mod(total, r)
 1263             f = f.power()
 1264             f = _aks_mod(f, r)
 1265             count = count >> 1
 1266         total_poly = total.coefficients_to_dict()
 1267         if total_poly != {0:b, n % r:1}:
 1268             return False
 1269     return True
 1270 
 1271 def _aks_mod(polynomial, r):
 1272     """
 1273     This function is used in aks.
 1274     polynomial modulo (x^r - 1)
 1275     """
 1276     total = polynomial.coefficients[:r]
 1277     aks_mod = polynomial.coefficients[r:]
 1278     while len(aks_mod) - 1 >= r:
 1279         for i in range(r):
 1280             total[i] += aks_mod[i]
 1281         aks_mod = aks_mod[r:]
 1282     for i in range(len(aks_mod)):
 1283         total[i] += aks_mod[i]
 1284     return array_poly.ArrayPolyMod(total, polynomial.mod)