fork(6) download
copy 
  1. # programming with prime numbers
  2.  
  3. from fractions import gcd
  4.  
  5. def primes(n):
  6.     if type(n) != int and type(n) != long:
  7.         raise TypeError('must be integer')
  8.     if n < 2:
  9.         raise ValueError('must be greater than one')
  10.     m = (n-1) // 2
  11.     b = [True] * m
  12.     i, p, ps = 0, 3, [2]
  13.     while p*p < n:
  14.         if b[i]:
  15.             ps.append(p)
  16.             j = 2*i*i + 6*i + 3
  17.             while j < m:
  18.                 b[j] = False
  19.                 j = j + 2*i + 3
  20.         i += 1; p += 2
  21.     while i < m:
  22.         if b[i]:
  23.             ps.append(p)
  24.         i += 1; p += 2
  25.     return ps
  26.  
  27. def td_prime(n, limit=1000000):
  28.     if type(n) != int and type(n) != long:
  29.         raise TypeError('must be integer')
  30.     if n % 2 == 0:
  31.         return n == 2
  32.     d = 3
  33.     while d * d <= n:
  34.         if limit < d:
  35.             raise OverflowError('limit exceeded')
  36.         if n % d == 0:
  37.             return False
  38.         d += 2
  39.     return True
  40.  
  41. def td_factors(n, limit=1000000):
  42.     if type(n) != int and type(n) != long:
  43.         raise TypeError('must be integer')
  44.     fs = []
  45.     while n % 2 == 0:
  46.         fs += [2]
  47.         n /= 2
  48.     if n == 1:
  49.         return fs
  50.     f = 3
  51.     while f * f <= n:
  52.         if limit < f:
  53.             raise OverflowError('limit exceeded')
  54.         if n % f == 0:
  55.             fs += [f]
  56.             n /= f
  57.         else:
  58.             f += 2
  59.     return fs + [n]
  60.  
  61. def is_prime(n):
  62.     if type(n) != int and type(n) != long:
  63.         raise TypeError('must be integer')
  64.     if n < 2:
  65.         return False
  66.     ps = [2,3,5,7,11,13,17,19,23,29,31,37,41,
  67.          43,47,53,59,61,67,71,73,79,83,89,97]
  68.     def is_spsp(n, a):
  69.         d, s = n-1, 0
  70.         while d%2 == 0:
  71.             d /= 2; s += 1
  72.         t = pow(a,d,n)
  73.         if t == 1:
  74.             return True
  75.         while s > 0:
  76.             if t == n-1:
  77.                 return True
  78.             t = (t*t) % n
  79.             s -= 1
  80.         return False
  81.     if n in ps: return True
  82.     for p in ps:
  83.         if not is_spsp(n,p):
  84.             return False
  85.     return True
  86.  
  87. def rho_factors(n, limit=1000000):
  88.     if type(n) != int and type(n) != long:
  89.         raise TypeError('must be integer')
  90.     def gcd(a,b):
  91.         while b: a, b = b, a%b
  92.         return abs(a)
  93.     def rho_factor(n, c, limit):
  94.         f = lambda(x): (x*x+c) % n
  95.         t, h, d = 2, 2, 1
  96.         while d == 1:
  97.             if limit == 0:
  98.                 raise OverflowError('limit exceeded')
  99.             t = f(t); h = f(f(h)); d = gcd(t-h, n)
  100.         if d == n:
  101.             return rho_factor(n, c+1, limit)
  102.         if is_prime(d):
  103.             return d
  104.         return rho_factor(d, c+1, limit)
  105.     if -1 <= n <= 1: return [n]
  106.     if n < -1: return [-1] + rho_factors(-n, limit)
  107.     fs = []
  108.     while n % 2 == 0:
  109.         n = n // 2; fs = fs + [2]
  110.     if n == 1: return fs
  111.     while not is_prime(n):
  112.         f = rho_factor(n, 1, limit)
  113.         n = n / f
  114.         fs = fs + [f]
  115.     return sorted(fs + [n])
  116.  
  117. print primes(100)
  118. print len(primes(1000000))
  119. print td_prime(600851475143)
  120. print td_factors(600851475143)
  121. print is_prime(600851475143)
  122. print is_prime(2305843009213693951)
  123. print rho_factors(600851475143)
Success #stdin #stdout 0.43s 5884KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
78498
False
[71, 839, 1471, 6857L]
False
True
[71L, 839L, 1471L, 6857L]
https://ideone.com/OBOZQ
language:
Python (cpython 2.7.16)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!