Find sum of all numbers whose largest prime factor is less than or equal to P and has K number of prime factors.

Question

Given a prime P and an integer K denoting the number of prime factors(not necessarily distinct), Find sum of all numbers whose largest prime factor is less than or equal to P and has exactly K prime factors.

For example P=5 and K=4,
number 40(ie = 2*2*2*5) is one such candidate.
number 100(ie = 2*2*5*5) is another such candidate

I need to find the sum of all such numbers in as efficent way as possible.

To address the demands to show my efforts, here you go:

Here is my brute force for the above example

import sympy

P = 5
K = 4

s = 0
for i in range(2, P**K + 1):
    all_facs = sympy.factorint(i, multiple=True)
    if all_facs[-1] > P:
        continue
    if len(all_facs) != K:
        continue
    print(i, all_facs)
    s += i
print("ans=", s)

# 16 [2, 2, 2, 2]
# 24 [2, 2, 2, 3]
# 36 [2, 2, 3, 3]
# 40 [2, 2, 2, 5]
# 54 [2, 3, 3, 3]
# 60 [2, 2, 3, 5]
# 81 [3, 3, 3, 3]
# 90 [2, 3, 3, 5]
# 100 [2, 2, 5, 5]
# 135 [3, 3, 3, 5]
# 150 [2, 3, 5, 5]
# 225 [3, 3, 5, 5]
# 250 [2, 5, 5, 5]
# 375 [3, 5, 5, 5]
# 625 [5, 5, 5, 5]
# ans= 2261

And here is my mathematics for a better $O(P.K)$ algorithm. Let S(P,K) be the desided answer and Let SS(P,K) be sum when the biggest prime is strictly equal to P and $P_{prev}$ be the prime just before prime P

$$S(P,K) = SS(P,K) + S(P_{prev},K) --------------------(1)$$

Lets eunmerate over each case and add the sum to $SS(P,K)$

case1:place all Ps for all K slots $$ P^K $$

case2:place all Ps for all last K-1 slots $$ P^{K-1} . S(P_{prev}, 1) $$

case3:place all Ps for all last K-2 slots $$ P^{K-2} . S(P_{prev}, 2) $$

$...$

$...$

$...$ caseK:place Ps for all last 1 slots $$ P^{1} . S(P_{prev}, K-1) $$

Adding all cases:

$$SS(P,K) = P^K + P^{K-1} . S(P_{prev},1) + P^{K-2} . S(P_{prev},2) + ... + P^1 . S(P_{prev},K-1) -------------------------(2)$$ putting (2) into (1),

$$ S(P, K) = P^K + P^{K-1} . S(P_{prev}, 1) + P^{K-2} . S(P_{prev}, 2) + ... + P^1 . S(P_{prev}, K-1) + S(P_{prev}, K) $$

The above recurence can be solved in $O(K . no\_of\_primes\_under\_P)$

So from $O(P^{5K/4})$ we are down to $O(PK)$

Anyone having better suggestions are welcome

Answer 1

Edit 2: OK, Wow. A quick search on OEIS turns up this lovely sequence for $k=5$. Further searching suggests that the sequences $F(p,k)$ for all $p$ have the generating function

$$\prod^{p_i \in \mathbb{P}}_{2 \le p_i \le p} \frac1 {(1-p_ik)}$$

So, for example, $F(5, k)$ has the OGF: $1/(1-2x)(1-3x)(1-5x)$, as noted in the linked OEIS sequence. So... that's the real answer for an efficient method.

---

Not a full answer, but it's too long for a comment.

Define the recursive function $F(p, k, m)$, where $m$ is the smallest prime factor used in the numbers being examined. Then $F(p,k) = \sum_{2 \le m \le p}^{m \in \mathbb{P}} F(p, k, m)$. Now we can write the recursion

$$F(p, k, m) = m \sum_{m \le i \le p}^{i \in \mathbb{P}} F(p, k-1, i)$$

And so

$$F(p, k) = \sum_{2 \le m \le p}^{m \in \mathbb{P}} \left( m \sum_{m \le i \le p}^{i \in \mathbb{P}} F(p, k-1, i) \right)$$

For example, using $p=5$ as you did:

$$ \small \begin{array}{l l l l} F(5, 1, 2) = 2 & F(5,1,3) = 3 & F(5,1,5) = 5 & F(5,1) = 10 \\ F(5, 2, 2) = 2 \cdot 10 = 20 & F(5,2,3) = 3 \cdot 8 = 24 & F(5,2,5) = 5 \cdot 5 = 25 & F(5,2) = 69 \\ F(5, 3, 2) = 2 \cdot 69 = 138 & F(5,3,3) = 3 \cdot 49 = 147 & F(5,3,5) = 5 \cdot 25 = 125 & F(5,2) = 410 \end{array} $$

I don't know algorithmic complexity calculations very well, but this requires you hold in memory only (A) the list of primes $(2, 3, \cdots p)$ and (B) the list of values of $F(p, k-1,m)$. Each new loop requires $\pi(p)$ multiplications. As $p$ gets large, the total multiplications for a given $k,p$ tends toward $kp/\log p$. I suspect SageMath's builtin functions and python's list comprehensions will make the calculations significantly faster.

Edit: Some quick code.

import copy
p = 997
k = 1000
Output = []
PList = [i for i in [2..p] if i in Primes()]
OldFList = [0 for i in PList]
NewFList = [0 for i in PList]
for m in [1..k]:
    for n in [0..len(PList)-1]:
        if OldFList[n] == 0:
            NewFList[n] = PList[n]
        else:
            NewFList[n] = PList[n] * sum(OldFList[n:])
    OldFList = copy.deepcopy(NewFList)
print("F("+str(p)+","+str(m)+") = "+str(sum(NewFList)))

This gave $F(997, 100)$ in $0.4$ s, $F(997, 1000)$ in $8.5$ s, $F(19997, 100)$ in $37$ s using SageMath/Python. The middle one has $3000+$ digits.