Project Euler

Project Euler Problem 405: A Rectangular Tiling

We tile a rectangle whose length is twice its width using a recursive substitution rule. T(0) is the tiling consisting of a single 2x1 rectangle. For n > 0, T(n) is obtained from T(n-1) by replacin...

Source sync Apr 19, 2026

Problem #0405

Level Level 18

Solved By 740

Languages C++, Python

Answer 237696125

Length 645 words

modular_arithmeticlinear_algebrasequence

We wish to tile a rectangle whose length is twice its width.

Let $T(0)$ be the tiling consisting of a single rectangle.

For $n > 0$, let $T(n)$ be obtained from $T(n - 1)$ by replacing all tiles in the following manner:

Problem illustration

The following animation demonstrates the tilings $T(n)$ for $n$ from $0$ to $5$:

Problem animation

Let $f(n)$ be the number of points where four tiles meet in $T(n)$.

For example, $f(1) = 0$, $f(4) = 82$ and $f(10^9) \mod 17^7 = 126897180$.

Find $f(10^k)$ for $k = 10^{18}$, give your answer modulo $17^7$.

Project Euler Problem 405: A Rectangular Tiling

Solution

1. Computing Small Values by Simulation

By explicitly constructing the tilings T(0) through T(7) and counting four-corner points:

n	f(n)
0	0
1	0
2	2
3	16
4	82
5	368
6	1554
7	6384

2. Linear Recurrence

From the computed values, we find that f(n) satisfies the recurrence:

f(n) = 5*f(n-1) - 2*f(n-2) - 8*f(n-3) + 6

with initial values f(1) = 0, f(2) = 2, f(3) = 16.

The characteristic polynomial of the homogeneous part is:

x^3 - 5x^2 + 2x + 8 = (x - 4)(x - 2)(x + 1)

with roots 4, 2, and -1.

3. Closed-Form Solution

The general solution (including the constant particular solution D = 1) is:

f(n) = A * 4^n + B * 2^n + C * (-1)^n + 1

Solving with the initial values:

A = 2/5,   B = -4/3,   C = -1/15

This gives the closed form:

f(n) = (6 * 4^n - 20 * 2^n - (-1)^n + 15) / 15

4. Matrix Exponentiation

For modular computation, rewrite the recurrence in matrix form:

| f(n+1) |   | 5  -2  -8  6 |   | f(n)   |
| f(n)   | = | 1   0   0  0 | * | f(n-1) |
| f(n-1) |   | 0   1   0  0 |   | f(n-2) |
|   1    |   | 0   0   0  1 |   |   1    |

Then f(n) is the first component of M^(n-3) applied to the vector [16, 2, 0, 1]^T.

5. Computing f(10^(10^18)) mod 17^7

Step A: Order of M modulo 17^7.

The order of M modulo 17 is 8 (verified by brute force). By Hensel lifting, the order modulo 17^e is 8 * 17^(e-1) for e >= 1. Thus the order modulo 17^7 is:

ord = 8 * 17^6 = 193,100,552

Step B: Reduce the exponent.

We need (10^(10^18) - 3) mod 193100552. Factor: 193100552 = 2^3 * 17^6.

By CRT:

mod 8: 10^(10^18) = 0 mod 8 (since 10^k = 0 mod 8 for k >= 3).
mod 17^6: Using Euler’s theorem with phi(17^6) = 22717712, 10^(10^18) mod 17^6 = 10^(10^18 mod 22717712) mod 17^6 = 7544618.
CRT gives: 10^(10^18) mod 193100552 = 152370032.
Final exponent: (152370032 - 3) mod 193100552 = 152370029.

Step C: Matrix power.

Compute M^152370029 mod 17^7 and apply to [16, 2, 0, 1]^T.

Step D: Closed-form verification.

Using the closed form directly:

N = 10^(10^18) is even, so (-1)^N = 1.
Compute 2^N mod 17^7 and 4^N mod 17^7 via modular exponentiation with exponent reduced modulo phi(17^7) = 386201104.
f(N) = (6 * 4^N - 20 * 2^N - 1 + 15) * 15^(-1) mod 17^7.

Both methods yield the same result.

Answer

f(10^(10^18)) mod 17^7 = 237696125

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity

Closed-form method: O(log N * log(mod)) for modular exponentiation, which is O(log(10^18) * log(17^7)) = extremely fast.
Matrix method: O(4^3 * log(exponent)) = also very fast after exponent reduction.

Key Insights

The substitution rule alternates tile orientations, creating a fractal-like pattern.
Four-corner counts follow a linear recurrence with characteristic roots 4, 2, -1.
The factor of 4 (dominant root) reflects the 4x scaling of interior area per step.
Modular computation requires finding the multiplicative order of the matrix (or the period of 2 and 4) modulo 17^7, using Hensel lifting from mod 17.

Answer

\boxed{237696125}

C++ project_euler/problem_405/solution.cpp

/*
 * Project Euler Problem 405: A Rectangular Tiling
 *
 * A rectangle (length = 2 * width) is tiled recursively via substitution.
 * f(n) = number of interior points where exactly 4 tiles meet in T(n).
 *
 * Recurrence:  f(n) = 5*f(n-1) - 2*f(n-2) - 8*f(n-3) + 6
 * Closed form: f(n) = (6*4^n - 20*2^n - (-1)^n + 15) / 15
 *
 * Task: f(10^(10^18)) mod 17^7
 * Answer: 237696125
 *
 * Approach: closed-form modular arithmetic with Euler's theorem
 * for exponent reduction.
 *
 * Compile: g++ -O2 -std=c++17 -o solution solution.cpp
 */

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef vector<vector<ll>> Matrix;

// -----------------------------------------------------------------------
// Modular arithmetic helpers
// -----------------------------------------------------------------------

ll mod(ll a, ll m) {
    return ((a % m) + m) % m;
}

ll power(ll base, ll exp, ll m) {
    ll result = 1;
    base = mod(base, m);
    while (exp > 0) {
        if (exp & 1) result = result * base % m;
        base = base * base % m;
        exp >>= 1;
    }
    return result;
}

ll mod_inv(ll a, ll m) {
    return power(a, m - 2, m);  // Works when m is prime; for prime power, use extended gcd
}

// Extended GCD for modular inverse when m is not prime
ll ext_gcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

ll mod_inv_general(ll a, ll m) {
    ll x, y;
    ll g = ext_gcd(a % m, m, x, y);
    assert(g == 1);
    return mod(x, m);
}

// -----------------------------------------------------------------------
// 4x4 matrix operations modulo m
// -----------------------------------------------------------------------

Matrix mat_mul(const Matrix &A, const Matrix &B, ll m) {
    int n = A.size();
    Matrix C(n, vector<ll>(n, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++) if (A[i][k])
            for (int j = 0; j < n; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % m;
    return C;
}

Matrix mat_pow(Matrix M, ll p, ll m) {
    int n = M.size();
    Matrix result(n, vector<ll>(n, 0));
    for (int i = 0; i < n; i++) result[i][i] = 1;
    while (p > 0) {
        if (p & 1) result = mat_mul(result, M, m);
        M = mat_mul(M, M, m);
        p >>= 1;
    }
    return result;
}

// -----------------------------------------------------------------------
// f(n) mod m via matrix exponentiation
// -----------------------------------------------------------------------

ll f_mod_matrix(ll n, ll m) {
    if (n <= 1) return 0;
    if (n == 2) return 2 % m;
    if (n == 3) return 16 % m;
    Matrix M = {
        {mod(5, m), mod(-2, m), mod(-8, m), mod(6, m)},
        {1, 0, 0, 0},
        {0, 1, 0, 0},
        {0, 0, 0, 1}
    };
    Matrix Mp = mat_pow(M, n - 3, m);
    vector<ll> v = {16, 2, 0, 1};
    ll ans = 0;
    for (int j = 0; j < 4; j++)
        ans = (ans + Mp[0][j] * v[j]) % m;
    return mod(ans, m);
}

// -----------------------------------------------------------------------
// Solve PE 405: f(10^(10^18)) mod 17^7
// -----------------------------------------------------------------------

ll solve() {
    const ll MOD = 1;
    // 17^7
    ll mod17_7 = 1;
    for (int i = 0; i < 7; i++) mod17_7 *= 17;  // 410338673

    // phi(17^7) = 16 * 17^6
    ll mod17_6 = mod17_7 / 17;  // 24137569
    ll phi_17_7 = 16 * mod17_6; // 386201104

    // phi(17^6) = 16 * 17^5
    ll mod17_5 = mod17_6 / 17;  // 1419857
    ll phi_17_6 = 16 * mod17_5; // 22717712

    // N = 10^(10^18).  We need N mod phi(17^7).
    // phi(17^7) = 2^4 * 17^6.
    //
    // mod 16: 10^(10^18) mod 16 = 0  (since 10^k mod 16 = 0 for k >= 4)
    ll r_16 = 0;

    // mod 17^6: 10^(10^18) mod 17^6
    //         = 10^(10^18 mod phi(17^6)) mod 17^6
    //         10^18 mod phi(17^6):
    ll exp_inner = 1;
    {
        // Compute 10^18 mod phi_17_6 = 10^18 mod 22717712
        // 10^18 fits in unsigned 64-bit, just do it directly
        ll val = 1;
        for (int i = 0; i < 18; i++) {
            val = val * 10 % phi_17_6;
        }
        exp_inner = val;
    }
    ll r_17_6 = power(10, exp_inner, mod17_6);

    // CRT: x = 0 mod 16, x = r_17_6 mod 17^6
    ll inv16 = mod_inv_general(16, mod17_6);
    ll t = r_17_6 * inv16 % mod17_6;
    ll N_mod_phi = 16 * t % phi_17_7;

    // Compute f(N) using closed form:
    // f(N) = (6*4^N - 20*2^N - (-1)^N + 15) * 15^{-1}  mod 17^7
    // N is even => (-1)^N = 1

    ll pow2N = power(2, N_mod_phi, mod17_7);
    ll twoN_mod_phi = 2 * N_mod_phi % phi_17_7;
    ll pow4N = power(2, twoN_mod_phi, mod17_7);

    ll inv15 = mod_inv_general(15, mod17_7);
    ll numerator = mod(6LL * pow4N - 20LL * pow2N - 1 + 15, mod17_7);
    ll answer = numerator % mod17_7 * inv15 % mod17_7;

    return answer;
}

// -----------------------------------------------------------------------
// Verification helpers
// -----------------------------------------------------------------------

void verify_small_values() {
    cout << "--- Small values via recurrence ---\n";
    vector<ll> f = {0, 0, 2, 16};  // f[0]..f[3]
    for (int n = 4; n <= 10; n++) {
        ll fn = 5 * f[n-1] - 2 * f[n-2] - 8 * f[n-3] + 6;
        f.push_back(fn);
    }
    for (int n = 0; n <= 10; n++) {
        cout << "  f(" << n << ") = " << f[n] << "\n";
    }
}

void verify_closed_form() {
    cout << "\n--- Closed form check ---\n";
    for (int n = 1; n <= 10; n++) {
        // f(n) = (6*4^n - 20*2^n - (-1)^n + 15) / 15
        ll val = 6;
        for (int i = 0; i < n; i++) val *= 4;
        ll p2 = 1;
        for (int i = 0; i < n; i++) p2 *= 2;
        val -= 20 * p2;
        val -= (n % 2 == 0 ? 1 : -1);
        val += 15;
        val /= 15;
        cout << "  f(" << n << ") = " << val << "\n";
    }
}

void verify_f_10_9() {
    ll mod17_7 = 1;
    for (int i = 0; i < 7; i++) mod17_7 *= 17;

    ll result = f_mod_matrix(1000000000LL, mod17_7);
    cout << "\n--- Modular check ---\n";
    cout << "  f(10^9) mod 17^7 = " << result;
    cout << (result == 126897180 ? "  [OK]" : "  [FAIL]") << "\n";
}

// -----------------------------------------------------------------------
// Main
// -----------------------------------------------------------------------

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    cout << string(60, '=') << "\n";
    cout << "Project Euler Problem 405: A Rectangular Tiling\n";
    cout << string(60, '=') << "\n\n";

    verify_small_values();
    verify_closed_form();
    verify_f_10_9();

    ll answer = solve();
    cout << "\n" << string(60, '=') << "\n";
    cout << "  ANSWER: f(10^(10^18)) mod 17^7 = " << answer << "\n";
    cout << string(60, '=') << "\n";

    return 0;
}

Python project_euler/problem_405/solution.py

"""
Project Euler Problem 405: A Rectangular Tiling

We tile a rectangle (length = 2 * width) by recursively substituting every tile:
  - Horizontal 2x1 -> 4x2 block (vert | horiz | horiz | vert)
  - Vertical   1x2 -> 2x4 block (horiz | vert | vert | horiz)

f(n) = number of interior points where exactly four tiles meet in T(n).

Recurrence: f(n) = 5*f(n-1) - 2*f(n-2) - 8*f(n-3) + 6
Closed form: f(n) = (6*4^n - 20*2^n - (-1)^n + 15) / 15

Answer: f(10^(10^18)) mod 17^7 = 237696125
"""

from collections import defaultdict

# ---------------------------------------------------------------------------
# Part 1: Tile substitution simulation
# ---------------------------------------------------------------------------

def replace_tile(x1, y1, x2, y2):
    """Replace a 2:1 tile at double scale with four tiles (alternating orientation)."""
    w, h = x2 - x1, y2 - y1
    nx1, ny1 = 2 * x1, 2 * y1
    nw, nh = 2 * w, 2 * h
    if w > h:  # horizontal tile
        return [
            (nx1, ny1, nx1 + h, ny1 + nh),            # left vertical
            (nx1 + h, ny1, nx1 + nw - h, ny1 + h),    # top-centre horizontal
            (nx1 + h, ny1 + h, nx1 + nw - h, ny1 + nh),  # bottom-centre horizontal
            (nx1 + nw - h, ny1, nx1 + nw, ny1 + nh),  # right vertical
        ]
    else:  # vertical tile
        return [
            (nx1, ny1, nx1 + nw, ny1 + w),            # top horizontal
            (nx1, ny1 + w, nx1 + w, ny1 + nh - w),    # left vertical
            (nx1 + w, ny1 + w, nx1 + nw, ny1 + nh - w),  # right vertical
            (nx1, ny1 + nh - w, nx1 + nw, ny1 + nh),  # bottom horizontal
        ]

def iterate(tiles):
    """Apply one substitution step to every tile."""
    result = []
    for t in tiles:
        result.extend(replace_tile(*t))
    return result

def count_four_corners(tiles):
    """Count interior points where exactly 4 tile corners meet."""
    corner_count = defaultdict(int)
    for x1, y1, x2, y2 in tiles:
        for pt in [(x1, y1), (x2, y1), (x1, y2), (x2, y2)]:
            corner_count[pt] += 1

    min_x = min(t[0] for t in tiles)
    max_x = max(t[2] for t in tiles)
    min_y = min(t[1] for t in tiles)
    max_y = max(t[3] for t in tiles)

    return sum(
        1 for (x, y), c in corner_count.items()
        if min_x < x < max_x and min_y < y < max_y and c == 4
    )

# ---------------------------------------------------------------------------
# Part 2: Recurrence / closed-form helpers
# ---------------------------------------------------------------------------

def f_recurrence(n):
    """Compute f(n) using the linear recurrence (exact integers)."""
    if n <= 0:
        return 0
    a, b, c = 0, 0, 0  # f(0), f(-1), f(-2) -- but we only need from f(1)
    vals = {0: 0, 1: 0, 2: 2, 3: 16}
    if n in vals:
        return vals[n]
    f1, f2, f3 = 16, 2, 0  # f(3), f(2), f(1)
    for _ in range(4, n + 1):
        fn = 5 * f1 - 2 * f2 - 8 * f3 + 6
        f3, f2, f1 = f2, f1, fn
    return f1

def f_closed(n):
    """Compute f(n) via the closed-form expression (exact integers)."""
    return (6 * 4**n - 20 * 2**n - (-1)**n + 15) // 15

# ---------------------------------------------------------------------------
# Part 3: Matrix exponentiation modulo m
# ---------------------------------------------------------------------------

def mat_mul(A, B, mod):
    """Multiply two 4x4 matrices modulo mod."""
    n = len(A)
    return [
        [sum(A[i][k] * B[k][j] for k in range(n)) % mod for j in range(n)]
        for i in range(n)
    ]

def mat_pow(M, p, mod):
    """Raise a 4x4 matrix to power p modulo mod via binary exponentiation."""
    n = len(M)
    result = [[int(i == j) for j in range(n)] for i in range(n)]
    base = [row[:] for row in M]
    while p > 0:
        if p & 1:
            result = mat_mul(result, base, mod)
        base = mat_mul(base, base, mod)
        p >>= 1
    return result

def f_mod_matrix(n, mod):
    """Compute f(n) mod m using 4x4 matrix exponentiation."""
    if n <= 1:
        return 0
    if n == 2:
        return 2 % mod
    if n == 3:
        return 16 % mod
    M = [[5, -2, -8, 6],
         [1,  0,  0, 0],
         [0,  1,  0, 0],
         [0,  0,  0, 1]]
    Mp = mat_pow(M, n - 3, mod)
    v = [16, 2, 0, 1]
    return sum(Mp[0][j] * v[j] for j in range(4)) % mod

def f_mod_closed(n_exp_reducer, mod):
    """
    Compute f(N) mod m where N = 10^(10^18) using the closed form.

    f(N) = (6*4^N - 20*2^N - (-1)^N + 15) / 15  mod m

    Parameters
    ----------
    n_exp_reducer : callable
        Function that, given (base, mod), returns base^N mod m
        using the tower of exponents with Euler's theorem.
    mod : int
        The modulus.
    """
    inv15 = pow(15, -1, mod)
    pow4N = n_exp_reducer(4, mod)
    pow2N = n_exp_reducer(2, mod)
    # N = 10^(10^18) is even, so (-1)^N = 1
    neg1N = 1
    return (6 * pow4N - 20 * pow2N - neg1N + 15) * inv15 % mod

# ---------------------------------------------------------------------------
# Part 4: Solve PE 405
# ---------------------------------------------------------------------------

def solve_pe405():
    """
    Compute f(10^(10^18)) mod 17^7.

    Uses the closed form with modular exponentiation, reducing exponents
    via Euler's theorem and CRT.
    """
    MOD = 17**7  # 410338673
    phi_mod = 16 * 17**6  # phi(17^7) = 386201104

    # We need base^N mod MOD where N = 10^(10^18).
    # By Euler's theorem: base^N mod MOD = base^(N mod phi(MOD)) mod MOD
    # provided gcd(base, MOD) = 1 (true for base=2,4 since MOD is a power of 17).
    #
    # N mod phi_mod, where phi_mod = 2^4 * 17^6:
    #   mod 16:   10^(10^18) = 0 mod 16  (since 10^k = 0 mod 16 for k >= 4)
    #   mod 17^6: 10^(10^18) mod 17^6
    #             = 10^(10^18 mod phi(17^6)) mod 17^6
    #             phi(17^6) = 16 * 17^5 = 22717712
    #             10^18 mod 22717712 is computed directly.
    #             Then pow(10, that, 17^6).

    mod_17_6 = 17**6
    phi_17_6 = 16 * 17**5

    exp_inner = pow(10, 18) % phi_17_6        # 10^18 mod phi(17^6)
    r_17_6 = pow(10, exp_inner, mod_17_6)     # 10^(10^18) mod 17^6

    # CRT: x = 0 mod 16,  x = r_17_6 mod 17^6
    inv16 = pow(16, -1, mod_17_6)
    t = (r_17_6 * inv16) % mod_17_6
    N_mod_phi = (16 * t) % phi_mod            # N mod phi(17^7)

    # Now compute f(N) via closed form
    pow2N = pow(2, N_mod_phi, MOD)
    twoN_mod_phi = (2 * N_mod_phi) % phi_mod
    pow4N = pow(2, twoN_mod_phi, MOD)         # 4^N = 2^(2N)
    neg1N = 1                                  # N is even

    inv15 = pow(15, -1, MOD)
    numerator = (6 * pow4N - 20 * pow2N - neg1N + 15) % MOD
    return (numerator * inv15) % MOD

# ---------------------------------------------------------------------------
# Part 5: Visualization
# ---------------------------------------------------------------------------

def draw_tiling(tiles, ax, title="", cmap_name="Set3"):
    """Draw a set of rectangular tiles on the given Axes."""
    cmap = plt.colormaps[cmap_name]
    min_x = min(t[0] for t in tiles)
    max_x = max(t[2] for t in tiles)
    min_y = min(t[1] for t in tiles)
    max_y = max(t[3] for t in tiles)

    for i, (x1, y1, x2, y2) in enumerate(tiles):
        color = cmap(i % 12 / 12.0)
        rect = patches.Rectangle(
            (x1, y1), x2 - x1, y2 - y1,
            linewidth=1.2, edgecolor="black", facecolor=color
        )
        ax.add_patch(rect)

    ax.set_xlim(min_x - 0.2, max_x + 0.2)
    ax.set_ylim(min_y - 0.2, max_y + 0.2)
    ax.set_aspect("equal")
    ax.invert_yaxis()
    ax.set_title(title, fontsize=11)
    ax.axis("off")

def create_visualization(save_path="visualization.png"):
    """Generate a multi-panel figure showing T(0)..T(4) and f(n) growth."""

Project Euler Problem 405: A Rectangular Tiling

Problem Statement

Project Euler Problem 405: A Rectangular Tiling

Solution

1. Computing Small Values by Simulation

2. Linear Recurrence

3. Closed-Form Solution

4. Matrix Exponentiation

5. Computing f(10^(10^18)) mod 17^7

Answer

Correctness

Complexity

Key Insights

Answer

Code