Positive number which has no prime factor greater than k

Question

I tried to create a function which takes two variables n and k.

The function returns the number of positive integers that have prime factors all less than or equal to k. The number of positive integers is limited by n which is the largest positive integer.

For example, if k = 4 and n = 10; the positive integers which have all prime factors less than or equal to 4 are 1, 2, 3, 4, 6, 8, 9, 12...(1 is always part for some reason even though its not prime) but since n is 10, 12 and higher numbers are ignored.

So the function will return 7. The code I wrote works for smaller values of n while it just keeps on running for larger values.

How can I optimize this code? Should I start from scratch and come up with a better algorithm?

int generalisedHammingNumbers(int n, int k)
{
    vector<int>store;
    vector<int>each_prime = {};

    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j <= i; ++j)
        {
            if (i%j == 0 && is_prime(j))
            {
                each_prime.push_back(j);   //temporary vector of prime factors for each integer(i) 
            }
        }
        for (int m = 0; m<each_prime.size(); ++m)
        {
            while(each_prime[m] <= k && m<each_prime.size()-1)  //search for prime factor greater than k
            {
                ++m;
            }
            if (each_prime[m] > k);  //do nothing for prime factor greater than k
            else store.push_back(i);  //if no prime factor greater than k, i is valid, store i
        }
        each_prime = {};
    }

    return (store.size()+1);
}

bool is_prime(int x)
{
    vector<int>test;
    if (x != 1)
    {
        for (int i = 2; i < x; ++i)
        {
            if (x%i == 0)test.push_back(i);
        }
        if (test.size() == 0)return true;
        else return false;

    }
    return false;

}


int main() 
{
    long n;
    int k;
    cin >> n >> k;
    long result = generalisedHammingNumbers(n, k);
    cout << result << endl;


}

Answer 1

Should I start from scratch and come up with a better algorithm?

Yes... I think so.

This seems to me a work for the Sieve of Eratosthenes.

So I propose to

1) create a std::vector<bool> to detect, through Eratosthenes, the primes to n

2) remove primes starting from k+1, and their multiples, from the pool of your numbers (another std::vector<bool>)

3) count the true remained values in the pool vector

The following is a full working example

#include <vector>
#include <iostream>
#include <algorithm>

std::size_t foo (std::size_t n, std::size_t k)
 {
   std::vector<bool> primes(n+1U, true);
   std::vector<bool> pool(n+1U, true);

   std::size_t const sqrtOfN = std::sqrt(n);

   // first remove the not primes from primes list (Sieve of Eratosthenes)
   for ( auto i = 2U ; i <= sqrtOfN ; ++i )
      if ( primes[i] )
         for ( auto j = i << 1 ; j <= n ; j += i )
            primes[j] = false;

   // then remove from pool primes, bigger than k, and multiples
   for ( auto i = k+1U ; i <= n ; ++i )
      if ( primes[i] )
         for ( auto j = i ; j <= n ; j += i )
            pool[j] = false;

   // last count the true value in pool (excluding the zero)
   return std::count(pool.begin()+1U, pool.end(), true);
 }

int main ()
 {
   std::cout << foo(10U, 4U) << std::endl;
 }

Answer 2

Here's an Eulerian version in Python (seems about 1.5 times faster than Paul Hankin's). We generate only the numbers themselves by multiplying a list by each prime and its powers in turn.

import time

start = time.time()

n = 1000000
k = 100

total = 1
a = [None for i in range(0, n+1)]
s = []
p = 1
while (p < k):
  p = p + 1

  if a[p] is None:
    #print("\n\nPrime: " + str(p))
    a[p] = True
    total = total + 1
    s.append(p)

    limit = n / p
    new_s = []

    for i in s:
      j = i
      while j <= limit:
        new_s.append(j)
        #print j*p
        a[j * p] = True
        total = total + 1
        j = j * p
    s = new_s

print("\n\nGilad's answer: " + str(total))
end = time.time()
print(end - start)

# Paul Hankin's solution
def limited_prime_factors(n, k):
    ps = [False] * (k+1)
    r = [True] * 2 + [False] * n
    for p in xrange(2, k+1):
        if ps[p]: continue
        for i in xrange(p, k+1, p):
            ps[i] = True
        for i in xrange(p, n+1, p):
            r[i] = r[i//p]
    return len([i for i, b in enumerate(r) if b]) - 1

start = time.time()
print "\nPaul's answer:" + str(limited_prime_factors(1000000, 100))
end = time.time()
print(end - start)

Answer 3

Generate the primes using a sieve of Erastothenes, and then use a modified coin-change algorithm to find numbers which are products of only those primes. In fact, one can do both simultaneously like this (in Python, but is easily convertible to C++):

def limited_prime_factors(n, k):
    ps = [False] * (k+1)
    r = [True] * 2 + [False] * n
    for p in xrange(2, k+1):
        if ps[p]: continue
        for i in xrange(p, k+1, p):
            ps[i] = True
        for i in xrange(p, n+1, p):
            r[i] = r[i//p]
    return [i for i, b in enumerate(r) if b]

print limited_prime_factors(100, 3)

The output is:

[0, 1, 2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 27, 32, 36, 48, 54, 64, 72, 81, 96]

Here, each time we find a prime p, we strike out all multiples of p in the ps array (as a standard Sieve of Erastothenes), and then in the r array, mark all multiples of any number that's a multiple of p whether their prime factors are all less than or equal to p.

It runs in O(n) space and O(n log log k) time, assuming n>k.

A simpler O(n log k) solution tests if all the factors of a number are less than or equal to k:

def limited_prime_factors(n, k):
    r = [True] * 2 + [False] * n
    for p in xrange(2, k+1):
        for i in xrange(p, n+1, p):
            r[i] = r[i//p]
    return [i for i, b in enumerate(r) if b]