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

Posted by programmingpraxis
Filed in Exercises
2 Comments »

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

Leave a comment