Counting integer tuples with specific property

Question

Given two numbers M_max,M_min count all the r-tuples of prime numbers (M_1,...,M_r) such that:

1. M_min < M_i < M_max
2. every M_i-1 is divisible by at least one square of an odd prime
3. lcm(M_1-1,...,M_r-1)=(M_1-1)*(M_2-1)...(M_r-1)/2^(r-1)

Rephrased Question: How many r-tuples among the k numbers M_1,M_2,...,M_k are coprime?

Answer 1

Here is working code. In fairly unoptimized Python.

# Stupid simple Euclidean primes filter.
def primes ():
    found_primes = []
    factor_filter = [False, False]
    block_size = 2

    candidate = 2
    while True:
        factor_filter.extend([False for _ in factor_filter])
        for p in found_primes:
            if 2*block_size <= p*p:
                break
            i = ((block_size + p - 1) // p) * p
            while i < 2 * block_size:
                factor_filter[i] = 1
                i = i + p

        block_size += block_size

        while (candidate < block_size):
            if 0 == factor_filter[candidate]:
                yield candidate
                found_primes.append(candidate)
            candidate = candidate + 1

def factor_range (m, n):
    factored = []
    remainder = []
    for i in range(n - m + 1):
        factored.append({})
        remainder.append(i+m)

    for p in primes():
        if n < p*p:
            break
        else:
            count = 1
            q = p
            while q <= n:
                # Find the first thing divisible by q
                i = ((m + q - 1) // q) * q
                while i <= n:
                    factored[i-m][p] = count
                    remainder[i-m] = remainder[i-m] // p
                    i += q
                q = p*q
                count += 1

    for i in range(len(remainder)):
        if 1 < remainder[i]:
            factored[i][remainder[i]] = 1

    return factored

def int_tuples (m, n, j):
    factored = factor_range(m, n)

    available = []
    for i in range(len(factored)):
        if 0 < i and i + m in factored[i]:
            for p, count in factored[i-1].items():
                if 2 < p and 1 < count:
                    available.append(i+m)
                    break

    prime_divides = {}
    for i in available:
        if m < i:
            for p, count in factored[i-1-m].items():
                if 1 < count or p != 2:
                    if p not in prime_divides:
                        prime_divides[p] = set([i])
                    else:
                        prime_divides[p].add(i)

    options_stack = [(-1, available, False)]
    chosen = []
    while 0 < len(options_stack):
        #print(("...", chosen, options_stack))
        i, options, is_pop = options_stack.pop()
        i = i+1
        if len(chosen) == j:
            yield tuple(chosen)
            if is_pop:
                chosen.pop()
        elif len(chosen) + len(options) - i < j:
            # Cannot generate an answer.
            if is_pop:
                chosen.pop()
        elif len(options) <= i:
            # i reached the end.
            if is_pop:
                chosen.pop()
        else:
            options_stack.append((i, options, is_pop))
            chosen.append(options[i])
            options = options[i:]
            for p, count in factored[options[0]-1].items():
                if p == 2 and 1 == count:
                    continue
                if p in prime_divides:
                    options = [option for option in options if option not in prime_divides[p]]
            options_stack.append((-1, options, True))

The idea is to find tuples such that the (M-1)/2 reduced series are pairwise relatively prime to each other. Which I'm doing by what amounts to a recursive search. However I kept explicit stacks so that I could make a generator by simply doing a yield.

My first try involved keeping bitvectors for efficiency. But a heavier weight data structure followed by determining "don't search farther" conclusions earlier works out better.

---

I decided to try recursive after all, and I'm glad I did. The use of caching is yet another step towards having a heavyweight data structure to avoid repeated calculations.

# Stupid simple Euclidean primes filter.
def primes ():
    found_primes = []
    factor_filter = [False, False]
    block_size = 2

    candidate = 2
    while True:
        # We get here when we need to double the length of our sieve.
        # First add new elements.
        factor_filter.extend([False for _ in factor_filter])
        # Cross off the ones that divisible by known primes.
        for p in found_primes:
            if 2*block_size <= p*p:
                break
            i = ((block_size + p - 1) // p) * p
            while i < 2 * block_size:
                factor_filter[i] = True
                i = i + p

        # Record that we did so.
        block_size += block_size

        # Now that we have sieved, start returning primes.
        while (candidate < block_size):
            if 0 == factor_filter[candidate]:
                yield candidate
                found_primes.append(candidate)
            candidate = candidate + 1

# Takes a range of numbers, factors all of them at once.
def factor_range (m, n):
    # This will hold the factorization.
    factored = []
    # This will hold the remainder we have not factorized.
    remainder = []
    # We start with no factors and the number as remainder.
    for i in range(n - m + 1):
        factored.append({})
        remainder.append(i+m)

    # We run through the primes up to sqrt(n)    
    for p in primes():
        if n < p*p:
            break
        else:
            # For each power of p, we add it to the factorization and divide it form the remainder.
            count = 1
            q = p
            while q <= n:
                # Find the first thing divisible by q
                i = ((m + q - 1) // q) * q
                while i <= n:
                    factored[i-m][p] = count
                    remainder[i-m] = remainder[i-m] // p
                    i += q
                q = p*q
                count += 1

    # Any left-over remainders must be primes.    
    for i in range(len(remainder)):
        if 1 < remainder[i]:
            factored[i][remainder[i]] = 1

    return factored

# This function does all of the work.
def _int_tuples (m, n, j):
    # BUG FIX HERE:
    # If our minimum is odd, extend the boundary down one.
    # This does not introduce new solutions, but does mean we
    # factorize all of the numbers that we need.
    if 1 == m // 2:
        m = m-1

    # Factorize all numbers.
    factored = factor_range(m, n)

    # available will be all of our candidate primes.
    available = []
    for i in range(len(factored)):
        if 0 < i and i + m in factored[i]:
            for p, count in factored[i-1].items():
                if 2 < p and 1 < count:
                    available.append(i+m)
                    break

    # Construct a data structure to answer, "which other candidate primes
    # have p-1 share a non-trivial factor with this one"
    prime_divides = {}
    for i in available:
        if m < i:
            for p, count in factored[i-1-m].items():
                if 1 < count or p != 2:
                    if p not in prime_divides:
                        prime_divides[p] = set([i])
                    else:
                        prime_divides[p].add(i)

    # PERFORMANCE HACK: We will recurse over lists of remaining candidates.
    # Putting candidates that eliminate lots of other candidates first
    # makes those lists smaller and that recursion more efficient.  We
    # also improve the odds that we will be looking at a list we have
    # already seen and already have an answer for.  And, finally, the
    # candidates at the end of our list will largely not filter each
    # other out.  So we are spending time on things that are more likely to
    # give us more solutions.
    count_related = {}
    for i in available:
        total = 0
        for p, count in factored[i-1-m].items():
            if 1 < count or p != 2:
                total += len(prime_divides[p]) - 1
        count_related[i] = total
    available = sorted(available, key = lambda i: -count_related[i])

    # To spend less time filtering, forget about the primes that only
    # divide one thing.
    for p in list(prime_divides.keys()):
        if 1 == len(prime_divides[p]):
            i = list(prime_divides[p])[0] # Get the element
            prime_divides.pop(p)
            factored[i-m-1].pop(p)

    # We will keep a cache of already computed results.
    # That way we do not have to recalculate them.  This trick
    # is called "memoizing".
    cache = {}

    # Recursive search here.  options are remaining candidate primes
    # and needed are how many more of them we need.
    # It returns (count_solutions, data_structure_encoding_solutions)
    def search (options, needed):
        # This is the key we will cache results under.
        key = (tuple(options), needed)
        # Return trivial base cases.
        if 0 == needed:
            return [1, None]
        elif len(options) < needed:
            return [0, None]
        elif key not in cache:
            # It is not trivial, and we have never calculated it.
            # Actually do work.
            answer = []
            total = 0
            # For each first candidate we can choose.
            for i in range(len(options) - needed + 1):
                # Remaining candidates are farther along.
                next_options = options[i+1:]
                # Filter out ones whose p-1 shares a factor with this one.
                for p, count in factored[options[i] - m-1].items():
                    if p == 2 and 1 == count:
                        # All are even, don't filter on p=2 unless divisible by 4.
                        continue
                    # Use list comprehension to do the filtering.
                    next_options = [option for option in next_options if option not in prime_divides[p]]
                # Now the recursive call.
                count, result = search(next_options, needed - 1)
                # Did we find any?
                if 0 < count:
                    total += count
                    answer.append([options[i], result])
            # cache it
            cache[key] = [total, answer]
        # and return the cached answer
        return cache[key]

    # And actually do the recursive search.    
    return search(available, j)

# The count is the first thing returned.
def count_int_tuples (m, n, j):
    return _int_tuples(m, n, j)[0]

# This will yield all of the tuples
def int_tuples (m, n, j):
    # the encoded data structure is the second thing returned.
    # It has the structure:
    # [
    #    [p1, recursive_encoding of the rest],
    #    [p2, recursive_encoding of the rest],
    #    ...
    # ]
    #
    data = _int_tuples(m, n, j)[1]
    # chosen is our chosen primes.  Since processing always pops first,
    # before sticking the next candidate in, it needs to start with
    # something, anything.  That will be thrown away.
    chosen = [None]
    # Think of stack as a recursive call stack.  Each call has a
    # (position, data).
    # We start before the beginning of all of our data.
    stack = [(-1, data)]
    # While there is work to do...
    while 0 < len(stack):
        # get the next call
        i, data = stack.pop()
        # clear the previous choice
        chosen.pop()
        # advance
        i = i+1

        # If we have another (candidate, more_data) pair to process
        if i < len(data):
            # Choose our candidate
            chosen.append(data[i][0])
            # Add walking through current data to the call stack.
            stack.append([i, data])
            # Did we get a full solution?
            if j == len(chosen):
                # Return it here, then continue.
                # If that is confusing, look up Python iterators.
                yield chosen
            else:
                # Add a recursive call to this data structure.
                # Because stacks are Last In First Out (LIFO), this
                # will be processed before undoing our choice.
                stack.append([-1, data[i][1]])
                chosen.append(None)

If I run print(count_int_tuples(2, 40000, 3)) on my laptop, it takes just over a second to decide that there are 2,198,808 tuples of length 3 meeting your condition.