Perfect Power Sequence

February 16, 2018

I discovered the Numberphile channel on YouTube a few months ago, and have greatly enjoyed its videos at the intersection of mathematics and programming. Their recent video discusses the Catalan Conjecture:

In the infinite sequence of perfect powers 1, 4, 8, 9, 16, 25, 32, 36, 49, 64, 81, 100, …, the only two adjacent numbers are 2³ = 8 and 3² = 9.

The conjecture has recently been proved, as discussed at Numberphile.

Your task is to write a program that generates the perfect power sequence. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Advertisements

Pages: 1 2

12 Responses to “Perfect Power Sequence”

  1. Zack said

    Here is my implementation in Julia.

    function IsPerfectPower(x::Int64)
    y = ceil(Int64, sqrt(x))

    for z = 2:y
        if mod(x ^ (1 / z), 1) == 0; return true; end
    end
    
    return false
    

    end

    function genseq(n::Int64)
    y = zeros(Int64, n)
    y[1] = 1

    for x = 2:n
        if IsPerfectPower(x); y[x] = x; end
    end
    
    return y[y .> 0]
    

    end

    Performance benchmark: on an Intel i7 laptop running on Win7, it generates the sequence of the PP numbers that are less or equal to 10000 in about 0.1 sec, using approx. 85 MB of RAM.

  2. Zack said

    Typo in my previous commend: 85 Kilobytes, not MB… Sorry about that!

  3. matthew said
    import Data.List.Ordered
    f n = n*n : merge powers (f (n+1)) where powers = n*n*n : map (n*) powers
    main = print (take 20 (nub (1 : (f 2))))
    -- [1,4,8,9,16,25,27,32,36,49,64,81,100,121,125,128,144,169,196,216]
    
  4. Paul said

    Two solutions in Python. I expected the first to be the fastest: merging already sorted sequences. But it turned out to be otherwise. The second is a little faster.

    from ma.mymath import isqrt  # integer square root
    from itertools import takewhile, count
    from functools import partial
    from operator import ge
    
    def genpows(n):
        """generate powers n^2, n^3, ...."""
        x = n
        while True:
            x *= n
            yield x
    
    def perpow(limit):
        """generate all powers (up to limit), merge them and handle duplicates"""
        maxn, last = isqrt(limit), 1
        its = [genpows(i) for i in range(2, maxn+1)]
        yield 1
        for i in takewhile(partial(ge, limit), merge(*its)):
            if i != last:
                yield i
                last = i
    
    def perpow2(limit):
        def nums(i):
            return takewhile(partial(ge, limit), 
                map(partial(pow, i), count(2)))
        powers = set((n for i in range(2, isqrt(limit)+1) for n in nums(i)))
        return [1] + sorted(powers)
    
  5. Scott said

    Here is some Java. It generates up to 10000 in around 60ms on a 3.3GHz i5 iMac:

    import java.util.PriorityQueue;
    
    public class PerfectPowerSequenceGenerator {
      private static class NumberSequence {
        public long number;
        public long nextPerfectPowerValue;
        
        public NumberSequence(long number) {
          this.number = number;
          this.nextPerfectPowerValue = number * number;
        }
      }
    
      private PriorityQueue<NumberSequence> sequences = new PriorityQueue<>((a, b) -> {
        if (b.nextPerfectPowerValue > a.nextPerfectPowerValue) {
          return -1;
        } else if (b.nextPerfectPowerValue < a.nextPerfectPowerValue) {
          return 1;
        }
        return 0;
      });
    
      private long nextSequence;    // next sequence to add to the heap
      private long lastResult = 0;  // last result we returned, to avoid dupes
      
      public PerfectPowerSequenceGenerator() {
        nextSequence = 2;
      }
    
      public long next() {
        long result = 1;
        if (lastResult != 0 && !sequences.isEmpty()) {
          NumberSequence seq;
          do {
            seq = sequences.poll();
            result = seq.nextPerfectPowerValue;
            seq.nextPerfectPowerValue *= seq.number;
            sequences.add(seq);
          } while (result == lastResult);
        }
            
        if (nextSequence * nextSequence > result) {
          sequences.add(new NumberSequence(nextSequence++));
        }
        
        lastResult = result;
      
        return result;
      }
      
      public static void main(String[] args) {
        long startTime = System.nanoTime();
        PerfectPowerSequenceGenerator gen = new PerfectPowerSequenceGenerator();
        while (true) {
          long pp = gen.next();
          System.out.println(pp);
          if (pp >= 10000) {
            break;
          }
        }
        System.out.println("Ran in " + (System.nanoTime() - startTime) / 1000000 + "ms");
      }
    }
    
  6. matthew said

    The queue solution works nicely in Python:

    import heapq
    def gen(i0,j0,f):
        global q
        q = [(f(i0,j0),i0,j0)]
        last = None
        while True:
            (n,i,j) = heapq.heappop(q)
            heapq.heappush(q,(f(i+1,j),i+1,j))
            if i == i0: heapq.heappush(q,(f(i,j+1),i,j+1))
            if n != last: yield n
            last = n
    
    g = gen(2,2,lambda i,j: i**j)
    s = [1]+[next(g) for _ in range(19)]
    print(s); print(q)
    # [1,4,8,9,16,25,27,32,36,49,64,81,100,121,125,128,144,169,196,216]
    # [(225,15,2),(243,3,5),(256,2,8),(256,4,4),(729,3,6),(2187,3,7),(343,7,3)]
    

    It’s quite important how the queue is extended – doing

            if j == j0: heapq.heappush(q,(f(i+1,j),i+1,j))
            heapq.heappush(q,(f(i,j+1),i,j+1))
    

    still works, but the queue gets much longer.

    Incidentally, the code on the solution page seems to be garbled – I think WordPress have messed something up with angle brackets in code.

  7. Paul said

    @matthew: I tried the priority method too and it worked nicely, but a little slow in Python. The method below takes about 0.1 sec for all powers below 10^9.

    from ma.mymath import isqrt  # integer square root
    from itertools import takewhile, count
    from functools import partial
    from operator import ge
    
    _small_pp = [1, 4, 8, 9, 16, 25, 27, 32, 36, 49, 64, 81, 100]
    
    def pp2(limit):
        def powers(i):
            return takewhile(partial(ge, limit), 
                map(partial(pow, i), count(2)))
        nmax = isqrt(limit)
        # remove the perfect powers from the range of bases to loop over
        if nmax > 120:  # 120 is the first pp after 100 (_small_pp)
            excl = pp2(nmax)
        else:
            excl = _small_pp
        bases = set(range(2, nmax+1)) - set(excl)
        allpowers = (n for base in bases for n in powers(base))
        return [1] + sorted(allpowers)
    
  8. matthew said

    @Paul: did you do a speed comparison with my priority queue function above?

  9. Paul said

    @matthew: I see now that your method works very fast (91 msec for yours and 121 ms for pp2 for all powers up to 10^9). I must have messed up something with my implementation.

  10. matthew said

    @Paul: that’s about what I made it – also the queue uses less memory, I’m running the program just printing out all the numbers, up to 25528390307779969 after about 30 mins.

    For the queue, as mentioned above, it’s significant how the queue is extended.

  11. Daniel said

    Here’s an implementation in C that uses the approach from the solution.

    #include <stdbool.h>
    #include <stdio.h>
    #include <stdlib.h>
    
    typedef unsigned long u_long;
    
    typedef struct {
      u_long base;
      u_long exponent;
      u_long value;
    } item_t;
    
    typedef struct {
      item_t* array;
      size_t count;
      size_t capacity;
    } heap_t;
    
    typedef struct {
      heap_t* heap_p;
      size_t n_generated;
    } generator_t;
    
    void init_heap(heap_t** heap_pp) {
      *heap_pp = malloc(sizeof(heap_t));
      heap_t* heap_p = *heap_pp;
      heap_p->array = malloc(sizeof(item_t));
      heap_p->count = 0;
      heap_p->capacity = 1;
    }
    
    void clear_heap(heap_t** heap_pp) {
      free((*heap_pp)->array);
      free(*heap_pp);
      *heap_pp = NULL;
    }
    
    void swap(item_t* a, item_t* b) {
      item_t tmp = *a;
      *a = *b;
      *b = tmp;
    }
    
    void heap_insert(heap_t* heap_p, item_t item) {
      size_t new_count = heap_p->count + 1;
      if (new_count > heap_p->capacity) {
        size_t new_capacity = 2 * (heap_p->capacity + 1) - 1;
        heap_p->array = realloc(heap_p->array, new_capacity * sizeof(item_t));
        heap_p->capacity = new_capacity;
      }
      item_t* array = heap_p->array;
      // reheapify upward
      size_t idx = heap_p->count;
      array[idx] = item;
      while (1) {
        if (idx == 0) break;
        size_t parent_idx = (idx-1)/2;
        u_long value = item.value;
        u_long parent_value = array[parent_idx].value;
        if (value > parent_value) break;
        swap(&array[idx], &array[parent_idx]);
        idx = parent_idx;
      }
      heap_p->count = new_count;
    }
    
    item_t heap_peak(heap_t* heap_p) {
      return (heap_p->array)[0];
    }
    
    // Popping from an empty heap is not supported
    item_t heap_pop(heap_t* heap_p) {
      item_t popped_item = heap_peak(heap_p);
      size_t new_count = heap_p->count - 1;
      if (new_count > 0) {
        item_t* array = heap_p->array;
        array[0] = array[new_count];
        size_t idx = 0;
        while (1) {
          u_long value = array[idx].value;
          size_t left_idx = 2 * idx + 1;
          size_t right_idx = left_idx + 1;
          if (left_idx > new_count - 1) break;
          u_long min_child_value = array[left_idx].value;
          size_t min_child_idx = left_idx;
          if (right_idx <= new_count - 1) {
            u_long right_child_value = array[right_idx].value;
            if (right_child_value < min_child_value) {
              min_child_value = right_child_value;
              min_child_idx = right_idx;
            }
          }
          if (value <= min_child_value) break;
          swap(&array[idx], &array[min_child_idx]);
          idx = min_child_idx;
        }
      }
      heap_p->count = new_count;
      return popped_item;
    }
    
    void init_generator(generator_t** generator_pp) {
      *generator_pp = malloc(sizeof(generator_t));
      generator_t* generator_p = *generator_pp;
      init_heap(&(generator_p->heap_p));
      item_t item = {.base = 2, .exponent = 2, .value = 4};
      heap_insert(generator_p->heap_p, item);
      generator_p->n_generated = 0;
    }
    
    void clear_generator(generator_t** generator_pp) {
      clear_heap(&(*generator_pp)->heap_p);
      free(*generator_pp);
      *generator_pp = NULL;
    }
    
    u_long next(generator_t* generator) {
      ++(generator->n_generated);
      if (generator->n_generated == 1) return 1;
      heap_t* heap_p = generator->heap_p;
      u_long value = heap_peak(heap_p).value;
      while (heap_peak(heap_p).value == value) {
        item_t item = heap_pop(heap_p);
        if (item.exponent == 2) {
          item_t new_item;
          new_item.base = item.base + 1;
          new_item.exponent = 2;
          new_item.value = new_item.base * new_item.base;
          heap_insert(heap_p, new_item);
        }
        ++(item.exponent);
        item.value *= item.base;
        heap_insert(heap_p, item);
      }
      return value;
    }
    
    int main(void) {
      generator_t* generator;
      init_generator(&generator);
      u_long limit = 1764;
      while (1) {
        u_long value = next(generator);
        printf("%lu", value);
        if (value >= limit) break;
        printf(",");
      }
      clear_generator(&generator);
      return 0;
    }
    

    Output:

    1,4,8,9,16,25,27,32,36,49,64,81,100,121,125,128,144,169,196,216,225,243,256,289,324,343,361,400,441,484,512,529,576,625,676,729,784,841,900,961,1000,1024,1089,1156,1225,1296,1331,1369,1444,1521,1600,1681,1728,1764
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s

%d bloggers like this: