List the positive integers in Cantor unpacking order of their prime factorizations.

cantor_order_integers.py

from numpy import prod

# Can we reorder the positive integers based on a logical ordering of their prime factorizations?

# If we could iterate through all finite-sized tuples of nonnegative integers, we could increment
# the first element of each to uniquely cover all possible prime factorizations.

# We can extend and invert the Cantor pairing function to generate all k-tuples for any positive k,
# but now we're only iterating through the integers with highest prime factor p(k).

# So let's use the inverse pairing function again to generate all 2-tuples! For each integer we
# get back a 2-tuple (x, y), which we'll use to generate the xth y-tuple and transform that into
# a prime factorization.

# TODO: For now this uses the Szudzik multidimensional unpacking function instead of Cantor.

# 1, 2, 3, 4, 9, 5, 25, 8, 6, 15, 7, 49, 35, 16, 18, 75, 245, 11, 121, 77, 847, 32, 27, 10, 21, 55, 13, 169, 143, 1859...

def test():
    primegen = gen_primes()
    primes = list()
    integers = [1]
    for i in range(0, 29):
        # choose index and dimension
        (n, k) = multidimensional_szudzik_unpairing(2, i)
        k += 1
        # generate prime factorization
        exponents = multidimensional_szudzik_unpairing(k, n)
        exponents[-1] += 1
        while k > len(primes):
            primes.append(next(primegen))
        factorization = zip(primes, exponents)
        # multiply prime factors to get the corresponding integer
        integer = prod([pair[0] ** pair[1] for pair in factorization])
        integers.append(integer)
    print(f'First {len(integers)} integers in unpacking order:')
    print(integers)
    first_missing = next(i + 1 for i,v in enumerate(sorted(integers)) if i != v - 1)
    print('First integer absent from sequence:')
    print(first_missing)

# thanks to David Eppstein: https://code.activestate.com/recipes/117119/
def gen_primes():
    D = {}  
    q = 2  
    while True:
        if q not in D:
            yield q        
            D[q * q] = [q]
        else:
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]
        q += 1

# the rest thanks to David R. Hagen: https://github.com/drhagen/pairing
def multidimensional_szudzik_unpairing(n, index):
    shell = floored_root(index, n)

    def recursive_indexes(dim, remaining):
        if dim == n - 1:
            # If this is reached, that means that no index so far has been on the shell. By construction, this
            # final index must be on the shell or else the point itself will not be on the shell.
            return [shell]
        else:
            # Number of dimensions of a slice perpendicular to the current axis
            slice_dims = n - dim - 1

            # Number of elements in a slice (only those on the current shell)
            subshell_count = (shell + 1) ** slice_dims - shell ** slice_dims

            index_i = min(remaining // subshell_count, shell)
            if index_i == shell:
                # Once an index on the shell is encountered, the remaining indexes follow
                # multidimensional box unpairing
                return [shell] + multidimensional_box_unpairing([shell + 1] * slice_dims,
                                                                remaining - subshell_count * shell)
            else:
                # Compute the next index the same way by
                # recursing to the next dimension
                return [index_i] + recursive_indexes(dim + 1, remaining - subshell_count * index_i)

    # Subtract out the elements from before this shell and recursively
    # find the index at each dimension from what remains
    return recursive_indexes(0, index - shell ** n)

def floored_root(x, n):
    """Integer component of nth root of x.

    Uses binary search to find the integer y such that y ** n <= x < (y + 1) ** n.

    Adapted from https://stackoverflow.com/a/356206/1485877
    """
    high = 1
    while high ** n <= x:
        high *= 2
    low = high // 2
    while low < high:
        mid = (low + high) // 2
        if low < mid and mid ** n < x:
            low = mid
        elif high > mid and mid ** n > x:
            high = mid
        else:
            return mid
    return mid + 1

def multidimensional_box_unpairing(lengths, index):
    # Should probably assert that index is less than the product of lengths
    n = len(lengths)
    indexes = [0] * n  # Preallocate list
    dimension_product = 1
    # Compute indexes from last to first because that is the order the product is grown
    for dimension in reversed(range(n)):
        indexes[dimension] = index // dimension_product % lengths[dimension]
        dimension_product *= lengths[dimension]

    return indexes
    
test()