How to divisor sieve numbers of form kn^p? [closed]

Question

Closed. This question is off-topic. It is not currently accepting answers.

---

Want to improve this question? Update the question so it's on-topic for Stack Overflow.

Closed 8 years ago.

Improve this question

Creating a divisor-counting sieve for all numbers 1-n is well known

func sieve(n)
    counts = array of size n+1, all elements set to 1, counts[0] = 0 arbitrary

    for 2 <= i <= n
        for i <= j <= n, stepwise i
            counts[j]++;
    return counts

However what if, instead of creating a sieve for numbers of form 1 * n, I wanted the divisor counts for numbers of form 6n^2 instead?

So instead of finding the divisor counts for 1, 2, 3, 4, 5, etc It might be instead looking for divisor counts of 6, 24, 54, 96, 150, etc

But really just numbers of form kn^p, in an efficient way so I am not actually storing an array of size kn^p at its largest. Seems like I should only need array of size N as before, only each spot represents the number of divisors of kn^p

Answer 1

Instead of counting the divisors directly, you can use this formula (taken from Wikipedia):

This will help you find the number of divisors for all square/cube/to the power of k numbers (n^k) while using O(n) memory.

// Stores the number of divisors for n^k
count[n + 1]

// Count number of divisors for 1^k, 2^k, ..., n^k
function countDivisors(n, k) {
    // Note: You can actually merge sieve function with the loop below
    prime[] = sieve(n)

    count[0] = 0
    count[1..n] = 1 // Set count[i] to count[n] to 1

    for all primes p <= n {
        for (i = 1; p^i <= n; i++) {
            // You can cache the value of p^i for next loop
            pi = p^i 

            for (m = pi; m <= n; m += pi) {
                // We "replace" the previous count with current count
                count[m] = count[m] / ((i - 1) * k + 1) * (i * k + 1)
            }
        }
    }
}

To extend it to a*n^k, find the prime factorization of a, and record the prime numbers along with the corresponding power.

• If any prime is larger than n, you can take them away and reduce them to a multiplier according to the divisor function.

• If any prime is smaller than n, supply it to the countDivisors function above, and modify the function a bit:

  for all primes p <= n {
      // Note: You are free to optimize this in actual code
      let e be maximal power of p in a
  
      if (e > 0) {
          // Update all number with the factors from a
          for (i = 1; i <= n; i++) {
              count[i] *= (e + 1)
          }
      }
  
      for (i = 1; p^i <= n; i++) {
          // You can cache the value of p^i for next loop
          pi = p^i 
  
          for (m = pi; m <= n; m += pi) {
              // We "replace" the previous count with current count
              count[m] = count[m] / ((i - 1) * k + e + 1) * (i * k + e + 1)
          }
      }
  }
  

As you can see, the time and space complexity does not depend on the power k, and only depends on n, the number of numbers whose number of divisors you want to find.