Problem 312: Cyclic Paths on Sierpiński Graphs

Project Euler Problem 312 Solution

EulerSolve provides an optimized solution for Project Euler Problem 312, Cyclic Paths on Sierpiński Graphs, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $C(n)$ be the number of Hamiltonian cycles on the Sierpinski graph $S_n$ described in the problem. We are asked to compute $$C(C(C(10000))) \pmod{13^8}.$$ The direct values explode immediately, so the only viable route is: $$\text{graph recurrence} \Longrightarrow \text{closed form} \Longrightarrow \text{modular exponent reduction}.$$ Mathematical Approach 1) Decompose $S_{n+1}$ into three copies of $S_n$. The graph $S_{n+1}$ consists of three corner-sharing copies of $S_n$. A Hamiltonian cycle on $S_{n+1}$ cannot stay closed inside one copy; it must enter each copy through one shared corner and leave through the other shared corner. Therefore each copy contributes a Hamiltonian path between two corner vertices . Let $P(n)$ denote the number of such corner-to-corner Hamiltonian paths in $S_n$ for a fixed ordered copy position. Then the three copies are independent once their corner roles are fixed, so $$C(n+1)=P(n)^3.$$ 2) Why $P(n)=3C(n)$. Inside one copy of $S_n$, a Hamiltonian path is obtained from a Hamiltonian cycle by choosing which one of the three outer corners plays the role of the “unused outer corner” when the copy is embedded into the next level. There are exactly three symmetric choices, and each choice converts bijectively between the cycle state and the required corner-to-corner path state....

Detailed mathematical approach

Problem Summary

Let $C(n)$ be the number of Hamiltonian cycles on the Sierpinski graph $S_n$ described in the problem. We are asked to compute

$$C(C(C(10000))) \pmod{13^8}.$$

The direct values explode immediately, so the only viable route is:

$$\text{graph recurrence} \Longrightarrow \text{closed form} \Longrightarrow \text{modular exponent reduction}.$$

Mathematical Approach

1) Decompose $S_{n+1}$ into three copies of $S_n$.

The graph $S_{n+1}$ consists of three corner-sharing copies of $S_n$. A Hamiltonian cycle on $S_{n+1}$ cannot stay closed inside one copy; it must enter each copy through one shared corner and leave through the other shared corner. Therefore each copy contributes a Hamiltonian path between two corner vertices.

Let $P(n)$ denote the number of such corner-to-corner Hamiltonian paths in $S_n$ for a fixed ordered copy position. Then the three copies are independent once their corner roles are fixed, so

$$C(n+1)=P(n)^3.$$

2) Why $P(n)=3C(n)$.

Inside one copy of $S_n$, a Hamiltonian path is obtained from a Hamiltonian cycle by choosing which one of the three outer corners plays the role of the “unused outer corner” when the copy is embedded into the next level. There are exactly three symmetric choices, and each choice converts bijectively between the cycle state and the required corner-to-corner path state. Hence

$$P(n)=3C(n).$$

Substituting into the previous relation gives the recurrence used by the solver:

$$C(n+1)=(3C(n))^3,\qquad C(3)=8.$$

As a quick check,

$$C(4)=(3\cdot 8)^3=24^3=13824,$$

which matches the iterative checkpoint in the code.

3) Closed form.

Write

$$C(n)=8\cdot 12^{e_n}\qquad (n\ge 3).$$

Then

$$C(n+1)=(3C(n))^3=(24\cdot 12^{e_n})^3=8\cdot 12^{3e_n+3},$$

so the exponents satisfy

$$e_{n+1}=3e_n+3,\qquad e_3=0.$$

This linear recurrence solves to

$$e_n=\frac{3^{n-2}-3}{2}.$$

Therefore

$$C(n)=8\cdot 12^{(3^{n-2}-3)/2}.$$

For example, for $n=5$,

$$e_5=\frac{3^3-3}{2}=12,\qquad C(5)=8\cdot12^{12}=71328803586048,$$

again exactly the checkpoint used by the implementation.

4) Reduce the outer power modulo $13^k$.

To compute $C(n)\bmod 13^k$, we need the multiplicative order of $12$ modulo $13^k$. Since

$$12\equiv -1 \pmod{13},$$

its order modulo $13$ is $2$. Also

$$12^2-1=143=11\cdot 13,$$

so only one factor of $13$ divides $12^2-1$. The standard lifting rule for orders over odd prime powers then gives

$$\operatorname{ord}_{13^k}(12)=2\cdot 13^{k-1}.$$

Hence the exponent only matters modulo

$$2\cdot 13^{k-1}.$$

5) Why the code computes $3^{n-2}$ modulo $4\cdot 13^{k-1}$.

The exponent is

$$e_n=\frac{3^{n-2}-3}{2}.$$

We must divide by $2$ after reducing modulo the exponent period. To do this safely, compute

$$3^{n-2}\pmod{4\cdot 13^{k-1}},$$

then subtract $3$, which is always even, and only then divide by $2$. This yields

$$e_n \pmod{2\cdot 13^{k-1}}.$$

6) How much of $n$ do we really need?

Now the inner task is to compute $3^{n-2}$ modulo $4\cdot 13^{k-1}$. Since $\gcd(3,4\cdot 13^{k-1})=1$, the exponent $n-2$ may be reduced modulo the Carmichael period of that modulus.

For $k\ge 2$,

$$\lambda(4\cdot 13^{k-1})=\operatorname{lcm}(\lambda(4),\lambda(13^{k-1}))=\operatorname{lcm}(2,12\cdot 13^{k-2})=12\cdot 13^{k-2}.$$

So to compute $C(n)\bmod 13^k$, it is enough to know

$$n \pmod{12\cdot 13^{k-2}},$$

plus whether $n\le 2$. For the huge nested values here, $n$ is certainly large, so only the residue matters.

7) Why CRT appears.

For every $n\ge 4$, the recurrence shows that $C(n)$ is divisible by $12$, because

$$C(n+1)=(3C(n))^3$$

has an obvious factor $3^3$, and from $C(4)=24^3$ onward the factor $4$ is also permanent. Therefore when the next level needs $n\bmod (12\cdot 13^t)$, we already know

$$n\equiv 0 \pmod{12}.$$

The solver separately computes

$$n\equiv r \pmod{13^t}$$

and combines the two congruences by the Chinese remainder theorem:

$$n\equiv 0 \pmod{12},\qquad n\equiv r \pmod{13^t}\quad \Longrightarrow \quad n\bmod (12\cdot 13^t).$$

8) Nested evaluation chain.

Define

$$a=C(10000),\qquad b=C(a),\qquad c=C(b).$$

To compute $b \bmod 13^6$, we only need

$$a \bmod (12\cdot 13^4).$$

So the code first computes

$$a \bmod 13^4,$$

then reconstructs $a \bmod (12\cdot 13^4)$ using CRT and $a\equiv 0\pmod{12}$.

Similarly, to compute $c \bmod 13^8$, we only need

$$b \bmod (12\cdot 13^6).$$

So the same trick is applied one more time. This completely avoids ever forming the astronomical integers $a$ and $b$.

Worked Checks

The implementation verifies the following checkpoints:

$$C(1)=1,\qquad C(2)=1,\qquad C(3)=8,$$

$$C(4)=13824,\qquad C(5)=71328803586048,$$

$$C(10000)\bmod 10^8 = 37652224,$$

$$C(10000)\bmod 13^8 = 617720485.$$

The final answer is

$$C(C(C(10000)))\bmod 13^8 = 324681947.$$

Algorithm

1) Implement fast modular multiplication and binary exponentiation.

2) Implement the closed-form evaluator

$$C(n)\bmod 13^k=8\cdot 12^{((3^{n-2}-3)/2)\bmod (2\cdot 13^{k-1})}\pmod{13^k},$$

using the $4\cdot 13^{k-1}$ trick before dividing by $2$.

3) Use CRT to lift residues from $13^t$ to $12\cdot 13^t$.

4) Evaluate the three nested levels with the minimal modulus needed by the next stage.

Complexity Analysis

Everything reduces to a small number of modular exponentiations. Each one costs $O(\log M)$ modular multiplications for modulus $M$. So the total runtime is tiny; the difficulty is entirely mathematical, not computational.

Problem 312 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % static_cast<u128>(mod));
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }
    base %= mod;
    u64 result = 1ULL % mod;
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

u64 pow_u64(u64 base, int exp) {
    u64 result = 1ULL;
    while (exp-- > 0) {
        result *= base;
    }
    return result;
}

u64 inverse_mod(u64 a, const u64 mod) {
    a %= mod;
    // mod is prime-power of 13 in this problem, so a^(phi-1) gives inverse for gcd(a,mod)=1.
    const u64 phi = mod - mod / 13ULL;
    return pow_mod(a, phi - 1ULL, mod);
}

u64 crt_coprime(const u64 a1, const u64 m1, const u64 a2, const u64 m2) {
    const u64 inv_m1 = inverse_mod(m1 % m2, m2);
    const u64 delta = (a2 >= a1 % m2) ? (a2 - (a1 % m2)) : (a2 + m2 - (a1 % m2));
    const u64 t = mul_mod(delta, inv_m1, m2);
    return a1 + m1 * t;
}

u64 c_mod_iterative(const u64 n, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }
    if (n <= 2ULL) {
        return 1ULL % mod;
    }

    u64 c = 8ULL % mod;  // C(3)
    for (u64 i = 4ULL; i <= n; ++i) {
        const u64 v = (3ULL * c) % mod;
        c = mul_mod(mul_mod(v, v, mod), v, mod);
    }
    return c;
}

u64 c_mod_13_power_from_n_residue(const u64 n_residue, const bool n_is_large, const int k) {
    const u64 mod = pow_u64(13ULL, k);

    if (!n_is_large && n_residue <= 2ULL) {
        return 1ULL % mod;
    }

    const u64 ord12 = 2ULL * pow_u64(13ULL, k - 1);        // ord_{13^k}(12)
    const u64 two_ord = 2ULL * ord12;                      // 4 * 13^{k-1}
    const u64 lambda = (k >= 2) ? (12ULL * pow_u64(13ULL, k - 2)) : 2ULL;

    const u64 exp = (n_residue + lambda - 2ULL) % lambda;  // n - 2 mod lambda(two_ord)
    const u64 x = pow_mod(3ULL, exp, two_ord);              // x = 3^{n-2} mod (2*ord)

    const u64 numer = (x + two_ord - 3ULL) % two_ord;       // x - 3 is even in this modulus
    const u64 e = (numer / 2ULL) % ord12;                   // ((3^{n-2}-3)/2) mod ord

    return mul_mod(8ULL % mod, pow_mod(12ULL, e, mod), mod);
}

u64 solve() {
    const u64 n0 = 10000ULL;

    // a = C(10000). For n >= 4, C(n) is divisible by 12.
    const u64 a_mod_13_4 = c_mod_13_power_from_n_residue(n0, false, 4);
    const u64 a_mod_12x13_4 = crt_coprime(0ULL, 12ULL, a_mod_13_4, pow_u64(13ULL, 4));

    // b = C(a). Need b mod (12 * 13^6) to evaluate C(b) mod 13^8.
    const u64 b_mod_13_6 = c_mod_13_power_from_n_residue(a_mod_12x13_4, true, 6);
    const u64 b_mod_12x13_6 = crt_coprime(0ULL, 12ULL, b_mod_13_6, pow_u64(13ULL, 6));

    // c = C(b) mod 13^8.
    return c_mod_13_power_from_n_residue(b_mod_12x13_6, true, 8);
}

bool run_checkpoints() {
    // Given checks from statement.
    if (c_mod_iterative(1ULL, 1000000000000000ULL) != 1ULL) {
        std::cerr << "Checkpoint failed: C(1)\n";
        return false;
    }
    if (c_mod_iterative(2ULL, 1000000000000000ULL) != 1ULL) {
        std::cerr << "Checkpoint failed: C(2)\n";
        return false;
    }
    if (c_mod_iterative(5ULL, 1000000000000000ULL) != 71328803586048ULL) {
        std::cerr << "Checkpoint failed: C(5)\n";
        return false;
    }
    if (c_mod_iterative(10000ULL, 100000000ULL) != 37652224ULL) {
        std::cerr << "Checkpoint failed: C(10000) mod 1e8\n";
        return false;
    }
    if (c_mod_iterative(10000ULL, pow_u64(13ULL, 8)) != 617720485ULL) {
        std::cerr << "Checkpoint failed: C(10000) mod 13^8\n";
        return false;
    }

    // Internal consistency: closed-form modulo 13^8 matches iterative value for n=10000.
    if (c_mod_13_power_from_n_residue(10000ULL, false, 8) != 617720485ULL) {
        std::cerr << "Checkpoint failed: closed-form/iterative mismatch\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve() << '\n';
    return 0;
}

Python

def inverse_mod(a, mod):
    a %= mod
    phi = mod - mod // 13
    return pow(a, phi - 1, mod)

def crt_coprime(a1, m1, a2, m2):
    inv_m1 = inverse_mod(m1 % m2, m2)
    delta = a2 - (a1 % m2)
    t = (delta * inv_m1) % m2
    return a1 + m1 * t

def c_mod_13_power_from_n_residue(n_residue, n_is_large, k):
    mod = 13 ** k
    
    if not n_is_large and n_residue <= 2:
        return 1 % mod
        
    ord12 = 2 * (13 ** (k - 1))
    two_ord = 2 * ord12
    lambda_val = 12 * (13 ** (k - 2)) if k >= 2 else 2
    
    exp = (n_residue + lambda_val - 2) % lambda_val
    x = pow(3, exp, two_ord)
    
    numer = (x + two_ord - 3) % two_ord
    e = (numer // 2) % ord12
    
    return ((8 % mod) * pow(12, e, mod)) % mod

def solve():
    n0 = 10000
    
    a_mod_13_4 = c_mod_13_power_from_n_residue(n0, False, 4)
    a_mod_12x13_4 = crt_coprime(0, 12, a_mod_13_4, 13 ** 4)
    
    b_mod_13_6 = c_mod_13_power_from_n_residue(a_mod_12x13_4, True, 6)
    b_mod_12x13_6 = crt_coprime(0, 12, b_mod_13_6, 13 ** 6)
    
    result = c_mod_13_power_from_n_residue(b_mod_12x13_6, True, 8)
    return str(result)

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

Java

import java.math.BigInteger;

public class Euler312 {
    static long pow_u64(long base, int exp) {
        long result = 1;
        while (exp-- > 0)
            result *= base;
        return result;
    }

    static long mul_mod(long a, long b, long mod) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long pow_mod(long base, long exp, long mod) {
        if (mod == 1)
            return 0;
        return BigInteger.valueOf(base).modPow(BigInteger.valueOf(exp), BigInteger.valueOf(mod)).longValue();
    }

    static long inverse_mod(long a, long mod) {
        a %= mod;
        long phi = mod - mod / 13;
        return pow_mod(a, phi - 1, mod);
    }

    static long crt_coprime(long a1, long m1, long a2, long m2) {
        long inv_m1 = inverse_mod(m1 % m2, m2);
        long delta = (a2 >= a1 % m2) ? (a2 - (a1 % m2)) : (a2 + m2 - (a1 % m2));
        long t = mul_mod(delta, inv_m1, m2);
        return a1 + m1 * t;
    }

    static long c_mod_13_power_from_n_residue(long n_residue, boolean n_is_large, int k) {
        long mod = pow_u64(13, k);

        if (!n_is_large && n_residue <= 2) {
            return 1 % mod;
        }

        long ord12 = 2 * pow_u64(13, k - 1);
        long two_ord = 2 * ord12;
        long lambda_val = (k >= 2) ? (12 * pow_u64(13, k - 2)) : 2;

        long exp = (n_residue + lambda_val - 2) % lambda_val;
        long x = pow_mod(3, exp, two_ord);

        long numer = (x + two_ord - 3) % two_ord;
        long e = (numer / 2) % ord12;

        return mul_mod(8 % mod, pow_mod(12, e, mod), mod);
    }

    public static String solve() {
        long n0 = 10000;

        long a_mod_13_4 = c_mod_13_power_from_n_residue(n0, false, 4);
        long a_mod_12x13_4 = crt_coprime(0, 12, a_mod_13_4, pow_u64(13, 4));

        long b_mod_13_6 = c_mod_13_power_from_n_residue(a_mod_12x13_4, true, 6);
        long b_mod_12x13_6 = crt_coprime(0, 12, b_mod_13_6, pow_u64(13, 6));

        long result = c_mod_13_power_from_n_residue(b_mod_12x13_6, true, 8);
        return String.valueOf(result);
    }

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