# Generation of products of power

Given positive integers a,b,c and limit, I want to generate all products having the form $$p^aq^br^c \leq limit$$ less than limit (where p,q,r are distinct primes). a,b and c are not necessary distinct.

I tried something like:

primes= [ 2,3,5,...] # primes up to 10**8
[ (p**a)*(q**b)*(r**c) for p in primes for q in primes for r in primes if (p**a)*(q**b)*(r**c) <= limit ]


But it is very slow because len(primes)(=5761455) is high.

Then I tried a very ugly code for printing values. It generates all values for p < q < r (p, q, r primes)

def generate3(a,b,c,limit):
global primes
i1 = 0
i2 = 1
i3 = 2
while i1 < len(primes) and i2 < len(primes) and i3 < len(primes) and (primes[i1]**a)*(primes[i2]**b)*primes[i3]**c < limit:
print(str(i1)+" "+str(i2)+str(i3))
while i1 < len(primes) and i2 < len(primes) and i3 < len(primes) and (primes[i1]**a)*(primes[i2]**b)*primes[i3]**c< limit:
while i1 < len(primes) and i2 < len(primes) and i3 < len(primes) and (primes[i1]**a)*(primes[i2]**b)*primes[i3]**c< limit:
print(str(primes[i1])+" "+str(primes[i2])+" "+str(primes[i3])+" "+str((primes[i1]**a)*(primes[i2]**b)*(primes[i3]**c)))
i3 = i3 + 1
i2 = i2 +1
i3 = i2 +1
i1 = i1 + 1
i2 = i1 + 1
i3 = i2 + 1


Is there a more pythonic way to do this? Is there a better method for generating these products of powers quickly?

• Is the constraint p < q < r  a coincidence, or one of the do they just need to be distinct? – Maarten Fabré Nov 16 '18 at 14:49
• They just need to be distinct but the second example I gave only generates cases for which p < q < r. – montardon Nov 16 '18 at 16:14

# global

there is no need to make primes a global. You only read from it, but don't assign to it, you can use it as is. It will be even faster if you make primes a local variable by passing it in as a parameter, so Python uses the LOAD_FAST bytecode instead of LOAD_GLOBAL. Since primes is called, indexed and sliced a lot, this can make a difference.

# while condition

i1 < len(primes) and i2 < len(primes) and i3 < len(primes) and .... Since i1<i2<i3, only i3 < len(primes) is needed. If you use len(primes) so often, it pays to make it a local variable.

# return, don't print

your method immediately prints the results. In general it is better to split the calculation and presentation, so to let the method return or yield the values, and another method do the presentation

# looping

I suggest you watch the talk Loop like a Pro by David Baumgold. It's recommended material for every python programmer.

Instead of looping over the indices, you can loop over the primes-list immediately.

For p, and the first index (i), you can loop over enumerate(primes)

Then you can use this index i to slice primes to only include the elements with an index larger than p and go on to r, so you arrive at the following, naive implementation:

def generate4_naive(a, b, c, limit):
for i, p in enumerate(primes):
for j, q in enumerate(primes[i+1:], i+1):
for r in primes[j+1:]:
product = p**a * q**b * r**c
if product < limit:
yield p, q, r, product


If you include the early breaks, you arrive at something like this:

def generate4(a, b, c, limit):
for i, p in enumerate(primes):
for j, q in enumerate(primes[i+1:], i+1):
sub_product = p**a * q**b
for r in primes[j+1:]:
product = sub_product * r**c
if product > limit:
break
yield p, q, r, product
if sub_product * primes[j+1]**c > limit:
break
if p ** a * primes[i+1] **b * primes[i+2] ** c > limit:
return


then I also used primes as a local variable by changing the signature: def generate4b(a, b, c, limit, primes):

# Performance:

limit = 1_000_000
%timeit tuple(generate3(2,3,4, limit=limit))

16.6 µs ± 166 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit tuple(generate3b(2,3,4, limit=limit, primes=primes))

16.2 µs ± 138 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit tuple(generate4(2,3,4, limit=limit))

14.2 µs ± 643 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit tuple(generate4b(2,3,4, limit=limit, primes=primes))

7.72 ms ± 516 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit list(generate4_naive(2,3,4, 1_000)) # primes also only to 1000

940 ms ± 84.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)