Prime power mathematical problem in Python

Question

I have a mathematical problem about prime power. A number \$x\$ can be considered a prime power if \$x = p^k\$ where \$p\$ is a prime and \$k\$ is a non-negative integer. For example, \$81\$ is a prime power because \$81 = 3^4\$. Now, form a sequence of numbers in the following way. Start by taking a random number.

1. If the number is a prime power, end the sequence.
2. If the number is not a prime power, minus that number with the biggest prime power that's not more than the number itself. Do it over and over again until a prime power is met, then stop.

For example, 34 is not a prime power. Hence we minus it with the highest prime power that's not more than it which is \$32\$, hence \$34-32=2\$. \$2\$ is a prime power, so we stop. In this case, the length of the operation is 2 (because 34 -> 2).

Another example is \$95\$. \$95\$ isn't a prime power, so \$95-89=6\$. \$6\$ is also not a prime power, so we minus it again with the highest prime power, \$6-5=1\$. \$1\$ is a prime power (because \$2^0 = 1\$), so we stop. In this case, the length of the operation is 3 (because 95 -> 6 -> 1).

Known that:

• 1 is the smallest initial number that can be formed with a length of 1 operation.
• 6 is the smallest initial number that can be formed with a length of 2 operations.
• 95 is the smallest initial number that can be formed with a length of 3 operations.
• 360748 is the smallest initial number that can be formed with a length of 4 operations.

I made an application using Python to find the smallest initial number that can be formed with a length of 5, 6 and 7 operations. It is using looping to solve the problem. However, to find the smallest initial number that can be formed with a length of 4 or more operations, the program takes a very very long time to run.

Is there any way I could improve my code to make the program to run much faster? I just want to focus on getting the number with 5, 6 and 7 operations.

from sympy.ntheory import factorint
def pp(q): #to check whether it is a prime power
    fact=factorint(q)
    p=int(list(fact.keys())[0])
    n=int(list(fact.values())[0])
    if q!=p**n:
        return False
    else:
        return True

a=[1,6][-1]

b=a*2
d=1
while d!=0:
    c=b-a
    if pp(b)==False and pp(c)==True:
        for i in range(b-1,c-1,-1):
            if i==c:
                d=0
            elif pp(i)==True:
                b=i+a
                break
    elif pp(b)==False:
        for i in range(b-1,c-1,-1):
            if i==c:
                b=b+a+1
            elif pp(i)==True:
                b=i+a
                break
    elif pp(b)==True:
        b=b+a
b

Answer 1

Welcome, there are some aspects to comment about.

• First is that, you should use specific names for your functions pp is ok and sounds quite interesting, but it would be better if you use is_prime_power (as mentioned by Konrad Rudolph in the comments).

• Now, for the first function, I would change only two things:

  def is_prime_power(q):
      fact = factorint(q)
      p = int(list(fact.keys())[0])
      n = int(list(fact.values())[0])
      return q == p**n # this line asks 'is prime power?'
  
  # could also be:
  def is_prime_power(q):
      return q == int(list(fact.keys())[0])**int(list(fact.values())[0])
  # but is quite strange because p,n are not specified.
  
  # this inverses the boolean answer, as you see
      if q != p**n:
          return False
      else:
          return True
  

• For the procedure, you should avoid using comparisons with True and False. If the value is a boolean, Python deduces if it is true or false, in the case of compare with False use not:

  a = [1,6][-1]
  
  b = a*2
  d = 1
  while d != 0:
      c = b-a
      if not is_prime_power(b) and is_prime_power(c):
          for i in range(b-1, c-1, -1):
              if i == c:
                  d = 0
              elif is_prime_power(i):
                  b = i+a
                  break
      elif not is_prime_power(b):
          for i in range(b-1, c-1, -1):
              if i == c:
                  b = b+a+1
              elif is_prime_power(i):
                  b = i + a
                  break
      elif is_prime_power(b):
          b = b+a
  b
  

Well, thanks for writing easy to read code, it is so good to see something like it.

Answer 2

The question seems to be related to prime gaps.

Some observations:

1. Since we are trying to find the smallest number formed with some length operations, that means it is the first number formed with that length operations if we start our search from 1.

2. For the number 6, the biggest prime power less than 6 is 5 i.e. 6-5 = 1

   For the number 95, the biggest prime power less than 95 is 89 i.e. 95-89 = 6

   For the number 360748, the biggest prime power less than 360748 is 360653 i.e. 360748-360653=95.

This shows that for the next number that will be formed using 5 length operations, the gap between that number and biggest prime power less that number must be 360748. The problem is simplified to finding the smallest number with has a prime power gap of 360748 and so on.

---

Now consider the list of prime numbers and list of prime power numbers. It is obvious that prime power numbers are more frequent that prime numbers. So it is safe to assume that the gaps between prime power numbers will be less than gaps between prime numbers.

Straight from Wikipedia article prime gaps,

As of August 2018 the largest known maximal prime gap has length 1550, found by Bertil Nyman. It is the 80th maximal gap, and it occurs after the prime 18361375334787046697.

From this post, we can make a reasonable estimate for the prime number with gap 360748 i.e. exp(sqrt(360748)) ~ 7*10^260. (There are better estimations available in literature).

Too big and not suitable for our Personal Computers. In general sieves are faster than checking primes directly. It is easy to make a sieve for prime power numbers, however it is impractical for current problem. To make a boolean array of length 10^260, we will need around 10^247 TeraByte.

Even if we somehow manage to determine the result for length 5, for the next number in sequence, it should have a prime power gap of 10^260. We will probably have to wait for Personal Quantum computers for that.

---

EDIT: Following is the code to generate a prime power sieve.

def prime_power_sieve(limit):
    sieve = [1]*limit
    sieve[0] = 0
    p = 2
    while p<limit:
        if sieve[p]>0:
            # mark multiples of p
            for j in range(p*p, limit, p):
                sieve[j] = 0
            # mark powers of p
            j = p
            while j*p<limit :
                j *= p
                sieve[j] = 1 
        p += 1
    return sieve

To find largest prime power lower than or equal to that number, we can apply a simple modification by tracking the previously encountered prime power.

def prime_power_sieve_with_previous(limit):
    sieve = list(range(0,limit))
    p = 2
    prev = p
    while p<limit:
        if sieve[p]>0:
            prev = p
            # mark multiples of p
            for j in range(p*p, limit, p):
                sieve[j] = 0
            # mark powers of p
            j = p
            while j*p<limit :
                j *= p
                sieve[j] = j
        elif sieve[p] == 0:
            sieve[p] = prev
        p += 1
    return sieve

With this, largest prime power less than or equal to p will be sieve[p].

print(prime_power_sieve(20))
print(prime_power_sieve_with_previous(20))

Output:
[0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1]
[0, 1, 2, 3, 4, 5, 5, 7, 8, 9, 9, 11, 11, 13, 13, 13, 16, 17, 17, 19]