Factorise

Simple Python module for prime factorisation. Uses the 6n plus/minus 1 phenomenom. Reasonably fast.

Last modified
Lines 162
Indexable Yes

Parent directory Download CGIread sitemap Main page

Quick links: factorise main make_offsets

  1. #!/usr/bin/python
  2. def make_offsets(primes):
  3.     '''
  4.     `primes` MUST be a list of all primes from 2 to stop.
  5.     
  6.     All primes except 2 and 3 are 6n - 1 or 6n + 1 because
  7.     2*3 = 6 (the period) and 6n+1 and 6n-1 are the only possible
  8.     primes when n > 0.
  9.     
  10.     This function extends this beyond just 2, 3, 6n+1 and 6n+5.
  11.         
  12.     Returns: period, offsets
  13.     
  14.     Notice: `offsets` will include 1 and will not include any prime
  15.     in `primes`.  Exceptions need to be made when n == 0.
  16.     
  17.     Example
  18.     =======
  19.     
  20.         >>> period, offsets = table(2, 3, 5)
  21.         >>> period
  22.         30
  23.         >>> offsets
  24.         [1, 7, 11, 13, 17, 19, 23, 29]
  25.     '''
  26.     # Calculate the period.
  27.     period = 1
  28.     for p in primes:
  29.         period *= p
  30.     # Calculate valid offsets
  31.     offsets = []
  32.     i = 0
  33.     while i < period:
  34.         if all(map(lambda x: i%x, primes)):
  35.             offsets.append(i)
  36.         i += 1
  37.     # Done.
  38.     return period, offsets
  39. def factorise(num):
  40.     '''
  41.     Return a list of (prime, power).
  42.     '''
  43.     # Handle bad values of `num`.    
  44.     factors = []
  45.     if num != int(num):
  46.         raise TypeError('`num` MUST be an integer.')    
  47.     elif num < -1:
  48.         factors.append((-1, 1))
  49.         num = -num
  50.     elif num == -1:
  51.         return [(-1, 1)]
  52.     elif num == 0:
  53.         return [(0, 1)]
  54.     elif num == 1:
  55.         return [(1, 1)]
  56.     
  57.     # There is no need to only use primes here.
  58.     # Consider this example: 2, 3, 4, 5
  59.     # There will be no 4s because the 2s have already been
  60.     # divided away.
  61.     # 
  62.     # Checking whether or not the number is a prime will take
  63.     # a longer time than realising it won't divide `num`, so
  64.     # why bother?
  65.     # 
  66.     # The "only" thing that could make this faster is to not use
  67.     # so damn many composite numbers in the first place, that's
  68.     # where the `make_offsets` function comes to use.
  69.     #
  70.     # Use the `make_offsets` function to exclude some composite numbers.
  71.     # Primes are {2, 3, 5} when n == 0 and {6n+1, 6n+5} when n > 0
  72.     # Extend this beyond just 2*3, go to 2*3*5*7*11*13*17
  73.     # 
  74.     #global primes      # Won't be modified, global declarations only
  75.     #global period      # here for clarity.
  76.     #global offsets
  77.     
  78.     table = primes + offsets[1:]    # Offset table to use when n == 0.
  79.     
  80.     n_times_period = 0      # n alone is only needed for checking if n is
  81.                             # zero or not.
  82.     
  83.     sqrt_num = num**0.5     # If the potential divident > sqrt(num) then
  84.                             # there must have been a previous divident, but
  85.                             # there wasn't. Therefore`, `num` must be prime.
  86.     
  87.     # Can't handle num == 1 properly.
  88.     while num > 1:
  89.         for offset in table:
  90.             divident = n_times_period + offset
  91.             # Check if `num` is prime.
  92.             if divident > sqrt_num:
  93.                 factors.append((num, 1))
  94.                 num = 1
  95.                 break
  96.             # Divide away the selected factor.
  97.             power = 0
  98.             while num % divident == 0:
  99.                 power += 1
  100.                 num /= divident
  101.             if power:
  102.                 factors.append((divident, power))
  103.                 # `num` does not magically become 1 -> continue in same block.
  104.                 if num == 1:
  105.                     break
  106.                 # `num` has changed, change `sqrt_num`
  107.                 sqrt_num = num**0.5
  108.         # Offset table to use when n > 0.
  109.         # Change only once.
  110.         if n_times_period == 0:
  111.             table = offsets
  112.         # Next iteration:
  113.         n_times_period += period
  114.     
  115.     return factors
  116. def main():
  117.     import sys
  118.     try:
  119.         while True:
  120.             # Get input
  121.             sys.stderr.write('\nEnter number to factorise: ')
  122.             sys.stderr.flush()
  123.             line = sys.stdin.readline()
  124.             if not line:
  125.                 break
  126.             if line == '\n':
  127.                 # Shut up.
  128.                 continue
  129.             # Parse input
  130.             try:
  131.                 num = int(line)
  132.             except ValueError as e:
  133.                 sys.stderr.write(str(e) + '\n')
  134.                 continue
  135.             # Compute input
  136.             try:
  137.                 factors = factorise(num)
  138.             except ValueError as e:
  139.                 sys.stderr.write(str(e) + '\n')
  140.                 continue
  141.             # Print out.
  142.             s = ' * '.join(map(lambda x: '{}**{}'.format(*x), factors))
  143.             sys.stdout.write('{} = {}\n'.format(num, s))
  144.     except KeyboardInterrupt:
  145.         pass
  146.     sys.stderr.write('\n')
  147. # Pre-initialize the offsets table.
  148. primes = [2, 3, 5, 7, 11, 13, 17]   # Becomes slow at 19.
  149. period, offsets = make_offsets(primes)
  150. if __name__ == '__main__':
  151.     main()