Programming Praxis

2 Responses to “Almost Primes”

1. Paul said

   April 22, 2019 at 7:29 AM

   In Python. almost_prime simply iterates over all numbers, factorises them and checks if they comply. This is slow for higher k and n. almost_prime3 builds numbers from factors and is a lot faster.

   from heapq import heappop, heappush
   from itertools import count
   from functools import partial, reduce
   from operator import mul
   
   from primesieve import n_primes  # primesieve from pypi.org
   
   from ma.primee import factors_auto  # td_factors and rho_factors
   
   def almost_prime(k, n, squarefree=False):
       """generate n k-almost_primes
          simply try all numbers and see if they check
       """
       res = []
       for i in count(2):
           facs = factors_auto(i)
           if len(facs) == k:
               if not squarefree or len(set(facs)) == k:
                   res.append(i)
                   if len(res) == n:
                       return res  
   
   # these are slow (especially skap)
   kap, skap = (partial(almost_prime, squarefree=b) for b in (False, True))
   
   def almost_prime3(k, n, squarefree=False):
       """generate n k-almost_primes
          use a priority queue (heap) with the value of the number as key
          pop the lowest number and increase the primes to get new candidate numbers
       """
       res = []
       prs = n_primes(k+n)
       next_prime = {k: v for k, v in zip(prs, prs[1:])}
       first = prs[:k] if squarefree else [2] * k
       heap = [(reduce(mul, first), first)]
       lastn = 0
       while heap:
           number, factors = heappop(heap)
           if number == lastn: # sometimes combinations are repeated
               continue
           lastn = number
           res.append(number)
           if len(res) == n:
               return res
           lastf = 0
           for i, f in enumerate(factors):
               if f == lastf:  # avoid as much as possible repetitions
                   continue
               lastf = f
               p = next_prime[f]
               if p and (not squarefree or p not in factors):
                   factors2 = list(factors)
                   factors2[i] = p
                   factors2.sort()
                   heappush(heap, (reduce(mul, factors2), factors2))
   
   kap3, skap3 = (partial(almost_prime3, squarefree=b) for b in (False, True))
   
   k, n = 20, 10
   print(kap3(k, n))
   print(skap3(k, n))
   
   [1048576, 1572864, 2359296, 2621440, 3538944, 3670016, 3932160, 5308416, 5505024, 5767168]
   [557940830126698960967415390, 573657473228859495079173570, 607905680585806330606288410, 620807402535341097414448110, 652240688739662165637964470, 657870531044913700245161430, 667699681954902035256087270, 676402095299700001660518090, 690333569478797019502056330, 691180431350985280004410110]
   
2. Ernie said

   April 25, 2019 at 3:25 PM

   Here is a faster algorithm for square free

   Select the smallest k-almost-prime (the product of the first k primes)
   Create z, the list of candidates from the smallest k-almost-prime. A candidate is an array
   of primes of length k
   Create a list, y, of next possible candidates from min of z.
   Example: say the current min for k = 3 is {0,3,4} where the values in min are indices
   into a list of primes. Increment each value. Then y is {1,3,4}, {0,4,5}, {0,3,5}. Note
   that {0,4,4} is not allowed because we require square free.
   Select the smallest of these candidates, add its product to result, remove it from y and
   then add the remaining y[i] to z
   Go to 3. Stop when result.Count = 10

   For k = 5 there is little difference between brute force and this algorithm but for k > 5 the difference in execution time
   becomes pronounced; for k = 7, brute force = 550 msec; above algorithm = 10 msec

   Still working on not square free