WARNING: The text below provides a guidance for a Project Euler problem.

  • If you pride yourself in tackling problems entirely on your own, you are encouraged to exit this post right now.

  • If you are simply looking for a hint for the problem, please visit the »Project Euler’s official forum«. You may find some hints here, but if you don’t want your problem spoiled, scroll cautiously.

  • If you are looking for a readily available copy/paste answer, you will not find it here. All of the code snippets in the post have been stripped of some crucial lines (these are clearly marked). It is my aim that anyone able to replicate the answer based on such snippets will have to understand the solution conceptually. At that point, by virtue of learning something new, I believe he deserves it anyway.

  • Guidance for future problems will not be published before 100 people have solved the problem and there are at least two more recent problems.

Riffle Shuffles

For a positive even number, the »riffle/out shuffle« is defined as the permutation:

The function is the (group-theoretic) order of the permutation and the Project Euler problem tasks us with computing the sum of inverses of , that is:

\[ \text{What is } S(k) := \sum_{n \in s^{-1}(k)} n \ \ ? \]

For start, let’s try and see what actually is. It is not hard to see that we can understand to be the order of 2 in the cyclic group .

Proof. Upon applying , the 0-th and (n-1)-th (first and last) cards remain fixed. It is easy to see that any other card at position will be sent to mod . Therefore, the card number will reach its original position no sooner than after shuffles, where is the smallest number such that mod . But after such number of shuffles, all of the cards will be at their original positions, which gives .

def s(n):
    assert n % 2 == 0
    if n == 2:
        return -1
    curr, cnt = 2, 1
    while curr != 1:
        # <REMOVED>
    return cnt

The value hence satisfies the relation , which is obviously a necessary but not sufficient condition that must satisfy. This allows us to compute and .

def s_inverse(k):
    factors = list(chain(*[[prime] * exp for prime, exp in factorize(2**k - 1)]))
    # Compute all divisors of 2^k - 1. These are candidates d-1 for s(d) = k
    divisors = set([reduce(mul, subset, 1) for subset in powerset(factors)])
    return # <REMOVED>

def S(k):
    return sum(s_inverse(k))

Which turns out fast enough on the problem input with runtime around seconds. Should one not feel good about rewriting ordinary number theoretic functions in pure python as above, the whole solution can be compressed into a one-liner using SymPy.

from sympy.ntheory import n_order, divisors
S = lambda k: sum(x + 1 for x in divisors(2**k - 1) if # <REMOVED>

The most expensive part of our approach is checking whether . It turns out that this can be avoided using a small trick involving the Möbius inversion. Writing where is the »divisor function«, we deduce:

This can be inverted by Möbius, giving us

Using SymPy again we obtain a pretty fast solution.

from sympy.ntheory import mobius, divisor_sigma, divisors
S = lambda k: sum(mobius(m) * (divisor_sigma(2**(k//m) - 1, 0) + divisor_sigma(2**(k//m) - 1, 1)) # <REMOVED>