Counting integers $n \leq x$ with a given prime signature

Question

Given is a prime signature $S$ and an integer $x$. The task is to count how many integers $n$ exist such that $n \leq x$, and if $n = p_1^{k_1}p_2^{k_2}p_3^{k_3}p_4^{k_4}...$ then $S = (k_1,k_2,k_3,...)$. Let $C(S,x)$ be that count.

For example, if $S = (3,1,1)$ and $x=1000$ then there are $21$ such integers, being $$120, 168, 264, 280, 312, 408, 440, 456, 520, 552, 616, 680, 696, 728, 744, 760, 888, 920, 945, 952, 984$$

And so $C\big((3,1,1),1000\big )=21$. Notice that the order matters, so $270=2\cdot 3^3\cdot 5$ does not appear in this specific list. It would appear for $S=(1,3,1)$ for large enough $x$. The signature consists only of positive integers, and the primes are sorted in increasing order, $p_i<p_j$ for all $i<j$.

Currently I have thought of just $1$ approach for this:

1. First, use recursion to generate such integers, keeping track of the currently generated integer $n$, and the index of the largest prime used, $i$.
2. Once the last part of $S$ is reached, call it $k_m$, count the number of primes smaller than $\sqrt[k_m]{\frac xn}$ and make sure only primes larger than $p_i$ are included. That is, add $\pi \Big(\sqrt[k_m]{\frac xn}\Big) - \pi(p_i)=\pi \Big(\sqrt[k_m]{\frac xn}\Big) - i$ to the total count.

It may be assumed the prime counting function implementation used for this algorithm is very efficient. Still, I have to compute $C(S,x)$ for thousands of signatures, and for $x=10^{12}$, so either my implementation of the above algorithm is lacking, or there is simply a better way to compute $C(S,x)$.

This is part of a recreational maths problem, but is in itself a nice problem, so any approaches to calculate the function efficiently are welcome.

Answer 1

Here is a complicated approach that might offer a modest improvement for values of $x$ of the size mentioned in your comment, when you want to compute $C(S,x)$ for many different $S$ and the same $x$:

---

Let me focus first on the case where $S$ has the form $S=(\cdots,1,1)$.

Let $f(u,v)$ be the number of integers of the form $p_1 p_2$ satisfying $u \le p_1<p_2$ and $p_1 p_2 \le v$. We can precompute a table that lets us compute $f(u,v)$ very rapidly whenever $u\le 10^3$, $v \le 10^6$.

A naive approach is to store $f(u,v)$ for all $u,v$ satisfying $u \le 10^3$, $v \le 10^6$. This can be stored in an array with $10^3 \times 10^6 = 10^9$ elements. We can fill in the array using the recurrence $f(u,v) = f(u+1,v)$ if $u$ is not prime, or $f(u,v) = f(u+1,v) + 1$ if $u$ is prime and $v/u$ is prime and $v>u^2$. In particular, for each $v$, start by computing $f(1000,v)$, then $f(999,v)$, then $f(998,v)$, etc.

A more efficient approach is to store $f(u,v)$ only when $u \le 10^3$ is prime and $v \le 10^6$ is a product of two different primes. If we build a list of all products of two different primes $\le 10^6$, given any $u,v$, we can find the smallest prime $u'\ge u$ and the largest product of two primes $v' \ge v$ (via binary search in the list) and finally look up the value $f(u',v')$, using the fact that $f(u,v)=f(u',v')$. It will be possible to fill in the lookup table with a few million operations.

Once we've built the lookup table, it is now possible to offer an improved algorithm for computing $C(S,x)$ whenever $S$ has the form $S=(\cdots,1,1)$. In particular, use the same algorithm as in the question, but once the last two parts of $S$ are reached, if $n \ge x/10^6$ and $q\le 10^3$ where $q$ was the prime used in the third-to-last part, then compute $f(q,x/n)$, add $f(q,x/n)$ to the total count, and end the recursion. (If $n<10/x^6$, continue on with the recursion as in your algorithm.)

You can optimize the parameters $10^3,10^6$ above to speed the algorithm up, based on the amount of storage available and the distribution of sizes of $x/n$. (I doubt $10^3,10^6$ are optimal; I chose them more or less arbitrarily to try to make the approach easy to understand, not because I believe those are the best parameters.)

---

Now, if you have multiple signatures all of the form $S=(\cdots,1,1)$ varying in their first parts, then you can reuse the same table $f$ for all of them.

In a similar way you can build another table for all signatures of the form $S=(\cdots,1,2)$ or $S=(\cdots,2,1)$ or any other common suffix.

---

You can also take this a step further, by building an analogous lookup table to help with all $S$ where $S=(\cdots,1,1,1)$.

Let $g(t,v)$ be the number of integers of the form $p_1p_2p_3$ satisfying $t \le p_1 <p_2<p_3$ and $p_1p_2p_3 \le v$. We'll precompute a lookup table that stores $g(t,v)$ for all $t,v$ where $t \le 10^3$, $v \le 10^9$, where $t$ is prime and $v$ is a product of three distinct primes.

To fill in the lookup table, we'll use the fact that $g(t,v) = g(t^*,v) + f(t^*,v/t^*)$, where $t^*$ is the smallest prime with $t^*>t$. In this way, for each $v$, we can fill in $g(t,v)$ in the order $g(997,v), g(991,v), g(983,v), g(977,v), \dots$. The lookup table will have fewer than $10^9$ entries and can be built with less than $10^9$ steps, so it is probably feasible to build it in a few seconds if implemented appropriately.

Once we've built the lookup table, we can use it to improve the algorithm for computing $C(S,x)$ whenever $S$ has the form $S=(\cdots,1,1,1)$. In particular, use the same algorithm as above, but once the last three parts of $S$ are reached, if $n \ge x/10^9$ and $q\le 10^3$ where $q$ was the prime used in the fourth-to-last recursion, then compute $g(q,x/n)$, add $f(q,x/n)$ to the total count, and end the recursion. (Otherwise, continue the algorithm as above).

Now you can reuse the $g$ lookup table for all $S$ that have the form $S=(\cdots,1,1,1)$.

You can optimize the parameters $10^3,10^9$ for performance, and do something similar for any other common suffix.