## "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})
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)