Need a unique gift idea?
A Pastebin account makes a great Christmas gift
SHARE
TWEET

Untitled

a guest Jul 15th, 2018 67 Never
Upgrade to PRO!
ENDING IN00days00hours00mins00secs
 
  1. """Prime number generation and number factorization."""
  2.  
  3. import bisect, itertools, random, sys
  4.  
  5. _primes = [2]
  6.  
  7. _miller_rabin_limit = 48611    # 5000th prime
  8. _miller_rabin_security = 7
  9.  
  10. def modpow (a, b, c):
  11.     """Efficiently compute (a^b)%c where a, b and c are positive integers."""
  12.     if b == 1: return a % c
  13.     d = b//2
  14.     x = modpow (a, d, c)
  15.     x = x*x%c
  16.     if b % 2 == 1: x = x*a%c
  17.     return x
  18.  
  19. def miller_rabin (n, t = _miller_rabin_security):
  20.     """Apply Miller-Rabin primality test to detect whether n is prime or not.
  21.     The test is run t times."""
  22.     if n % 2 == 0: return False
  23.     r = (n-1)//2
  24.     for s in itertools.count (1):
  25.         if r % 2: break
  26.         r //= 2
  27.     for tt in xrange (t):
  28.         a = int (random.uniform (2, n-1))
  29.         y = modpow (a, r, n)
  30.         if y != 1 and y != n-1:
  31.             for j in xrange (1, s):
  32.                 y = y**2 % n
  33.                 if y == 1: return False
  34.                 if y == n-1: break
  35.             if y != n-1: return False
  36.     return True
  37.  
  38. def _is_prime (i):
  39.     if i > _miller_rabin_limit: return miller_rabin (i)
  40.     s = int (i**0.5)
  41.     for j in _primes:
  42.         if i % j == 0: return False
  43.         if j > s: return True
  44.     return True
  45.  
  46. def _gen_primes ():
  47.     for i in xrange (3, sys.maxint, 2):
  48.         if _is_prime (i):
  49.             _primes.append (i)
  50.             yield i
  51.  
  52. _primary = _gen_primes ()
  53.  
  54. def _gen_primes (minval):
  55.     i = 0
  56.     if minval:
  57.         l = _primes[-1]
  58.         if minval > l:
  59.             while minval > l: l =_primary.next ()
  60.             i = len (_primes) - 1
  61.         else:
  62.             i = bisect.bisect_left (_primes, minval)
  63.     for i in itertools.count (i):
  64.         if i == len (_primes): yield _primary.next ()
  65.         else: yield _primes[i]
  66.  
  67. def gen_primes (maxval = None, minval = None):
  68.     """Prime numbers generator. If minval is given, only primes greater or
  69.     equal to minval are returned. If maxval is given, only primes smaller
  70.     or equal to maxval are returned.
  71.  
  72.     >>> for p in gen_primes(5): print p
  73.     ...
  74.     2
  75.     3
  76.     5
  77.     7
  78.     11"""
  79.     if maxval:
  80.         return itertools.takewhile (lambda x: x<=maxval, _gen_primes (minval))
  81.     return _gen_primes (minval)
  82.  
  83. def gen_factors (n, duplicates = True):
  84.     """Generator for factors of n (n > 1). If duplicates is False, do not
  85.     send the same factor more than once."""
  86.     assert n > 1
  87.     if n > _miller_rabin_limit and miller_rabin (n):
  88.         yield n
  89.         return
  90.     s = int (n**0.5)
  91.     for i in gen_primes ():
  92.         if i > s:
  93.             yield n
  94.             return
  95.         if n % i == 0:
  96.             yield i
  97.             n //= i
  98.             while n % i == 0:
  99.                 if duplicates: yield i
  100.                 n //= i
  101.             if n == 1: return
  102.             if n > _miller_rabin_limit and miller_rabin (n):
  103.                 yield n
  104.                 return
  105.             s = int (n**0.5)
  106.  
  107. def factors (n):
  108.     """Factorize n (n > 1) into its prime factors. Return a dictionary where
  109.     keys are prime factors and values are powers.
  110.  
  111.     >>> factors(18)
  112.     {2: 1, 3: 2}"""
  113.     l = {}
  114.     for i in gen_factors (n):
  115.         try: l[i] += 1
  116.         except KeyError: l[i] = 1
  117.     return l
  118.  
  119. def factorslist (n):
  120.     """Return a list of prime factors of n (n > 1).
  121.  
  122.     >>> factorslist(18)
  123.     [2, 3, 3]"""
  124.     return list (gen_factors (n))
  125.  
  126. def is_prime (n):
  127.     """Check whether n is prime or not."""
  128.     if n < 2: return False
  129.     return gen_factors(n).next() == n
  130.  
  131. def ufactors (n):
  132.     """Return the list of unique prime factors of n (n > 1).
  133.  
  134.     >>> ufactors(100)
  135.     [2, 5]"""
  136.     return list (gen_factors(n, duplicates = False))
  137.  
  138. def nfactors (n):
  139.     """Return the number of unique prime factors of n."""
  140.     return len (ufactors (n))
  141.  
  142. def totient (n):
  143.     """Euler's totient function. Returns the number of integers between
  144.     1 and n-1 relatively prime to n (n > 0).
  145.  
  146.     >>> totient(6)
  147.     2
  148.     (6 is relatively prime to 1 and 5)"""
  149.     if n == 1: return 0
  150.     num = den = 1
  151.     for p in gen_factors (n, duplicates = False):
  152.         num *= (p-1)
  153.         den *= p
  154.     return n * num // den
  155.  
  156. _s = {1: 1}
  157. def s (n):
  158.     """Sum of divisors of n (n > 0).
  159.  
  160.     >>> s(6)
  161.     12
  162.     (divisors of 6 are 1, 2, 3 and 6, summing to 12)"""
  163.     if not _s.has_key (n):
  164.         t = 1
  165.         for (p, k) in factors(n).items():
  166.             t *= (p**(k+1)-1)//(p-1)
  167.         _s[n] = t
  168.     return _s[n]
  169.  
  170. if __name__ == '__main__':
  171.     # Quick test -- add primes up to 100,000 and compare to J result to:
  172.     #   +/p:i.(p:^:_1)100000x
  173.     assert sum (gen_primes (100000)) == 454396537
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top