Problem 440: GCD and Tiling

Project Euler Problem 440 Solution

EulerSolve provides an optimized solution for Project Euler Problem 440, GCD and Tiling, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $T(n)$ denote the number of tilings of a strip of length $n$ using $1\times 2$ dominoes and $1\times 1$ tiles colored with one of ten digits. The empty strip contributes one tiling, so $$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$ The task is to evaluate $$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$ for $L=2000$. A direct computation is impossible because there are $L^3$ triples and the indices $c^a$ become enormous almost immediately. Mathematical Approach Step 1: Rewrite the tiling sequence as a Lucas sequence Define a Lucas sequence of the first kind by $$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$ Then $$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$ so the tiling numbers are simply shifted by one position: $$T(n)=U_{n+1}.$$ This is the key structural observation, because Lucas sequences with parameter $Q=-1$ satisfy the strong divisibility law $$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$ Therefore, for arbitrary nonnegative integers $x$ and $y$, $$\gcd(T(x),T(y))=\gcd(U_{x+1},U_{y+1})=U_{\gcd(x+1,y+1)}=T\bigl(\gcd(x+1,y+1)-1\bigr).$$ So the original problem is reduced to the arithmetic of $\gcd(c^a+1,c^b+1)$. Step 2: Evaluate $\gcd(c^a+1,c^b+1)$ with 2-adic valuations Write $$a=2^r u,\qquad b=2^s v,$$ where $u$ and $v$ are odd....

Detailed mathematical approach

Problem Summary

Let $T(n)$ denote the number of tilings of a strip of length $n$ using $1\times 2$ dominoes and $1\times 1$ tiles colored with one of ten digits. The empty strip contributes one tiling, so

$$T(0)=1,\qquad T(1)=10,\qquad T(n)=10T(n-1)+T(n-2)\quad(n\ge 2).$$

The task is to evaluate

$$S(L)=\sum_{1\le a,b,c\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)\pmod{987898789}$$

for $L=2000$. A direct computation is impossible because there are $L^3$ triples and the indices $c^a$ become enormous almost immediately.

Mathematical Approach

Step 1: Rewrite the tiling sequence as a Lucas sequence

Define a Lucas sequence of the first kind by

$$U_0=0,\qquad U_1=1,\qquad U_{n+2}=10U_{n+1}+U_n.$$

Then

$$U_1=1,\quad U_2=10,\quad U_3=101,\quad U_4=1020,\dots$$

so the tiling numbers are simply shifted by one position:

$$T(n)=U_{n+1}.$$

This is the key structural observation, because Lucas sequences with parameter $Q=-1$ satisfy the strong divisibility law

$$\gcd(U_m,U_n)=U_{\gcd(m,n)}.$$

Therefore, for arbitrary nonnegative integers $x$ and $y$,

$$\gcd(T(x),T(y))=\gcd(U_{x+1},U_{y+1})=U_{\gcd(x+1,y+1)}=T\bigl(\gcd(x+1,y+1)-1\bigr).$$

So the original problem is reduced to the arithmetic of $\gcd(c^a+1,c^b+1)$.

Step 2: Evaluate $\gcd(c^a+1,c^b+1)$ with 2-adic valuations

Write

$$a=2^r u,\qquad b=2^s v,$$

where $u$ and $v$ are odd. The answer depends only on whether $r$ and $s$ are equal, i.e. on whether $\nu_2(a)=\nu_2(b)$.

If $r=s$, let $d=\gcd(a,b)$. Then we may write

$$a=d\alpha,\qquad b=d\beta,\qquad \gcd(\alpha,\beta)=1,$$

and both $\alpha$ and $\beta$ are odd. Setting $y=c^d$, we get

$$\gcd(c^a+1,c^b+1)=\gcd(y^\alpha+1,y^\beta+1)=y+1=c^d+1.$$

If $r\ne s$, factor out $d=\gcd(a,b)$ again, but now one of the quotients $a/d$, $b/d$ is odd and the other is even. Put $y=c^d$, and let $g$ be a common divisor of $y^{a/d}+1$ and $y^{b/d}+1$. Then

$$y^{a/d}\equiv -1 \pmod g,\qquad y^{b/d}\equiv -1 \pmod g.$$

Raise the first congruence to the even exponent $b/d$, and the second to the odd exponent $a/d$. This gives

$$y^{ab/d^2}\equiv 1 \pmod g,\qquad y^{ab/d^2}\equiv -1 \pmod g,$$

hence $2\equiv 0\pmod g$. So every common divisor must divide $2$. Therefore

$$\gcd(c^a+1,c^b+1)=\begin{cases} c^{\gcd(a,b)}+1,& \nu_2(a)=\nu_2(b),\\ 2,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is odd},\\ 1,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is even}. \end{cases}$$

Step 3: Translate that gcd back to tiling numbers

Substitute the previous identity into the Lucas-sequence formula. If $d=\gcd(a,b)$, then

$$\gcd\bigl(T(c^a),T(c^b)\bigr)=\begin{cases} T(c^d),& \nu_2(a)=\nu_2(b),\\ T(1)=10,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is odd},\\ T(0)=1,& \nu_2(a)\ne \nu_2(b)\text{ and }c\text{ is even}. \end{cases}$$

This is exactly the dichotomy exploited by the implementations. The huge indices disappear from the gcd itself; only the equal-valuation pairs still need values of the form $T(c^d)$.

Step 4: Rearrange the triple sum by counting exponent pairs once

For each $d\in\{1,\dots,L\}$, define the number of ordered pairs

$$C_d=\#\{(a,b):1\le a,b\le L,\ \nu_2(a)=\nu_2(b),\ \gcd(a,b)=d\}.$$

Also define the complementary count

$$C_{\mathrm{diff}}=L^2-\sum_{d=1}^{L} C_d,$$

which counts the ordered pairs with different 2-adic valuations. For a fixed base $c$, let

$$B(c)=\begin{cases} 10,& c\text{ odd},\\ 1,& c\text{ even}. \end{cases}$$

Then the entire contribution of this single $c$ is

$$\sum_{1\le a,b\le L}\gcd\bigl(T(c^a),T(c^b)\bigr)=C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d).$$

Summing over $c$ yields the final rearrangement

$$\boxed{S(L)=\sum_{c=1}^{L}\left(C_{\mathrm{diff}}\,B(c)+\sum_{d=1}^{L} C_d\,T(c^d)\right).}$$

Now the pair statistics $\{C_d\}$ are independent of $c$, so they can be precomputed once.

Step 5: Compute $T(c^d)\bmod M$ by period reduction

Let $M=987898789$. The recurrence can be encoded with the companion matrix

$$A=\begin{bmatrix}10&1\\ 1&0\end{bmatrix},\qquad \begin{bmatrix}T(n+1)\\ T(n)\end{bmatrix}=A^n\begin{bmatrix}10\\ 1\end{bmatrix}.$$

Modulo $M$, the state space is finite, so the sequence is periodic. The implementations compute the exact period from the matrix order. Since $M$ is prime and the discriminant

$$\Delta=10^2+4=104$$

is a quadratic residue modulo $M$, namely

$$\left(\frac{104}{M}\right)=1,$$

the eigenvalues of $A$ lie in $\mathbb{F}_M$, so the order of $A$ divides $M-1$. The algorithm starts from the candidate $M-1=987898788$, factors it, and removes prime factors whenever the smaller exponent already gives $A^k=I$. For this modulus no reduction occurs, so the period is

$$\pi=987898788.$$

Hence every huge index is reduced as

$$T(n)\equiv T(n\bmod \pi)\pmod M.$$

For each fixed $c$, the reduced powers are generated iteratively:

$$p_1\equiv c,\qquad p_{d+1}\equiv p_d\,c \pmod \pi,$$

so the program never forms $c^d$ as a gigantic integer.

Worked example: $L=2$

When $L=2$, the ordered pairs $(a,b)$ are $(1,1)$, $(1,2)$, $(2,1)$, $(2,2)$. Equal 2-adic valuation occurs only for $(1,1)$ and $(2,2)$, so

$$C_1=1,\qquad C_2=1,\qquad C_{\mathrm{diff}}=2.$$

For $c=1$, every term $T(1^d)$ equals $T(1)=10$, and $B(1)=10$. Thus

$$2\cdot 10 + 1\cdot 10 + 1\cdot 10 = 40.$$

For $c=2$, we have $B(2)=1$, $T(2)=101$, and $T(4)=10301$, so

$$2\cdot 1 + 1\cdot 101 + 1\cdot 10301 = 10404.$$

Therefore

$$S(2)=40+10404=10444,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical plan. They first precompute the 2-adic valuation of every integer from $1$ to $L$, then scan all ordered pairs $(a,b)$. If the valuations are equal, the pair is added to the bucket indexed by $\gcd(a,b)$; otherwise it contributes only to the residual count $C_{\mathrm{diff}}$.

Next they determine the recurrence period modulo $987898789$, build a table of binary powers of the companion matrix, and loop over $c=1,\dots,L$. For each fixed $c$, the reduced indices $c^1,c^2,\dots,c^L\pmod \pi$ are generated by repeated multiplication. Each reduced index is converted to $T(\cdot)\bmod M$ by matrix exponentiation, multiplied by its precomputed pair count, and added to the running sum. The C++ implementation additionally parallelizes the outer loop over $c$, but the arithmetic is identical in all three languages.

Complexity Analysis

The pair-count precomputation examines all ordered pairs $(a,b)$, so it costs $O(L^2)$ time. The main accumulation also has $L^2$ terms. Evaluating one reduced index uses a fixed table of matrix squares, so formally the cost is $O(\log \pi)$ per term; here $\pi=987898788$ is fixed by the modulus, so in practice this is a small constant. Thus the total running time is $O(L^2\log \pi)$, which is effectively $O(L^2)$ for this problem, and the memory usage is $O(L)$.

References

Problem page: https://projecteuler.net/problem=440
Lucas sequences: Wikipedia — Lucas sequence
Linear recurrences and companion matrices: Wikipedia — Linear recurrence with constant coefficients
Exponentiation by squaring: Wikipedia — Exponentiation by squaring
Legendre symbol: Wikipedia — Legendre symbol

Problem 440 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <utility>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr int MOD = 987898789;

struct Options {
    int limit = 2000;
    int threads = 0;
    bool run_checkpoints = true;
};

struct Mat2 {
    int a00 = 1;
    int a01 = 0;
    int a10 = 0;
    int a11 = 1;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit) ||
            parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1 && options.threads >= 0;
}

int choose_thread_count(int requested, std::size_t work_items) {
    if (work_items <= 1U) {
        return 1;
    }
    int threads = requested;
    if (threads <= 0) {
        threads = static_cast<int>(std::thread::hardware_concurrency());
    }
    if (threads <= 0) {
        threads = 4;
    }
    if (threads > static_cast<int>(work_items)) {
        threads = static_cast<int>(work_items);
    }
    if (threads < 1) {
        threads = 1;
    }
    return threads;
}

int mod_pow(int base, u64 exp, int mod) {
    u64 result = 1ULL;
    u64 cur = static_cast<u64>((base % mod + mod) % mod);
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * cur) % static_cast<u64>(mod));
        }
        cur = static_cast<u64>((static_cast<u128>(cur) * cur) % static_cast<u64>(mod));
        e >>= 1ULL;
    }
    return static_cast<int>(result);
}

Mat2 mul(const Mat2& x, const Mat2& y) {
    Mat2 r;
    r.a00 = static_cast<int>((static_cast<u128>(x.a00) * y.a00 + static_cast<u128>(x.a01) * y.a10) % MOD);
    r.a01 = static_cast<int>((static_cast<u128>(x.a00) * y.a01 + static_cast<u128>(x.a01) * y.a11) % MOD);
    r.a10 = static_cast<int>((static_cast<u128>(x.a10) * y.a00 + static_cast<u128>(x.a11) * y.a10) % MOD);
    r.a11 = static_cast<int>((static_cast<u128>(x.a10) * y.a01 + static_cast<u128>(x.a11) * y.a11) % MOD);
    return r;
}

Mat2 mat_pow(Mat2 base, u64 exp) {
    Mat2 result;
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = mul(result, base);
        }
        e >>= 1ULL;
        if (e > 0ULL) {
            base = mul(base, base);
        }
    }
    return result;
}

std::vector<std::pair<u64, int>> factorize(u64 n) {
    std::vector<std::pair<u64, int>> fac;
    u64 x = n;
    for (u64 p = 2; p * p <= x; ++p) {
        if (x % p != 0ULL) {
            continue;
        }
        int e = 0;
        while (x % p == 0ULL) {
            x /= p;
            ++e;
        }
        fac.push_back({p, e});
    }
    if (x > 1ULL) {
        fac.push_back({x, 1});
    }
    return fac;
}

u64 sequence_period() {
    // T(n) recurrence: T_n = 10*T_{n-1} + T_{n-2} (mod MOD).
    // Companion matrix A = [[10,1],[1,0]].
    const Mat2 A{10, 1, 1, 0};
    const int leg = mod_pow(104 % MOD, static_cast<u64>((MOD - 1) / 2), MOD);
    u64 candidate = (leg == 1) ? static_cast<u64>(MOD - 1) : static_cast<u64>(MOD + 1);

    u64 period = candidate;
    const auto fac = factorize(candidate);
    for (const auto& [q, _e] : fac) {
        while (period % q == 0ULL) {
            const u64 trial = period / q;
            const Mat2 p = mat_pow(A, trial);
            if (p.a00 == 1 && p.a01 == 0 && p.a10 == 0 && p.a11 == 1) {
                period = trial;
            } else {
                break;
            }
        }
    }
    return period;
}

std::vector<Mat2> build_powers_for_T() {
    std::vector<Mat2> pw;
    pw.reserve(64);
    pw.push_back(Mat2{10, 1, 1, 0});
    for (int i = 1; i < 64; ++i) {
        pw.push_back(mul(pw.back(), pw.back()));
    }
    return pw;
}

int T_mod_from_index(u64 n, const std::vector<Mat2>& pw) {
    if (n == 0ULL) {
        return 1;  // T(0)
    }

    int v0 = 10;  // T(1)
    int v1 = 1;   // T(0)
    u64 e = n - 1ULL;
    int bit = 0;
    while (e > 0ULL) {
        if (e & 1ULL) {
            const Mat2& m = pw[static_cast<std::size_t>(bit)];
            const int nv0 =
                static_cast<int>((static_cast<u128>(m.a00) * v0 + static_cast<u128>(m.a01) * v1) % MOD);
            const int nv1 =
                static_cast<int>((static_cast<u128>(m.a10) * v0 + static_cast<u128>(m.a11) * v1) % MOD);
            v0 = nv0;
            v1 = nv1;
        }
        e >>= 1ULL;
        ++bit;
    }
    return v0;
}

std::vector<i64> build_counts_same(int l) {
    std::vector<i64> cnt(static_cast<std::size_t>(l + 1), 0);
    std::vector<int> v2(static_cast<std::size_t>(l + 1), 0);
    for (int x = 1; x <= l; ++x) {
        v2[static_cast<std::size_t>(x)] = __builtin_ctz(static_cast<unsigned>(x));
    }

    for (int a = 1; a <= l; ++a) {
        const int ta = v2[static_cast<std::size_t>(a)];
        for (int b = 1; b <= l; ++b) {
            if (v2[static_cast<std::size_t>(b)] != ta) {
                continue;
            }
            const int g = std::gcd(a, b);
            ++cnt[static_cast<std::size_t>(g)];
        }
    }
    return cnt;
}

int solve_mod(int l, int requested_threads) {
    const u64 period = sequence_period();
    const std::vector<Mat2> pw = build_powers_for_T();
    const std::vector<i64> counts_same = build_counts_same(l);

    i64 same_total = 0;
    for (int d = 1; d <= l; ++d) {
        same_total += counts_same[static_cast<std::size_t>(d)];
    }
    const i64 total_pairs = static_cast<i64>(l) * static_cast<i64>(l);
    const i64 count_diff = total_pairs - same_total;

    std::vector<int> counts_same_mod(static_cast<std::size_t>(l + 1), 0);
    for (int d = 1; d <= l; ++d) {
        counts_same_mod[static_cast<std::size_t>(d)] =
            static_cast<int>(counts_same[static_cast<std::size_t>(d)] % MOD);
    }

    const int thread_count = choose_thread_count(requested_threads, static_cast<std::size_t>(l));
    std::vector<int> partial(static_cast<std::size_t>(thread_count), 0);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(thread_count));

    for (int t = 0; t < thread_count; ++t) {
        const int begin = l * t / thread_count + 1;
        const int end = l * (t + 1) / thread_count;

        workers.emplace_back([&, begin, end, t]() {
            int local = 0;

            for (int c = begin; c <= end; ++c) {
                const int base = (c & 1) ? 10 : 1;
                int sum_c = static_cast<int>((static_cast<u128>(count_diff % MOD) * base) % MOD);

                u64 p = 1ULL;
                for (int d = 1; d <= l; ++d) {
                    p = static_cast<u64>((static_cast<u128>(p) * static_cast<u64>(c)) % period);
                    const int tval = T_mod_from_index(p, pw);
                    sum_c = static_cast<int>((sum_c +
                                              static_cast<u128>(counts_same_mod[static_cast<std::size_t>(d)]) * tval) %
                                             MOD);
                }

                local += sum_c;
                if (local >= MOD) {
                    local -= MOD;
                }
            }

            partial[static_cast<std::size_t>(t)] = local;
        });
    }

    for (auto& worker : workers) {
        worker.join();
    }

    int ans = 0;
    for (int x : partial) {
        ans += x;
        if (ans >= MOD) {
            ans -= MOD;
        }
    }
    return ans;
}

bool run_checkpoints(int requested_threads) {
    if (solve_mod(2, requested_threads) != 10444) {
        std::cerr << "Checkpoint failed: S(2)\n";
        return false;
    }
    if (solve_mod(3, requested_threads) != 862246950) {
        std::cerr << "Checkpoint failed: S(3) mod 987898789\n";
        return false;
    }
    if (solve_mod(4, requested_threads) != 670616280) {
        std::cerr << "Checkpoint failed: S(4) mod 987898789\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints(options.threads)) {
        return 2;
    }

    std::cout << solve_mod(options.limit, options.threads) << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 987898789
    l = 2000

    def mod_pow_int(base, exp, mod=MOD):
        r = 1; base %= mod
        while exp > 0:
            if exp & 1: r = r * base % mod
            base = base * base % mod; exp >>= 1
        return r

    def mat_mul(x, y):
        return [
            (x[0]*y[0] + x[1]*y[2]) % MOD, (x[0]*y[1] + x[1]*y[3]) % MOD,
            (x[2]*y[0] + x[3]*y[2]) % MOD, (x[2]*y[1] + x[3]*y[3]) % MOD
        ]

    def mat_pow(base, exp):
        r = [1,0,0,1]
        while exp > 0:
            if exp & 1: r = mat_mul(r, base)
            exp >>= 1
            if exp > 0: base = mat_mul(base, base)
        return r

    # Sequence period
    A = [10,1,1,0]
    leg = mod_pow_int(104 % MOD, (MOD-1)//2)
    cand = MOD - 1 if leg == 1 else MOD + 1

    def factorize(n):
        fac = []; x = n
        p = 2
        while p*p <= x:
            if x % p == 0:
                e = 0
                while x % p == 0: x //= p; e += 1
                fac.append((p, e))
            p += 1
        if x > 1: fac.append((x, 1))
        return fac

    period = cand
    for q, _ in factorize(cand):
        while period % q == 0:
            trial = period // q
            p = mat_pow(A, trial)
            if p == [1,0,0,1]: period = trial
            else: break

    # Powers of matrix for fast T computation
    pw = [A[:]]
    for _ in range(63): pw.append(mat_mul(pw[-1], pw[-1]))

    def T_mod(n):
        if n == 0: return 1
        v0, v1 = 10, 1; e = n - 1; bit = 0
        while e > 0:
            if e & 1:
                m = pw[bit]
                v0, v1 = (m[0]*v0 + m[1]*v1) % MOD, (m[2]*v0 + m[3]*v1) % MOD
            e >>= 1; bit += 1
        return v0

    # Counts of pairs (a,b) with same v2 and gcd=d
    v2 = [0]*(l+1)
    for x in range(1, l+1):
        v2[x] = (x & -x).bit_length() - 1
    cnt = [0]*(l+1)
    for a in range(1, l+1):
        ta = v2[a]
        for b in range(1, l+1):
            if v2[b] != ta: continue
            cnt[math.gcd(a, b)] += 1

    same_total = sum(cnt[1:l+1])
    total_pairs = l * l
    count_diff = total_pairs - same_total

    ans = 0
    for c in range(1, l+1):
        base = 10 if c & 1 else 1
        sum_c = count_diff % MOD * base % MOD
        p = 1
        for d in range(1, l+1):
            p = p * c % period
            tv = T_mod(p)
            sum_c = (sum_c + cnt[d] % MOD * tv) % MOD
        ans = (ans + sum_c) % MOD

    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler440 {
    private static final int MOD = 987898789;
    private static final int LIMIT = 2000;

    private static class Mat2 {
        long a00, a01, a10, a11;

        Mat2(long a00, long a01, long a10, long a11) {
            this.a00 = a00;
            this.a01 = a01;
            this.a10 = a10;
            this.a11 = a11;
        }
    }

    private static long modPow(long base, long exp, long mod) {
        long result = 1;
        long cur = (base % mod + mod) % mod;
        long e = exp;
        while (e > 0) {
            if ((e & 1) != 0) {
                result = (result * cur) % mod;
            }
            cur = (cur * cur) % mod;
            e >>= 1;
        }
        return result;
    }

    private static Mat2 mul(Mat2 x, Mat2 y) {
        return new Mat2(
                (x.a00 * y.a00 + x.a01 * y.a10) % MOD,
                (x.a00 * y.a01 + x.a01 * y.a11) % MOD,
                (x.a10 * y.a00 + x.a11 * y.a10) % MOD,
                (x.a10 * y.a01 + x.a11 * y.a11) % MOD);
    }

    private static Mat2 matPow(Mat2 base, long exp) {
        Mat2 result = new Mat2(1, 0, 0, 1);
        long e = exp;
        while (e > 0) {
            if ((e & 1) != 0) {
                result = mul(result, base);
            }
            e >>= 1;
            if (e > 0) {
                base = mul(base, base);
            }
        }
        return result;
    }

    private static long sequencePeriod() {
        Mat2 A = new Mat2(10, 1, 1, 0);
        long leg = modPow(104 % MOD, (MOD - 1) / 2, MOD);
        long candidate = (leg == 1) ? (MOD - 1) : (MOD + 1);

        long period = candidate;
        long x = candidate;
        List<long[]> fac = new ArrayList<>();
        for (long p = 2; p * p <= x; p++) {
            if (x % p == 0) {
                int e = 0;
                while (x % p == 0) {
                    x /= p;
                    e++;
                }
                fac.add(new long[] { p, e });
            }
        }
        if (x > 1) {
            fac.add(new long[] { x, 1 });
        }

        for (long[] f : fac) {
            long q = f[0];
            while (period % q == 0) {
                long trial = period / q;
                Mat2 pMat = matPow(A, trial);
                if (pMat.a00 == 1 && pMat.a01 == 0 && pMat.a10 == 0 && pMat.a11 == 1) {
                    period = trial;
                } else {
                    break;
                }
            }
        }
        return period;
    }

    private static List<Mat2> buildPowersForT() {
        List<Mat2> pw = new ArrayList<>();
        pw.add(new Mat2(10, 1, 1, 0));
        for (int i = 1; i < 64; i++) {
            pw.add(mul(pw.get(i - 1), pw.get(i - 1)));
        }
        return pw;
    }

    private static long tModFromIndex(long n, List<Mat2> pw) {
        if (n == 0)
            return 1;

        long v0 = 10;
        long v1 = 1;
        long e = n - 1;
        int bit = 0;
        while (e > 0) {
            if ((e & 1) != 0) {
                Mat2 m = pw.get(bit);
                long nv0 = (m.a00 * v0 + m.a01 * v1) % MOD;
                long nv1 = (m.a10 * v0 + m.a11 * v1) % MOD;
                v0 = nv0;
                v1 = nv1;
            }
            e >>= 1;
            bit++;
        }
        return v0;
    }

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

    public static String solve() {
        int l = LIMIT;
        long period = sequencePeriod();
        List<Mat2> pw = buildPowersForT();

        long[] countsSame = new long[l + 1];
        int[] v2 = new int[l + 1];
        for (int x = 1; x <= l; x++) {
            v2[x] = Integer.numberOfTrailingZeros(x);
        }

        for (int a = 1; a <= l; a++) {
            int ta = v2[a];
            for (int b = 1; b <= l; b++) {
                if (v2[b] != ta)
                    continue;
                int g = gcd(a, b);
                countsSame[g]++;
            }
        }

        long sameTotal = 0;
        for (int d = 1; d <= l; d++) {
            sameTotal += countsSame[d];
        }
        long totalPairs = (long) l * l;
        long countDiff = totalPairs - sameTotal;

        long[] countsSameMod = new long[l + 1];
        for (int d = 1; d <= l; d++) {
            countsSameMod[d] = countsSame[d] % MOD;
        }

        long ans = 0;
        for (int c = 1; c <= l; c++) {
            long base = (c % 2 != 0) ? 10 : 1;
            long sumC = ((countDiff % MOD) * base) % MOD;

            long p = 1;
            for (int d = 1; d <= l; d++) {
                p = (p * c) % period;
                long tval = tModFromIndex(p, pw);
                sumC = (sumC + countsSameMod[d] * tval) % MOD;
            }
            ans = (ans + sumC) % MOD;
        }

        return String.valueOf(ans);
    }

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