NZMATH  1.2.0
About: NZMATH is a Python based number theory oriented calculation system.
  Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

prime.py
Go to the documentation of this file.
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 
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 
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 
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 
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 
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)
nzmath.prime.FactoredInteger.exact_division
def exact_division(self, other)
Definition: prime.py:788
nzmath.prime.Zeta.__getitem__
def __getitem__(self, position)
Definition: prime.py:621
nzmath.bigrandom
Definition: bigrandom.py:1
nzmath.prime.Zeta.__pos__
def __pos__(self)
Definition: prime.py:583
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.prime.Status.isDone
def isDone(self, key)
Definition: prime.py:910
nzmath.prime.fpsp
def fpsp(n, a, b)
Definition: prime.py:169
nzmath.prime.FactoredInteger.__pos__
def __pos__(self)
Definition: prime.py:714
nzmath.prime.FactoredInteger.__gt__
def __gt__(self, other)
Definition: prime.py:744
nzmath.prime.FactoredInteger.__mod__
def __mod__(self, other)
Definition: prime.py:761
nzmath.prime.primonial
def primonial(p)
Definition: prime.py:404
nzmath.prime.Zeta
Definition: prime.py:497
nzmath.prime.Zeta.z
z
Definition: prime.py:503
nzmath.prime._factor
def _factor(n, bound=0)
Definition: prime.py:454
nzmath.prime.Zeta.promote
def promote(self, size)
Definition: prime.py:624
nzmath.prime.millerRabin
def millerRabin(n, times=20)
Definition: prime.py:79
nzmath.prime._calc_bound
def _calc_bound(n, bound=0)
Definition: prime.py:477
nzmath.prime.properDivisors
def properDivisors(n)
Definition: prime.py:431
nzmath.prime.spsp
def spsp(n, base, s=None, t=None)
Definition: prime.py:42
nzmath.prime.FactoredInteger.__le__
def __le__(self, other)
Definition: prime.py:732
nzmath.prime.FactoredInteger.from_partial_factorization
def from_partial_factorization(cls, integer, partial)
Definition: prime.py:677
nzmath.prime.Zeta.__add__
def __add__(self, other)
Definition: prime.py:507
nzmath.prime.FactoredInteger.__mul__
def __mul__(self, other)
Definition: prime.py:695
nzmath.prime.FactoredInteger.__hash__
def __hash__(self)
Definition: prime.py:726
nzmath.prime.Zeta.__pow__
def __pow__(self, e, mod=0)
Definition: prime.py:548
nzmath.prime.FactoredInteger.is_divisible_by
def is_divisible_by(self, other)
Definition: prime.py:779
nzmath.prime.Status.sub8
def sub8(self, q, k, n, J)
Definition: prime.py:963
nzmath.prime.apr
def apr(n)
Definition: prime.py:1155
nzmath.gcd
Definition: gcd.py:1
nzmath.prime.generator_eratosthenes
def generator_eratosthenes(n)
Definition: prime.py:272
nzmath.prime.full_euler
def full_euler(n, divisors)
Definition: prime.py:209
nzmath.prime.JacobiSum.make
def make(self, q)
Definition: prime.py:1094
nzmath.prime.Lucas_chain
def Lucas_chain(n, f, g, x_0, x_1)
Definition: prime.py:107
nzmath.prime.by_primitive_root
def by_primitive_root(n, divisors)
Definition: prime.py:184
nzmath.prime.Status.subodd
def subodd(self, p, q, n, J)
Definition: prime.py:916
nzmath.prime.Zeta._square
def _square(self)
Definition: prime.py:577
nzmath.prime.FactoredInteger
Definition: prime.py:658
nzmath.prime.lpsp
def lpsp(n, a, b)
Definition: prime.py:158
nzmath.prime.FactoredInteger.factors
factors
Definition: prime.py:672
nzmath.prime.FactoredInteger.divisors
def divisors(self)
Definition: prime.py:825
nzmath.prime.Zeta.__init__
def __init__(self, size, pos=None, val=1)
Definition: prime.py:501
nzmath.prime.randPrime
def randPrime(n)
Definition: prime.py:342
nzmath.prime.FactoredInteger.__lt__
def __lt__(self, other)
Definition: prime.py:738
nzmath.prime.miller
def miller(n)
Definition: prime.py:62
nzmath.prime.FactoredInteger.copy
def copy(self)
Definition: prime.py:776
nzmath.multiplicative
Definition: multiplicative.py:1
nzmath.prime.FactoredInteger.proper_divisors
def proper_divisors(self)
Definition: prime.py:838
nzmath.prime.Zeta.__len__
def __len__(self)
Definition: prime.py:633
nzmath.prime.smallSpsp
def smallSpsp(n, s=None, t=None)
Definition: prime.py:366
nzmath.prime.TestPrime.et
et
Definition: prime.py:862
nzmath.prime.Status.sub2
def sub2(self, q, n)
Definition: prime.py:1024
nzmath.prime.Status
Definition: prime.py:885
nzmath.prime.TestPrime.primecache
primecache
Definition: prime.py:853
nzmath.prime.Zeta.weight
def weight(self)
Definition: prime.py:645
nzmath.prime.FactoredInteger.__init__
def __init__(self, integer, factors=None)
Definition: prime.py:662
nzmath.prime.Zeta.size
size
Definition: prime.py:502
nzmath.prime.Status.__init__
def __init__(self)
Definition: prime.py:889
nzmath.prime.ImplementLimit
Definition: prime.py:1207
nzmath.prime.Zeta.__lshift__
def __lshift__(self, offset)
Definition: prime.py:540
nzmath.prime.JacobiSum
Definition: prime.py:1075
nzmath.prime.Zeta._pow_mod
def _pow_mod(self, e, mod)
Definition: prime.py:564
nzmath.prime.FactoredInteger.__ne__
def __ne__(self, other)
Definition: prime.py:729
nzmath.prime.aks
def aks(n)
Definition: prime.py:1215
nzmath.prime.FactoredInteger.__str__
def __str__(self)
Definition: prime.py:717
nzmath.prime.Status.done
def done(self, key)
Definition: prime.py:898
nzmath.prime.Status.subrest
def subrest(self, p, n, et, J, ub=200)
Definition: prime.py:1033
nzmath.prime._lucas_test_sequence
def _lucas_test_sequence(n, a, b)
Definition: prime.py:125
nzmath.prime.Status.yet
def yet(self, key)
Definition: prime.py:892
nzmath.prime.JacobiSum.get
def get(self, group, p, q=None)
Definition: prime.py:1082
nzmath.arith1
Definition: arith1.py:1
nzmath.poly.array
Definition: array.py:1
nzmath.plugins
Definition: plugins.py:1
nzmath.prime.TestPrime.__init__
def __init__(self, t=12)
Definition: prime.py:855
nzmath.prime.Status.yet_keys
def yet_keys(self)
Definition: prime.py:904
nzmath.prime.JacobiSum.makefx
def makefx(q)
Definition: prime.py:1126
nzmath.prime.FactoredInteger.__iter__
def __iter__(self)
Definition: prime.py:689
nzmath.prime.Zeta.__mod__
def __mod__(self, mod)
Definition: prime.py:607
nzmath.prime.Zeta.__mul__
def __mul__(self, other)
Definition: prime.py:518
nzmath.prime.FactoredInteger.__long__
def __long__(self)
Definition: prime.py:756
nzmath.factor.ecm.add
def add(P, Q, M, n)
Definition: ecm.py:314
nzmath.prime.FactoredInteger.__rfloordiv__
def __rfloordiv__(self, other)
Definition: prime.py:772
nzmath.prime.primeq
def primeq(n)
Definition: prime.py:380
nzmath.prime.Zeta.__setitem__
def __setitem__(self, position, value)
Definition: prime.py:618
nzmath.prime.Status.sub4
def sub4(self, q, n, J)
Definition: prime.py:1010
nzmath.prime.FactoredInteger.__ge__
def __ge__(self, other)
Definition: prime.py:750
nzmath.prime.TestPrime.primes
primes
Definition: prime.py:852
nzmath.prime.Zeta.__hash__
def __hash__(self)
Definition: prime.py:642
nzmath.bigrange
Definition: bigrange.py:1
nzmath.prime.TestPrime
Definition: prime.py:851
nzmath.prime.FactoredInteger.integer
integer
Definition: prime.py:670
nzmath.prime.Status.d
d
Definition: prime.py:890
nzmath.prime.generator
def generator()
Definition: prime.py:256
nzmath.prime.TestPrime.t
t
Definition: prime.py:857
nzmath.prime.Zeta.mass
def mass(self)
Definition: prime.py:651
nzmath.prime.bigprimeq
def bigprimeq(z)
Definition: prime.py:93
nzmath.prime._aks_mod
def _aks_mod(polynomial, r)
Definition: prime.py:1271
nzmath.prime.Zeta._square_mod
def _square_mod(self, mod)
Definition: prime.py:580
nzmath.prime.nextPrime
def nextPrime(n)
Definition: prime.py:328
nzmath.prime.primitive_root
def primitive_root(p)
Definition: prime.py:484
nzmath.prime.FactoredInteger.__pow__
def __pow__(self, other)
Definition: prime.py:707
nzmath.prime.Zeta.__eq__
def __eq__(self, other)
Definition: prime.py:636
nzmath.prime.JacobiSum.shelve
shelve
Definition: prime.py:1080
nzmath.prime.FactoredInteger.prime_divisors
def prime_divisors(self)
Definition: prime.py:844
nzmath.prime._isprime
def _isprime(n, pdivisors=None)
Definition: prime.py:413
nzmath.prime.TestPrime.next
def next(self)
Definition: prime.py:873
nzmath.prime.trialDivision
def trialDivision(n, bound=0)
Definition: prime.py:23
nzmath.prime.JacobiSum.__init__
def __init__(self)
Definition: prime.py:1079
nzmath.prime.prime
def prime(s)
Definition: prime.py:239
nzmath.prime.FactoredInteger.__eq__
def __eq__(self, other)
Definition: prime.py:720
nzmath.config
Definition: config.py:1