Almost Primes

April 19, 2019

We determine if a number is k-almost-prime by factoring it. Any factoring algorithm can be used; here we use a 2,3,5-wheel:

(define (factors n)
  (let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)))
    (let loop ((n n) (f 2) (w 0) (fs (list)))
      (if (< n (* f f)) (reverse (cons n fs))
        (if (zero? (modulo n f))
            (loop (/ n f) f w (cons f fs))
            (loop n (+ f (vector-ref wheel w))
                  (if (= w 10) 3 (+ w 1)) fs))))))

Then a number is k-almost-prime if the length of the list of factors is k:

(define (kap? k n) (= (length (factors n)) k))

We can check if a k-almost-prime is squarefree by checking that all its factors are distinct:

(define (distinct? xs)
  (if (or (null? xs) (null? (cdr xs))) #t
    (if (= (car xs) (cadr xs)) #f
      (distinct? (cdr xs)))))
(define (sfkap? k n)
  (let ((fs (factors n)))
    (and (= (length fs) k) (distinct? fs))))

We can generate both sequences by passing the k-almost-prime test as a parameter:

(define (seq test? k len)
  (let loop ((k k) (m len) (n 2) (ks (list)) (kss (list)))
    (if (zero? k) kss
      (if (zero? m)
          (loop (- k 1) len 2 (list) (cons (reverse ks) kss))
          (if (test? k n)
              (loop k (- m 1) (+ n 1) (cons n ks) kss)
              (loop k m (+ n 1) ks kss))))))

And here is the output:

> (seq kap? 5 10)
((2 3 5 7 11 13 17 19 23 29)
  (4 6 9 10 14 15 21 22 25 26)
  (8 12 18 20 27 28 30 42 44 45)
  (16 24 36 40 54 56 60 81 84 88)
  (32 48 72 80 108 112 120 162 168 176))
> (seq sfkap? 5 10)
((2 3 5 7 11 13 17 19 23 29)
  (6 10 14 15 21 22 26 33 34 35)
  (30 42 66 70 78 102 105 110 114 130)
  (210 330 390 462 510 546 570 690 714 770)
  (2310 2730 3570 3990 4290 4830 5610 6006 6090 6270))

If you have an interest in number theory, you might want to search each of those sequences in OEIS, which will tell you much more about them, in fascinating detail.

You can run the program at https://ideone.com/lNRjAI.

Pages: 1 2

2 Responses to “Almost Primes”

  1. Paul said

    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

    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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: