Problem 622: Riffle Shuffles

Project Euler Problem 622 Solution

EulerSolve provides an optimized solution for Project Euler Problem 622, Riffle Shuffles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For an even deck size $n$, let $s(n)$ be the smallest number of perfect riffle shuffles needed to return the deck to its original order. The shuffle used here is the perfect out-shuffle: the deck is cut into two equal halves and the cards are interleaved so that the original top card stays on top. The task is to compute $$\sum_{s(n)=60} n.$$ The search space is infinite if treated naively, so the solution turns the shuffle into a modular permutation problem and then filters a finite divisor set. Mathematical Approach Number the positions from $0$ at the top to $n-1$ at the bottom. Because $n$ is even, the deck splits into two halves of size $n/2$, and a perfect out-shuffle interleaves those halves in a rigid, deterministic way. Step 1: Model One Shuffle as a Permutation If a card starts in the top half, at position $i$ with $0\le i<n/2$, then after the out-shuffle it moves to position $2i$. If a card starts in the bottom half, at position $i=n/2+j$ with $0\le j<n/2$, then it moves to position $2j+1$. Rewriting that in terms of $i$ gives $$2j+1=2\left(i-\frac{n}{2}\right)+1=2i-(n-1).$$ Therefore, for every position except the bottom card, the shuffle is described by $$f(i)\equiv 2i \pmod{n-1}\qquad (0\le i\le n-2),$$ while the bottom card stays fixed at position $n-1$. The top card is also fixed, since $f(0)=0$....

Detailed mathematical approach

Problem Summary

For an even deck size $n$, let $s(n)$ be the smallest number of perfect riffle shuffles needed to return the deck to its original order. The shuffle used here is the perfect out-shuffle: the deck is cut into two equal halves and the cards are interleaved so that the original top card stays on top.

The task is to compute

$$\sum_{s(n)=60} n.$$

The search space is infinite if treated naively, so the solution turns the shuffle into a modular permutation problem and then filters a finite divisor set.

Mathematical Approach

Number the positions from $0$ at the top to $n-1$ at the bottom. Because $n$ is even, the deck splits into two halves of size $n/2$, and a perfect out-shuffle interleaves those halves in a rigid, deterministic way.

Step 1: Model One Shuffle as a Permutation

If a card starts in the top half, at position $i$ with $0\le i<n/2$, then after the out-shuffle it moves to position $2i$.

If a card starts in the bottom half, at position $i=n/2+j$ with $0\le j<n/2$, then it moves to position $2j+1$. Rewriting that in terms of $i$ gives

$$2j+1=2\left(i-\frac{n}{2}\right)+1=2i-(n-1).$$

Therefore, for every position except the bottom card, the shuffle is described by

$$f(i)\equiv 2i \pmod{n-1}\qquad (0\le i\le n-2),$$

while the bottom card stays fixed at position $n-1$. The top card is also fixed, since $f(0)=0$.

Step 2: Repeated Shuffles Give a Multiplicative Order

Applying the same permutation $k$ times multiplies the position index by $2^k$ modulo $n-1$:

$$f^{(k)}(i)\equiv 2^k i \pmod{n-1}\qquad (0\le i\le n-2).$$

The deck returns to its original order exactly when every position is restored, so we need

$$2^k i\equiv i\pmod{n-1}$$

for all $i$. Checking $i=1$ already forces

$$2^k\equiv 1\pmod{n-1}.$$

Since $n$ is even, $n-1$ is odd, so $\gcd(2,n-1)=1$ and the multiplicative order is defined. Hence

$$\boxed{s(n)=\operatorname{ord}_{n-1}(2).}$$

Step 3: Restrict Candidates Using $2^{60}-1$

Set

$$m=n-1.$$

If $s(n)=60$, then $\operatorname{ord}_m(2)=60$. By the definition of multiplicative order, that implies

$$2^{60}\equiv 1\pmod m,$$

so $m$ must divide $2^{60}-1$:

$$m\mid(2^{60}-1).$$

Therefore every valid deck size has the form $n=d+1$, where $d$ is a divisor of $2^{60}-1$. This is the key reduction: instead of testing all even $n$, we only test the divisors of one fixed integer.

Step 4: Enforce Exact Order Rather Than a Smaller Divisor

Dividing $2^{60}-1$ is only a necessary condition. A divisor $d$ might satisfy $2^{60}\equiv 1\pmod d$ but still have order $1,2,3,4,5,6,10,12,15,20,$ or $30$.

The prime divisors of $60$ are $2$, $3$, and $5$. So $\operatorname{ord}_d(2)=60$ is equivalent to

$$2^{60}\equiv 1\pmod d,$$

$$2^{30}\not\equiv 1\pmod d,\qquad 2^{20}\not\equiv 1\pmod d,\qquad 2^{12}\not\equiv 1\pmod d.$$

This prime-divisor test is sufficient because if the true order were a proper divisor $r$ of $60$, then $60/r$ would have at least one prime factor $p\in\{2,3,5\}$. That would imply $r\mid 60/p$, and then $2^{60/p}\equiv 1\pmod d$, contradicting the test.

Step 5: Worked Example with Shuffle Order $8$

The implementations include a smaller checkpoint before the order-$60$ run, so it is useful to see the same method on $s(n)=8$.

First compute

$$2^8-1=255=3\cdot 5\cdot 17.$$

Its divisors are

$$1,\ 3,\ 5,\ 15,\ 17,\ 51,\ 85,\ 255.$$

Now apply the exact-order test. A valid divisor must satisfy $2^8\equiv 1\pmod d$, but also $2^4\not\equiv 1\pmod d$ and $2^2\not\equiv 1\pmod d$.

The divisors that pass are

$$17,\ 51,\ 85,\ 255.$$

Therefore the corresponding deck sizes are

$$18,\ 52,\ 86,\ 256,$$

and their sum is

$$18+52+86+256=412.$$

This is exactly the checkpoint verified by the implementation.

Step 6: Final Summation Formula

Combining the previous steps, the desired total is

$$\boxed{\sum_{s(n)=60} n=\sum_{\substack{d\mid(2^{60}-1)\\ \operatorname{ord}_d(2)=60}}(d+1).}$$

So the entire problem reduces to three finite tasks: factor $2^{60}-1$, enumerate its divisors, and keep only those divisors for which $2$ has exact order $60$.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. They first factor $2^{60}-1$ using fast modular arithmetic, a deterministic primality test suitable for 64-bit integers, and Pollard's rho splitting for composite factors. Once the prime factorization is known, they recursively generate every divisor.

For each divisor, the implementation performs fast modular exponentiation to test the order conditions above. Any divisor that satisfies the exact-order filter contributes $d+1$ to the running sum. Before computing the target case $60$, the code also checks the smaller order-$8$ example whose total is $412$, providing a direct sanity check that the divisor generation and order filter agree.

The three languages differ only in arithmetic details. The C++ version uses native 128-bit multiplication for safe modular products, Python uses arbitrary-precision integers, and Java uses an overflow-safe modular multiplication routine. The underlying algorithm is otherwise identical.

Complexity Analysis

For this specific problem, the target order $60$ is fixed, so the total running time is effectively constant. More structurally, the expensive step is factoring $2^{60}-1$; after that, the algorithm enumerates all divisors and performs a constant number of modular exponentiations for each one.

If $\tau(2^{60}-1)$ denotes the number of divisors of $2^{60}-1$, then the filtering phase costs $O(\tau(2^{60}-1)\log 60)$ modular multiplications after factorization, because each candidate is tested with exponents $60$, $30$, $20$, and $12$. Memory usage is $O(\tau(2^{60}-1))$ if the full divisor list is stored, and can be viewed as linear in the number of generated divisors.

Footnotes and References

Problem page: https://projecteuler.net/problem=622
Faro shuffle: Wikipedia — Faro shuffle
Multiplicative order: Wikipedia — Multiplicative order
Modular arithmetic: Wikipedia — Modular arithmetic
Pollard's rho algorithm: Wikipedia — Pollard's rho algorithm

Problem 622 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <map>
#include <numeric>
#include <random>
#include <string>
#include <utility>
#include <vector>

using u64 = unsigned long long;
using u128 = __uint128_t;

static u64 mod_mul(u64 a, u64 b, u64 mod) { return (__uint128_t)a * b % mod; }

static u64 mod_pow(u64 a, u64 e, u64 mod) {
    if (mod == 1) return 0;
    u64 r = 1 % mod;
    while (e) {
        if (e & 1) r = mod_mul(r, a, mod);
        a = mod_mul(a, a, mod);
        e >>= 1;
    }
    return r;
}

static bool is_prime(u64 n) {
    if (n < 2) return false;
    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (n % p == 0) return n == p;
    }

    u64 d = n - 1;
    int s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        ++s;
    }

    auto witness = [&](u64 a) -> bool {
        if (a % n == 0) return false;
        u64 x = mod_pow(a, d, n);
        if (x == 1 || x == n - 1) return false;
        for (int i = 1; i < s; ++i) {
            x = mod_mul(x, x, n);
            if (x == n - 1) return false;
        }
        return true;
    };

    for (u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
        if (witness(a)) return false;
    }
    return true;
}

static u64 pollard_rho(u64 n) {
    if ((n & 1ULL) == 0) return 2;
    static std::mt19937_64 rng(987654321);
    std::uniform_int_distribution<u64> dist(0, n - 1);

    while (true) {
        u64 c = dist(rng) % n;
        u64 x = dist(rng) % n;
        u64 y = x;
        u64 d = 1;

        auto f = [&](u64 v) { return (mod_mul(v, v, n) + c) % n; };

        while (d == 1) {
            x = f(x);
            y = f(f(y));
            u64 diff = (x > y) ? (x - y) : (y - x);
            d = std::gcd(diff, n);
        }
        if (d != n) return d;
    }
}

static void factor_rec(u64 n, std::map<u64, int> &out) {
    if (n == 1) return;
    if (is_prime(n)) {
        out[n]++;
        return;
    }
    u64 d = pollard_rho(n);
    factor_rec(d, out);
    factor_rec(n / d, out);
}

static void gen_divisors(const std::vector<std::pair<u64, int>> &pf, int idx, u64 cur, std::vector<u64> &out) {
    if (idx == (int)pf.size()) {
        out.push_back(cur);
        return;
    }
    auto [p, e] = pf[idx];
    u64 v = 1;
    for (int i = 0; i <= e; ++i) {
        gen_divisors(pf, idx + 1, cur * v, out);
        v *= p;
    }
}

static std::vector<int> prime_factors_distinct(int n) {
    std::vector<int> pf;
    for (int p = 2; 1LL * p * p <= n; ++p) {
        if (n % p != 0) continue;
        pf.push_back(p);
        while (n % p == 0) n /= p;
    }
    if (n > 1) pf.push_back(n);
    return pf;
}

static bool has_exact_order(u64 mod, int ord, const std::vector<int> &ord_primes) {
    if (mod == 1) return false;
    if (mod_pow(2, (u64)ord, mod) != 1) return false;
    for (int p : ord_primes) {
        if (mod_pow(2, (u64)(ord / p), mod) == 1) return false;
    }
    return true;
}

static u128 sum_n_with_shuffle_order(int ord) {
    const u64 M = (ord == 64) ? ~0ULL : ((1ULL << ord) - 1ULL);
    std::map<u64, int> f;
    factor_rec(M, f);

    std::vector<std::pair<u64, int>> pf;
    pf.reserve(f.size());
    for (auto [p, e] : f) pf.push_back({p, e});

    std::vector<u64> divs;
    divs.reserve(8192);
    gen_divisors(pf, 0, 1, divs);

    const std::vector<int> ord_primes = prime_factors_distinct(ord);
    u128 sum = 0;
    for (u64 d : divs) {
        if (has_exact_order(d, ord, ord_primes)) sum += (u128)(d + 1);
    }
    return sum;
}

static std::string to_string_u128(u128 x) {
    if (x == 0) return "0";
    std::string s;
    while (x) {
        int digit = (int)(x % 10);
        s.push_back(char('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

int main() {
    assert(sum_n_with_shuffle_order(8) == 412);
    std::cout << to_string_u128(sum_n_with_shuffle_order(60)) << "\n";
    return 0;
}

Python

import math
import random
from collections import defaultdict

def solve():
    def mod_mul(a, b, m): return a * b % m
    def mod_pow(a, e, m):
        r = 1 % m
        a %= m
        while e > 0:
            if e & 1: r = r * a % m
            a = a * a % m
            e >>= 1
        return r

    def is_prime(n):
        if n < 2: return False
        for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
            if n % p == 0: return n == p
        d = n - 1; s = 0
        while d % 2 == 0: d >>= 1; s += 1
        for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
            if a % n == 0: continue
            x = mod_pow(a, d, n)
            if x == 1 or x == n - 1: continue
            w = True
            for _ in range(1, s):
                x = x * x % n
                if x == n - 1: w = False; break
            if w: return False
        return True

    def pollard_rho(n):
        if n % 2 == 0: return 2
        while True:
            c = random.randint(1, n - 1)
            x = random.randint(0, n - 1)
            y = x; d = 1
            while d == 1:
                x = (x * x + c) % n
                y = (y * y + c) % n
                y = (y * y + c) % n
                d = math.gcd(abs(x - y), n)
            if d != n: return d

    def factor(n):
        if n <= 1: return {}
        if is_prime(n): return {n: 1}
        d = pollard_rho(n)
        f1 = factor(d)
        f2 = factor(n // d)
        for p, e in f2.items():
            f1[p] = f1.get(p, 0) + e
        return f1

    def gen_divs(pf, idx, cur, out):
        if idx == len(pf):
            out.append(cur)
            return
        p, e = pf[idx]
        v = 1
        for _ in range(e + 1):
            gen_divs(pf, idx + 1, cur * v, out)
            v *= p

    def distinct_pf(n):
        pf = []
        p = 2
        while p * p <= n:
            if n % p == 0:
                pf.append(p)
                while n % p == 0: n //= p
            p += 1
        if n > 1: pf.append(n)
        return pf

    def has_exact_order(mod, ord, ord_primes):
        if mod == 1: return False
        if mod_pow(2, ord, mod) != 1: return False
        for p in ord_primes:
            if mod_pow(2, ord // p, mod) == 1: return False
        return True

    random.seed(987654321)
    ord = 60
    M = (1 << ord) - 1
    f = factor(M)
    pf = list(f.items())
    divs = []
    gen_divs(pf, 0, 1, divs)
    ord_primes = distinct_pf(ord)

    total = sum(d + 1 for d in divs if has_exact_order(d, ord, ord_primes))
    return str(total)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;
import java.util.Map;
import java.util.Random;
import java.util.TreeMap;

public class Euler622 {
    static long modMul(long a, long b, long mod) {
        long res = 0;
        a %= mod;
        while (b > 0) {
            if ((b & 1) == 1) {
                res = (res + a) % mod;
            }
            a = (a << 1) % mod;
            b >>= 1;
        }
        return res;
    }

    static long modPow(long a, long e, long mod) {
        if (mod == 1)
            return 0;
        long r = 1 % mod;
        while (e > 0) {
            if ((e & 1) == 1)
                r = modMul(r, a, mod);
            a = modMul(a, a, mod);
            e >>= 1;
        }
        return r;
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n % p == 0)
                return n == p;
        }

        long d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }

        long[] bases = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
        for (long a : bases) {
            if (a % n == 0)
                continue;
            long x = modPow(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean composite = true;
            for (int i = 1; i < s; i++) {
                x = modMul(x, x, n);
                if (x == n - 1) {
                    composite = false;
                    break;
                }
            }
            if (composite)
                return false;
        }
        return true;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static Random rng = new Random(987654321L);

    static long pollardRho(long n) {
        if ((n & 1) == 0)
            return 2;

        while (true) {
            long c = (long) (rng.nextDouble() * n) % n;
            long x = (long) (rng.nextDouble() * n) % n;
            long y = x;
            long d = 1;

            while (d == 1) {
                x = (modMul(x, x, n) + c) % n;
                y = (modMul(y, y, n) + c) % n;
                y = (modMul(y, y, n) + c) % n;
                long diff = Math.abs(x - y);
                d = gcd(diff, n);
            }
            if (d != n)
                return d;
        }
    }

    static void factorRec(long n, Map<Long, Integer> out) {
        if (n == 1)
            return;
        if (isPrime(n)) {
            out.put(n, out.getOrDefault(n, 0) + 1);
            return;
        }
        long d = pollardRho(n);
        factorRec(d, out);
        factorRec(n / d, out);
    }

    static void genDivisors(List<long[]> pf, int idx, long cur, List<Long> out) {
        if (idx == pf.size()) {
            out.add(cur);
            return;
        }
        long p = pf.get(idx)[0];
        long e = pf.get(idx)[1];
        long v = 1;
        for (int i = 0; i <= e; i++) {
            genDivisors(pf, idx + 1, cur * v, out);
            v *= p;
        }
    }

    static List<Integer> primeFactorsDistinct(int n) {
        List<Integer> pf = new ArrayList<>();
        for (int p = 2; p * p <= n; p++) {
            if (n % p == 0) {
                pf.add(p);
                while (n % p == 0)
                    n /= p;
            }
        }
        if (n > 1)
            pf.add(n);
        return pf;
    }

    static boolean hasExactOrder(long mod, int ord, List<Integer> ordPrimes) {
        if (mod == 1)
            return false;
        if (modPow(2, ord, mod) != 1)
            return false;
        for (int p : ordPrimes) {
            if (modPow(2, ord / p, mod) == 1)
                return false;
        }
        return true;
    }

    static String sumNWithShuffleOrder(int ord) {
        long M = (1L << ord) - 1L;

        Map<Long, Integer> f = new TreeMap<>();
        factorRec(M, f);

        List<long[]> pf = new ArrayList<>();
        for (Map.Entry<Long, Integer> entry : f.entrySet()) {
            pf.add(new long[] { entry.getKey(), entry.getValue() });
        }

        List<Long> divs = new ArrayList<>();
        genDivisors(pf, 0, 1, divs);

        List<Integer> ordPrimes = primeFactorsDistinct(ord);
        long sum = 0;
        for (long d : divs) {
            if (hasExactOrder(d, ord, ordPrimes)) {
                sum += (d + 1);
            }
        }
        return Long.toString(sum);
    }

    public static String solve() {
        return sumNWithShuffleOrder(60);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}