Serg Iakovlev

Код
 	
 import math
 from datetime import datetime, timedelta
 
 PRIMES_LE_31 = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31)
 
 def gcd(a, b):
     """
     Return the greatest common divisor of 2 integers a and b.
     """
     a, b = abs(a), abs(b)
     while b:
         a, b = b, a % b
     return a
 
 def legendre(a, m):
     """
     This function returns the Legendre symbol (a/m).
     If m is an odd composite then this is the Jacobi symbol.
     """
     a = a % m
     symbol = 1
     while a != 0:
         while a & 1 == 0:
             a >>= 1
             if m & 7 == 3 or m & 7 == 5:
                 symbol = -symbol
         a, m = m, a
         if a & 3 == 3 and m & 3 == 3:
             symbol = -symbol
         a = a % m
     if m == 1:
         return symbol
     return 0
 
 def vp(n, p, k=0):
     """
     Return p-adic valuation and indivisible part of given integer.
 
     For example:
     >>> vp(100, 2)
     (2, 25)
 
     That means, 100 is 2 times divisible by 2, and the factor 25 of
     100 is indivisible by 2.
 
     The optional argument k will be added to the valuation.
     """
     q = p
     while not (n % q):
         k += 1
         q *= p
 
     return (k, n // (q // p))
 
 def log(n, base=2):
     """
     Return the integer part of logarithm of the given natural number
     'n' to the 'base'.  The default value for 'base' is 2.
     """
     if n < base:
         return 0
     if base == 10:
         return _log10(n)
     if base == 2 and n < POWERS_OF_2[-1]:
         # a shortcut
         return bisect.bisect_right(POWERS_OF_2, n) - 1
     fit = base
     result = 1
     stock = [(result, fit)]
     while fit < n:
         next = fit ** 2
         if next <= n:
             fit = next
             result += result
             stock.append((result, fit))
         else:
             break
     else: # just fit
         return result
     stock.reverse()
     for index, power in stock:
         prefit = fit * power
         if prefit == n:
             result += index
             break
         elif prefit < n:
             fit = prefit
             result += index
     return result
 
 def _log10(n):
     return len(str(n))-1
 
 
 def bigrange(start, stop=None, step=None):
     """
     Generate range like finite integer sequence but can generate more
     than sys.maxint elements.
     """
     if step is None:
         step = 1
     elif not isinstance(step, (int, long)):
         raise ValueError("non-integer step for range()")
     if not isinstance(start, (int, long)):
         raise ValueError("non-integer arg 1 for range()")
     if stop is None:
         start, stop = 0, start
     elif not isinstance(stop, (int, long)):
         raise ValueError("non-integer stop for range()")
 
     if step > 0:
         n = start
         while n < stop:
             yield n
             n += step
     elif step < 0:
         n = start
         while n > stop:
             yield n
             n += step
     else:
         raise ValueError("zero step for range()")
 
 def floorsqrt(a):
     """
     Return the floor of square root of the given integer.
     """
     if a < (1 << 59):
         return int(math.sqrt(a))
     else:
         # Newton method
         x = pow(10, (log(a, 10) >> 1) + 1) # compute initial value
         while True:
             x_new = (x + a//x) >> 1
             if x <= x_new:
                 return x
             x = x_new
 
 def product(iterable, init=None):
     """
     product(iterable) is a product of all elements in iterable.  If
     init is given, the multiplication starts with init instead of the
     first element in iterable. If the iterable is empty, then init or
     1 will be returned.
 
     If iterable is an iterator, it will be exhausted.
     """
     iterator = iter(iterable)
     if init is None:
         try:
             result = iterator.next()
         except StopIteration:
             return 1 # empty product
     else:
         result = init
     for item in iterator:
         result *= item
     return result
 
 PRIMONIAL_31 = product(PRIMES_LE_31)
 
 def spsp(n, base, s=None, t=None):
     """
     Strong Pseudo-Prime test.  Optional third and fourth argument
     s and t are the numbers such that n-1 = 2**s * t and t is odd.
     """
     if s is None or t is None:
         s, t = vp(n-1, 2)
     z = pow(base, t, n)
     if z != 1 and z != n-1:
         j = 0
         while j < s:
             j += 1
             z = pow(z, 2, n)
             if z == n-1:
                 break
         else:
             return False
     return True
 
 def full_euler(n, divisors):
     """
     Return True iff n is prime.
 
     The method proves the primality of n by the equality
     phi(n) = n - 1, where phi denotes the Euler totient.
     It requires a sequence of all prime divisors of n - 1.
 
     Brillhart,J. & Selfridge,J.L., Some Factorizations of $2^n\pm 1$
     and Related Results, Math.Comp., Vol.21, 87-96, 1967
     """
     m_order = n - 1
     primes = set(divisors)
     for g in bigrange(2, n):
         if pow(g, m_order, n) != 1:
             return False
         if 2 in primes:
             jacobi = legendre(g, n)
             if jacobi == 0:
                 return False
             elif jacobi == -1:
                 primes.remove(2)
         satisfied = [p for p in primes if p > 2 and pow(g, m_order // p, n) != 1]
         if satisfied:
             primes.difference_update(satisfied)
         if not primes:
             return True
     return True # only when n=2, flow reaches this line
 
 
 
 
 def smallSpsp(n, s=None, t=None):
     """
     4 spsp tests are sufficient to determine whether an integer less
     than 10**12 is prime or not.  Optional third and fourth argument
     s and t are the numbers such that n - 1 = 2**s * t and t is odd.
     """
     if s is None or t is None:
         s, t = vp(n - 1, 2)
     for p in (2, 13, 23, 1662803):
         if not spsp(n, p, s, t):
             return False
     return True
 
 
 
 def primonial(p):
     """
     Return 2*3*...*p for given prime p.
     """
     return product(generator_eratosthenes(p))
 
 
 # defs for APR algorithm
 
 def _isprime(n, pdivisors=None):
     """
     Return True iff n is prime.
 
     The check is done without APR.
     Assume that n is very small (less than 10**12) or
     prime factorization of n - 1 is known (prime divisors are passed to
     the optional argument pdivisors as a sequence).
     """
     if gcd(n, PRIMONIAL_31) > 1:
         return (n in PRIMES_LE_31)
     elif n < 10 ** 12:
         # 1369 == 37**2
         # 1662803 is the only prime base in smallSpsp which has not checked
         return n < 1369 or n == 1662803 or smallSpsp(n)
     else:
         return full_euler(n, pdivisors)
 
 def _factor(n, bound=0):
     """
     Trial division factorization for a natural number.
     Optional second argument is a search bound of primes.
     If the bound is given and less than the sqaure root of n,
     result is not proved to be a prime factorization.
     """
     factors = []
     if not (n & 1):
         v2, n = vp(n, 2)
         factors.append((2, v2))
     m = _calc_bound(n, bound)
     p = 3
     while p <= m:
         if not (n % p):
             v, n = vp(n, p)
             factors.append((p, v))
             m = _calc_bound(n, bound)
         p += 2
     if n > 1:
         factors.append((n, 1))
     return factors
 
 def _calc_bound(n, bound=0):
     if bound:
         m = min((bound, floorsqrt(n)))
     else:
         m = floorsqrt(n)
     return m
 
 def primitive_root(p):
     """
     Return a primitive root of p.
     """
     pd = FactoredInteger(p - 1).proper_divisors()
     for i in bigrange(2, p):
         for d in pd:
             if pow(i, (p - 1)//d, p) == 1:
                 break
         else:
             return i
 
 
 class Zeta(object):
     """
     Represent linear combinations of roots of unity.
     """
     def __init__(self, size, pos=None, val=1):
         self.size = size
         self.z = [0] * self.size
         if pos is not None:
             self.z[pos % self.size] = val
 
     def __add__(self, other):
         if self.size == other.size:
             m = self.size
             zr_a = Zeta(m)
             for i in range(m):
                 zr_a.z[i] = self.z[i] + other.z[i]
             return zr_a
         else:
             m = gcd.lcm(self.size, other.size)
             return self.promote(m) + other.promote(m)
 
     def __mul__(self, other):
         if not isinstance(other, Zeta):
             zr_m = Zeta(self.size)
             zr_m.z = [x*other for x in self.z]
             return zr_m
         elif self.size == other.size:
             zr_m = Zeta(self.size)
             other = +other
             for k in range(other.size):
                 if not other.z[k]:
                     continue
                 elif other.z[k] == 1:
                     zr_m = zr_m + (self << k)
                 else:
                     zr_m = zr_m + (self << k)*other.z[k]
             return zr_m
         else:
             m = gcd.lcm(self.size, other.size)
             return self.promote(m)*other.promote(m)
 
     __rmul__ = __mul__
 
     def __lshift__(self, offset):
         """
         The name is shift but the meaning of function is rotation.
         """
         new = Zeta(self.size)
         new.z = self.z[-offset:] + self.z[:-offset]
         return new
 
     def __pow__(self, e, mod=0):
         if mod:
             return self._pow_mod(e, mod)
 
         r = Zeta(self.size, 0)
         if e == 0:
             return r
         z = +self
         while True:
             if e & 1:
                 r = z*r
                 if e == 1:
                     return r
             e //= 2
             z = z._square()
 
     def _pow_mod(self, e, mod):
         r = Zeta(self.size, 0)
         if e == 0:
             return r
         z = self % mod
         while True:
             if e & 1:
                 r = z * r % mod
                 if e == 1:
                     return r
             e //= 2
             z = z._square_mod(mod)
 
     def _square(self):
         return self * self
 
     def _square_mod(self, mod):
         return self * self % mod
 
     def __pos__(self):
         m = self.size
         z_p = Zeta(m)
         if m & 1 == 0:
             mp = m >> 1
             for i in range(mp):
                 if self.z[i] > self.z[i+mp]:
                     z_p.z[i] = self.z[i] - self.z[i+mp]
                 else:
                     z_p.z[i+mp] = self.z[i+mp] - self.z[i]
         else:
             p = 3
             while m % p:
                 p += 2
             mp = m // p
             for i in range(mp):
                 mini = self.z[i]
                 for j in range(mp + i, m, mp):
                     if mini > self.z[j]:
                         mini = self.z[j]
                 for j in range(i, m, mp):
                     z_p.z[j] = self.z[j] - mini
         return z_p
 
     def __mod__(self, mod):
         """
         Return a new Zeta instance modulo 'mod'.
 
         All entries are reduced to the least absolute residues.
         """
         new = Zeta(self.size)
         half = (mod - 1) >> 1 # -1 to make 1 (mod 2) be 1, not -1.
         new.z = [(x + half) % mod - half for x in self.z]
         return new
 
     def __setitem__(self, position, value):
         self.z[position % self.size] = value
 
     def __getitem__(self, position):
         return self.z[position % self.size]
 
     def promote(self, size):
         if size == self.size:
             return +self
         new = Zeta(size)
         r = size // self.size
         for i in range(self.size):
             new.z[i*r] = self.z[i]
         return new
 
     def __len__(self):
         return self.size
 
     def __eq__(self, other):
         for i in range(self.size):
             if self.z[i] != other.z[i]:
                 return False
         return True
 
     def __hash__(self):
         return sum([hash(self.z[i]) for i in range(1, self.size)])
 
     def weight(self):
         """
         Return the number of nonzero entries.
         """
         return len(filter(None, self.z))
 
     def mass(self):
         """
         Return the sum of all entries.
         """
         return sum(self.z)
 
 
 class FactoredInteger(object):
     """
     Integers with factorization information.
     """
     def __init__(self, integer, factors=None):
         """
         FactoredInteger(integer [, factors])
 
         If factors is given, it is a dict of type {prime:exponent}
         and the product of prime**exponent is equal to the integer.
         Otherwise, factorization is carried out in initialization.
         """
         self.integer = int(integer)
         if factors is None:
             self.factors = dict(_factor(self.integer))
         else:
             self.factors = dict(factors)
 
     @classmethod
     def from_partial_factorization(cls, integer, partial):
         """
         Construct a new FactoredInteger object from partial
         factorization information given as dict of type
         {prime:exponent}.
         """
         partial_factor = 1
         for p, e in partial.iteritems():
             partial_factor *= p**e
         assert not integer % partial_factor, "wrong factorization"
         return cls(integer // partial_factor) * cls(partial_factor, partial)
 
     def __iter__(self):
         """
         Default iterator
         """
         return self.factors.iteritems()
 
     def __mul__(self, other):
         if isinstance(other, FactoredInteger):
             integer = self.integer * other.integer
             new_factors = self.factors.copy()
             for p in other.factors:
                 new_factors[p] = new_factors.get(p, 0) + other.factors[p]
             return self.__class__(integer, new_factors)
         else:
             return self * FactoredInteger(other)
 
     __rmul__ = __mul__
 
     def __pow__(self, other):
         new_integer = self.integer**other
         new_factors = {}
         for p in self.factors:
             new_factors[p] = self.factors[p] * other
         return self.__class__(new_integer, new_factors)
 
     def __pos__(self):
         return self.copy()
 
     def __str__(self):
         return str(self.integer)
 
     def __eq__(self, other):
         try:
             return self.integer == other.integer
         except AttributeError:
             return self.integer == int(other)
 
     def __hash__(self):
         return hash(self.integer)
 
     def __ne__(self, other):
         return not self.__eq__(other)
 
     def __le__(self, other):
         try:
             return self.integer <= other.integer
         except AttributeError:
             return self.integer <= int(other)
 
     def __lt__(self, other):
         try:
             return self.integer < other.integer
         except AttributeError:
             return self.integer < int(other)
 
     def __gt__(self, other):
         try:
             return self.integer > other.integer
         except AttributeError:
             return self.integer > int(other)
 
     def __ge__(self, other):
         try:
             return self.integer >= other.integer
         except AttributeError:
             return self.integer >= int(other)
 
     def __long__(self):
         return int(self.integer)
 
     __int__ = __long__
 
     def __mod__(self, other):
         # maybe you want self.is_divisible_by(other)
         try:
             if other.integer in self.factors:
                 return 0
             return self.integer % other.integer
         except AttributeError:
             if int(other) in self.factors:
                 return 0
             return self.integer % int(other)
 
     def __rfloordiv__(self, other):
         # assume other is an integer.
         return other // self.integer
 
     def copy(self):
         return self.__class__(self.integer, self.factors.copy())
 
     def is_divisible_by(self, other):
         """
         Return True if other divides self.
         """
         if int(other) in self.factors:
             # other is prime and divides
             return True
         return not self.integer % int(other)
 
     def exact_division(self, other):
         """
         Divide by a factor.
         """
         divisor = int(other)
         quotient = self.copy()
         if divisor in quotient.factors:
             if quotient.factors[divisor] == 1:
                 del quotient.factors[divisor]
             else:
                 quotient.factors[divisor] -= 1
         elif not isinstance(other, FactoredInteger):
             dividing = divisor
             for p, e in self.factors.iteritems():
                 while not dividing % p:
                     dividing //= p
                     if quotient.factors[p] == 1:
                         del quotient.factors[p]
                         assert dividing % p, dividing
                     else:
                         quotient.factors[p] -= 1
                 if dividing == 1:
                     break
             assert dividing == 1
         else:
             for p, e in other.factors.iteritems():
                 assert p in quotient.factors and quotient.factors[p] >= e
                 if quotient.factors[p] == e:
                     del quotient.factors[p]
                 else:
                     quotient.factors[p] -= e
         quotient.integer //= divisor
         return quotient
 
     # maybe this is what you want, isn't it?
     __floordiv__ = exact_division
 
     def divisors(self):
         """
         Return all divisors.
         """
         divs = [FactoredInteger(1)]
         for p, e in self.factors.iteritems():
             q = FactoredInteger(1)
             pcoprimes = list(divs)
             for j in range(1, e + 1):
                 q *= FactoredInteger(p, {p:1})
                 divs += [k * q for k in pcoprimes]
         return divs
 
     def proper_divisors(self):
         """
         Return the proper divisors (divisors of n excluding 1 and n).
         """
         return self.divisors()[1:-1]
 
     def prime_divisors(self):
         """
         Return the list of primes that divides the number.
         """
         return self.factors.keys()
 
 
 class TestPrime(object):
     primes = PRIMES_LE_31
     primecache = set(primes)
 
     def __init__(self, t=12):
         if isinstance(t, (int, long)):
             self.t = FactoredInteger(t)
         else:
             assert isinstance(t, FactoredInteger)
             self.t = t
         powerof2 = self.t.factors[2] + 2
         self.et = FactoredInteger(2 ** powerof2, {2:powerof2})
         for d in self.t.divisors():
             # d is an instance of FactoredInteger
             p = d.integer + 1
             if p & 1 and (p in self.primecache or _isprime(p, d.factors)):
                 self.et = self.et * FactoredInteger(p, {p:1})
                 if p in self.t.factors:
                     e = self.t.factors[p]
                     self.et = self.et * FactoredInteger(p**e, {p:e})
                 self.primecache.add(p)
 
     def next(self):
         eu = []
         for p in self.primes:
             if p in self.t.factors:
                 eu.append((p - 1) * p**(self.t.factors[p] - 1))
             else:
                 eu.append(p - 1)
                 break
         p = self.primes[eu.index(min(eu))]
         return self.__class__(self.t * FactoredInteger(p, {p:1}))
 
 
 class Status:
     """
     status collector for apr.
     """
     def __init__(self):
         self.d = {}
 
     def yet(self, key):
         """
         set key's status be 'yet'.
         """
         self.d[key] = 0
 
     def done(self, key):
         """
         set key's status be 'done'.
         """
         self.d[key] = 1
 
     def yet_keys(self):
         """
         Return keys whose stati are 'yet'.
         """
         return [k for k in self.d if not self.d[k]]
 
     def isDone(self, key):
         """
         Return whether key's status is 'done' or not.
         """
         return self.d[key]
 
     def subodd(self, p, q, n, J):
         """
         Return the sub result for odd key 'p'.
         If it is True, the status of 'p' is flipped to 'done'.
         """
         s = J.get(1, p, q)
         Jpq = J.get(1, p, q)
         m = s.size
         for x in range(2, m):
             if x % p == 0:
                 continue
             sx = Zeta(m)
             i = x
             j = 1
             while i > 0:
                 sx[j] = Jpq[i]
                 i = (i + x) % m
                 j += 1
             sx[0] = Jpq[0]
             sx = pow(sx, x, n)
             s = s*sx % n
         s = pow(s, n//m, n)
         r = n % m
         t = 1
         for x in range(1, m):
             if x % p == 0:
                 continue
             c = (r*x) // m
             if c:
                 tx = Zeta(m)
                 i = x
                 j = 1
                 while i > 0:
                     tx[j] = Jpq[i]
                     i = (i + x) % m
                     j += 1
                 tx[0] = Jpq[0]
                 tx = pow(tx, c, n)
                 t = t*tx % n
         s = +(t*s % n)
         if s.weight() == 1 and s.mass() == 1:
             for i in range(1, m):
                 if gcd(m, s.z.index(1)) == 1:
                     self.done(p)
                 return True
         return False
 
     def sub8(self, q, k, n, J):
         s = J.get(3, q)
         J3 = J.get(3, q)
         m = len(s)
         sx_z = {1:s}
         x = 3
         step = 2
         while m > x:
             z_4b = Zeta(m)
             i = x
             j = 1
             while i != 0:
                 z_4b[j] = J3[i]
                 i = (i + x) % m
                 j += 1
             z_4b[0] = J3[0]
             sx_z[x] = z_4b
             s = pow(sx_z[x], x, n) * s
             step = 8 - step
             x += step
 
         s = pow(s, n//m, n)
 
         r = n % m
         step = 2
         x = 3
         while m > x:
             c = r*x
             if c > m:
                 s = pow(sx_z[x], c//m, n) * s
             step = 8 - step
             x += step
         r = r & 7
         if r == 5 or r == 7:
             s = J.get(2, q).promote(m) * s
         s = +(s % n)
 
         if s.weight() == 1 and s.mass() == 1:
             if gcd(m, s.z.index(1)) == 1 and pow(q, (n-1) >> 1, n) == n-1:
                 self.done(2)
             return True
         elif s.weight() == 1 and s.mass() == n-1:
             if gcd(m, s.z.index(n-1)) == 1 and pow(q, (n-1) >> 1, n) == n-1:
                 self.done(2)
             return True
         return False
 
     def sub4(self, q, n, J):
         j2 = J.get(1, 2, q)**2
         s = q*j2 % n
         s = pow(s, n >> 2, n)
         if n & 3 == 3:
             s = s*j2 % n
         s = +(s % n)
         if s.weight() == 1 and s.mass() == 1:
             i = s.z.index(1)
             if (i == 1 or i == 3) and pow(q, (n-1) >> 1, n) == n-1:
                 self.done(2)
             return True
         return False
 
     def sub2(self, q, n):
         s = pow(n-q, (n-1) >> 1, n)
         if s == n-1:
             if n & 3 == 1:
                 self.done(2)
         elif s != 1:
             return False
         return True
 
     def subrest(self, p, n, et, J, ub=200):
         if p == 2:
             q = 5
             while q < 2*ub + 5:
                 q += 2
                 if not _isprime(q) or et % q == 0:
                     continue
                 if n % q == 0:
                     _log.info("%s divides %s.\n" % (q, n))
                     return False
                 k = vp(q-1, 2)[0]
                 if k == 1:
                     if n & 3 == 1 and not self.sub2(q, n):
                         return False
                 elif k == 2:
                     if not self.sub4(q, n, J):
                         return False
                 else:
                     if not self.sub8(q, k, n, J):
                         return False
                 if self.isDone(p):
                     return True
             else:
                 raise ImplementLimit("limit")
         else:
             step = p*2
             q = 1
             while q < step*ub + 1:
                 q += step
                 if not _isprime(q) or et % q == 0:
                     continue
                 if n % q == 0:
                     _log.info("%s divides %s.\n" % (q, n))
                     return False
                 if not self.subodd(p, q, n, J):
                     return False
                 if self.isDone(p):
                     return True
             else:
                 raise ImplementLimit("limit")
 
 
 class JacobiSum:
     """
     A repository of the Jacobi sums.
     """
     def __init__(self):
         self.shelve = {}
 
     def get(self, group, p, q=None):
         if q:
             assert group == 1
             if (group, p, q) not in self.shelve:
                 self.make(q)
             return self.shelve[group, p, q]
         else:
             assert group == 2 or group == 3
             if (group, p) not in self.shelve:
                 self.make(p)
             return self.shelve[group, p]
 
     def make(self, q):
         fx = self.makefx(q)
         qpred = q - 1
         qt = _factor(qpred)
         qt2 = [k for (p, k) in qt if p == 2][0]
         k, pk = qt2, 2**qt2
         if k >= 2:
             J2q = Zeta(pk, 1 + fx[1])
             for j in range(2, qpred):
                 J2q[j + fx[j]] = J2q[j + fx[j]] + 1
             self.shelve[1, 2, q] = +J2q
             if k >= 3:
                 J2 = Zeta(8, 3 + fx[1])
                 J3 = Zeta(pk, 2 + fx[1])
                 for j in range(2, qpred):
                     J2[j*3 + fx[j]] = J2[j*3 + fx[j]] + 1
                     J3[j*2 + fx[j]] = J3[j*2 + fx[j]] + 1
                 self.shelve[3, q] = +(self.shelve[1, 2, q]*J3)
                 self.shelve[2, q] = +(J2*J2)
         else:
             self.shelve[1, 2, q] = 1
         for (p, k) in qt:
             if p == 2:
                 continue
             pk = p**k
             Jpq = Zeta(pk, 1 + fx[1])
             for j in range(2, qpred):
                 Jpq[j + fx[j]] = Jpq[j + fx[j]] + 1
             self.shelve[1, p, q] = +Jpq
         del fx
 
     @staticmethod
     def makefx(q):
         """
         Return a dict called 'fx'.
         The dict stores the information that fx[i] == j iff
           g**i + g**j = 1 mod q
         for g, a primitive root of the prime q.
         """
         g = primitive_root(q)
         qpred = q - 1
         qd2 = qpred >> 1
         g_mf = [0, g]
         for _ in range(2, qpred):
             g_mf.append((g_mf[-1]*g) % q)
         fx = {}
         for i in range(1, qd2):
             if i in fx:
                 continue
             # search j s.t. g**j + g**i = 1 mod q
             j = g_mf.index(q + 1 - g_mf[i])
             fx[i] = j
             fx[j] = i
             fx[qpred - i] = (j - i + qd2) % qpred
             fx[fx[qpred - i]] = qpred - i
             fx[qpred - j] = (i - j + qd2) % qpred
             fx[fx[qpred - j]] = qpred - j
         del g_mf
         return fx
 
 
 def apr(n):
     """
     apr is the main function for Adleman-Pomerance-Rumery primality test.
     Assuming n has no prime factors less than 32.
     Assuming n is spsp for several bases.
     """
     L = Status()
 
     rb = floorsqrt(n) + 1
     el = TestPrime()
     while el.et <= rb:
         el = el.next()
 
     #print 't=%s e(t)=%s' % (el.t, el.et)
 
     plist = el.t.factors.keys()
     plist.remove(2)
     #print plist
     L.yet(2)
     for p in plist:
         if pow(n, p-1, p*p) != 1:
             L.done(p)
         else:
             L.yet(p)
     qlist = el.et.factors.keys()
     qlist.remove(2)
     #print qlist
     J = JacobiSum()
     for q in qlist:
         for p in plist:
             if (q-1) % p != 0:
                 continue
             if not L.subodd(p, q, n, J):
 		print 'subodd'
                 return False
         k = vp(q-1, 2)[0]
         if k == 1:
             if not L.sub2(q, n):
 		#print 'sub2'
                 return False
         elif k == 2:
             if not L.sub4(q, n, J):
 		print 'sub4'
                 return False
         else:
             if not L.sub8(q, k, n, J):
 		print 'sub8'
                 return False
     for p in L.yet_keys():
         if not L.subrest(p, n, el.et, J):
 	    print 'subrest'
             return False
     r = int(n)
     for _ in bigrange(1, el.t.integer):
         r = (r*n) % el.et.integer
         if n % r == 0 and r != 1 and r != n:
             _log.info("%s divides %s.\n" %(r, n))
             print "%s divides %s.\n" %(r, n)
             return False
     return True
 
 n=4806682957193450524071628095151496599650578449844083078134025220629010456661487588579545817257185675376893
 aaa = datetime.now()
 print apr(n)
 bbb = datetime.now()
 print 'diff %s' % (bbb-aaa).total_seconds()
Тест Адлемана — Померанса — Румели

Символ Лежандра

Символ Дирихле

Символ Якоби

Первообразный корень

Китайская теорема об остатках.

Реализация APR
