Problem 908: Clock Sequence II

Project Euler Problem 908 Solution

EulerSolve provides an optimized solution for Project Euler Problem 908, Clock Sequence II, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $T_n=\dfrac{n(n+1)}{2}$. On a clock with $s$ positions, taking successive steps of lengths $1,2,3,\dots$ lands on the residues $T_n \pmod{s}$. Define $$R_s=\left\{T_n \bmod s : n \ge 0\right\}, \qquad U(s)=|R_s|.$$ The whole problem is driven by these forced residue sets. For each modulus $s$, the triangular walk determines the mandatory positions $R_s$; larger clock sequences are obtained by adjoining extra positions that the walk never visits. The target quantity is the primitive cumulative count $C(10000)$ modulo $1111211113$. A direct search over all moduli and all candidate sequences would be wasteful. The implementations instead derive an explicit multiplicative formula for $U(s)$, enumerate only the moduli with $U(s)\le n$, and then assemble the final count through binomial contributions and Möbius inversion. Mathematical Approach The solution has two genuinely different parts. First, determine how many triangular residues each modulus has. Second, convert that arithmetic information into a count of primitive clock sequences. The Residue Set Forced by the Triangular Walk Fix a modulus $s$. The set $R_s$ is forced: every admissible sequence attached to this modulus must contain all residues reached by the triangular walk....

Detailed mathematical approach

Problem Summary

Let $T_n=\dfrac{n(n+1)}{2}$. On a clock with $s$ positions, taking successive steps of lengths $1,2,3,\dots$ lands on the residues $T_n \pmod{s}$. Define

$$R_s=\left\{T_n \bmod s : n \ge 0\right\}, \qquad U(s)=|R_s|.$$

The whole problem is driven by these forced residue sets. For each modulus $s$, the triangular walk determines the mandatory positions $R_s$; larger clock sequences are obtained by adjoining extra positions that the walk never visits. The target quantity is the primitive cumulative count $C(10000)$ modulo $1111211113$.

A direct search over all moduli and all candidate sequences would be wasteful. The implementations instead derive an explicit multiplicative formula for $U(s)$, enumerate only the moduli with $U(s)\le n$, and then assemble the final count through binomial contributions and Möbius inversion.

Mathematical Approach

The solution has two genuinely different parts. First, determine how many triangular residues each modulus has. Second, convert that arithmetic information into a count of primitive clock sequences.

The Residue Set Forced by the Triangular Walk

Fix a modulus $s$. The set $R_s$ is forced: every admissible sequence attached to this modulus must contain all residues reached by the triangular walk. If the final sequence has size $p$, then exactly $p-U(s)$ additional residues must be chosen from the $s-U(s)$ positions outside $R_s$. Therefore a single modulus contributes

$$\binom{s-U(s)}{p-U(s)}$$

to the number of sequences of exact size $p$. Summing over all moduli with $U(s)\le p$ gives the exact-size count

$$B(p)=\sum_{s:\,U(s)\le p}\binom{s-U(s)}{p-U(s)}.$$

The implementations then prefix-sum these exact counts:

$$A(n)=\sum_{p=1}^{n} B(p).$$

What the problem asks for is the primitive version, extracted at the end by Möbius inversion:

$$C(n)=\sum_{d=1}^{n}\mu(d)\,A\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\pmod{1111211113}.$$

This is the place where repeated copies of a shorter clock sequence are removed.

Odd Prime Powers Become a Square-Residue Problem

The key identity is

$$8T_n+1=(2n+1)^2.$$

For an odd prime power $p^e$, multiplication by $8$ and addition of $1$ are bijections modulo $p^e$. So a residue $x$ is triangular modulo $p^e$ exactly when $8x+1$ is a square modulo $p^e$. Counting triangular residues is therefore the same as counting square residues.

Every nonzero square modulo $p^e$ has even $p$-adic valuation $2t$. After factoring out $p^{2t}$, the remaining unit must be a square modulo $p^{e-2t}$. For odd $p$, exactly half of the units modulo $p^m$ are squares, so the contribution of valuation $2t$ is

$$\frac{p^{e-2t}-p^{e-2t-1}}{2}.$$

Including the zero residue gives

$$U(p^e)=1+\sum_{t=0}^{\lfloor (e-1)/2\rfloor}\frac{p^{e-2t}-p^{e-2t-1}}{2}\qquad(p \text{ odd}).$$

From this one gets the prime-power rules used by the implementations:

$$U(p)=\frac{p+1}{2},$$

$$U(p^e)=p\,U(p^{e-1})-\begin{cases} p-1, & e \text{ even},\\[4pt] \dfrac{p-1}{2}, & e \text{ odd}, \end{cases} \qquad (p \text{ odd},\ e\ge 2),$$

together with the equivalent closed form

$$U(p^e)=\left\lfloor\frac{p^{e+1}}{2p+2}+1\right\rfloor \qquad (p \text{ odd}).$$

Powers of Two Cover Every Residue

The $2$-power case is different. Here the triangular walk is actually surjective: $U(2^e)=2^e$. A clean proof is to look at the first $2^e$ triangular numbers. If $0\le a<b<2^e$ and $T_a\equiv T_b \pmod{2^e}$, then

$$T_b-T_a=\frac{(b-a)(a+b+1)}{2}$$

is divisible by $2^e$, so $(b-a)(a+b+1)$ is divisible by $2^{e+1}$. But the two factors have opposite parity, and both are strictly smaller than $2^{e+1}$. The only way their product can contain $2^{e+1}$ is the trivial case $a=b$, which is impossible. Hence $T_0,T_1,\dots,T_{2^e-1}$ are distinct modulo $2^e$, so every residue class occurs exactly once.

Multiplicativity Across Coprime Moduli

If $a$ and $b$ are coprime, the Chinese remainder theorem identifies a residue modulo $ab$ with one residue modulo $a$ and one residue modulo $b$. On odd prime-power factors, triangular residues are precisely the affine images of square residues; on the $2$-part, every residue is allowed. So the condition of being triangular splits independently across coprime factors, which implies

$$U(ab)=U(a)U(b)\qquad(\gcd(a,b)=1).$$

Therefore, for a factorization

$$s=\prod_i p_i^{e_i},$$

we have the multiplicative formula

$$U(s)=\prod_i U\!\left(p_i^{e_i}\right).$$

This is the number-theoretic core of the entire algorithm: once prime powers are understood, all composite moduli follow immediately.

From $U(s)$ to the Primitive Count

For a known state $(s,U(s))$, write

$$f=s-U(s).$$

The available extra positions are exactly these $f$ untouched residues, so the modulus contributes the binomial row

$$\binom{f}{0},\binom{f}{1},\binom{f}{2},\dots,\binom{f}{\min(f,n-U(s))}$$

to the exact sizes $U(s),U(s)+1,\dots$. After all moduli have been processed, a prefix sum turns $B$ into $A$, and Möbius inversion turns $A$ into the primitive count $C$.

Worked Example: $n=4$

The moduli with $U(s)\le 4$ are

$$s\in\{1,2,3,4,5,6,7,9\}.$$

Their triangular-residue counts are

$$U(1)=1,\ U(2)=2,\ U(3)=2,\ U(4)=4,\ U(5)=3,\ U(6)=4,\ U(7)=4,\ U(9)=4.$$

Two samples explain where these numbers come from. First, modulo $5$,

$$R_5=\{0,1,3\},$$

so $U(5)=3$. Second, modulo $9$, the prime-power recurrence gives

$$U(9)=3\,U(3)-(3-1)=3\cdot 2-2=4.$$

For $s=5$, there are $5-3=2$ untouched residues, so this modulus contributes

$$\binom{2}{0}=1$$

to size $3$, and

$$\binom{2}{1}=2$$

to size $4$. Summing all moduli yields

$$B(1)=1,\qquad B(2)=2,\qquad B(3)=2,\qquad B(4)=6,$$

hence

$$A(1)=1,\qquad A(2)=3,\qquad A(3)=5,\qquad A(4)=11.$$

Now apply Möbius inversion:

$$C(4)=\mu(1)A(4)+\mu(2)A(2)+\mu(3)A(1)+\mu(4)A(1)=11-3-1+0=7.$$

This is exactly one of the sanity checks used by the solution.

How the Code Works

The C++, Python, and Java implementations all follow the same arithmetic pipeline. The differences are mostly in engineering details, not in the underlying mathematics.

Enumerating Only Relevant States

The implementations first build the prime table needed for multiplicative state generation and the Möbius table needed for the final inversion. They then enumerate states of the form $(s,U(s))$, but only when $U(s)\le n$, because larger values can never contribute to $C(n)$.

Powers of $2$ are the natural starting states since $U(2^e)=2^e$ is immediate. From a current state with value $u=U(s)$, a new odd prime $p$ is worth trying only if its smallest possible factor satisfies

$$U(p)=\frac{p+1}{2}\le \frac{n}{u},$$

that is,

$$p\le 2\left\lfloor\frac{n}{u}\right\rfloor-1.$$

Higher powers $p^e$ are extended as long as the new product $u\,U(p^e)$ stays within the bound. Because odd primes are appended in increasing order, every admissible factorization is generated exactly once.

Accumulating the Exact-Size Contributions

For each state, the number of free residues is $f=s-U(s)$. That state contributes a whole run of binomial coefficients to the exact-size table. Consecutive terms are updated by

$$\binom{f}{k+1}=\binom{f}{k}\cdot\frac{f-k}{k+1},$$

performed modulo $1111211113$. Since the modulus is prime, the needed inverses of $1,2,\dots,n+1$ exist and are precomputed once. After all states have contributed, the exact-size array is prefix-summed into the cumulative array $A$.

Final Primitive Extraction and Checks

The last stage evaluates

$$C(n)=\sum_{d=1}^{n}\mu(d)\,A\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\pmod{1111211113}.$$

The C++ and Java implementations split the expensive state-accumulation work across several workers; the Python implementation performs the same arithmetic sequentially. The solution strategy is also sanity-checked on small cases: the triangular-residue formula matches direct enumeration on small moduli, and the counts $C(3)=3$, $C(4)=7$, and $C(10)=561$ agree with the derived formulas.

Complexity Analysis

Let $S(n)$ be the number of enumerated states $(s,U(s))$ with $U(s)\le n$. Building the prime sieve up to $2n$ and the Möbius array up to $n$ costs $O(n\log\log n)$ time and $O(n)$ memory.

State generation itself costs $O(S(n))$. A state $(s,u)$ contributes

$$1+\min(s-u,n-u)$$

binomial updates, so the total running time is

$$O\!\left(n\log\log n + S(n) + \sum_{(s,u)}\bigl(1+\min(s-u,n-u)\bigr)\right).$$

The memory usage is $O(n+S(n))$. For $n=10000$, the accumulation phase dominates, which is why the threaded implementations parallelize it.

Footnotes and References

Problem page: Project Euler 908
Triangular numbers: Wikipedia - Triangular number
Quadratic residues: Wikipedia - Quadratic residue
Chinese remainder theorem: Wikipedia - Chinese remainder theorem
Möbius inversion: Wikipedia - Möbius inversion formula

Problem 908 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <limits>
#include <pthread.h>
#include <unistd.h>
#include <utility>
#include <vector>

namespace {

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

constexpr i64 kMod = 1'111'211'113LL;

i64 mod_pow(i64 a, i64 e) {
    i64 r = 1 % kMod;
    i64 x = a % kMod;
    i64 p = e;
    while (p > 0) {
        if (p & 1) {
            r = static_cast<i64>((__int128)r * x % kMod);
        }
        x = static_cast<i64>((__int128)x * x % kMod);
        p >>= 1;
    }
    return r;
}

std::vector<int> primes_up_to(int n) {
    std::vector<bool> is_prime(n + 1, true);
    if (n >= 0) {
        is_prime[0] = false;
    }
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int p = 2; (i64)p * p <= n; ++p) {
        if (!is_prime[p]) {
            continue;
        }
        for (int q = p * p; q <= n; q += p) {
            is_prime[q] = false;
        }
    }
    std::vector<int> primes;
    for (int x = 2; x <= n; ++x) {
        if (is_prime[x]) {
            primes.push_back(x);
        }
    }
    return primes;
}

u64 triangular_residue_count_bruteforce(int s) {
    std::vector<char> seen(s, 0);
    int cnt = 0;
    for (int n = 1; n <= 2 * s; ++n) {
        const int r = (int)((1LL * n * (n + 1) / 2) % s);
        if (!seen[r]) {
            seen[r] = 1;
            ++cnt;
        }
    }
    return (u64)cnt;
}

u64 u_prime_power(int p, int e) {
    if (p == 2) {
        return u64{1} << e;
    }
    u64 u = (u64)(p + 1) / 2;
    for (int k = 2; k <= e; ++k) {
        const u64 d = (k % 2 == 0) ? (u64)(p - 1) : (u64)(p - 1) / 2;
        u = (u64)p * u - d;
    }
    return u;
}

u64 u_from_factorization(int s, const std::vector<int>& primes) {
    int x = s;
    u64 u = 1;
    for (int p : primes) {
        if ((i64)p * p > x) {
            break;
        }
        if (x % p != 0) {
            continue;
        }
        int e = 0;
        while (x % p == 0) {
            x /= p;
            ++e;
        }
        u *= u_prime_power(p, e);
    }
    if (x > 1) {
        u *= u_prime_power(x, 1);
    }
    return u;
}

std::vector<int> mobius_array(int n) {
    std::vector<int> mu(n + 1, 0);
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            const i64 v = 1LL * i * p;
            if (v > n || p > spf[i]) {
                break;
            }
            spf[(int)v] = p;
            if (i % p == 0) {
                mu[(int)v] = 0;
                break;
            }
            mu[(int)v] = -mu[i];
        }
    }
    return mu;
}

struct State908 {
    u64 integer = 0;
    int count = 0;
    int prime_index = 0;
};

std::vector<State908> enumerate_states(int n_limit) {
    const std::vector<int> primes = primes_up_to(2 * n_limit);
    std::vector<State908> states;
    states.reserve((std::size_t)n_limit * 20);

    for (u64 x = 1; x <= (u64)2 * n_limit; x <<= 1) {
        states.push_back(State908{x, (int)x, 1});
    }

    for (std::size_t head = 0; head < states.size(); ++head) {
        const State908 st = states[head];
        if (st.count > n_limit) {
            continue;
        }
        const int max_new_value_factor = n_limit / st.count;
        if (max_new_value_factor < 2) {
            continue;
        }
        const int max_prime = 2 * max_new_value_factor - 1;
        for (int new_prime_index = st.prime_index; new_prime_index < (int)primes.size(); ++new_prime_index) {
            const int prime = primes[new_prime_index];
            if (prime > max_prime) {
                break;
            }
            u64 new_integer = st.integer * (u64)prime;
            u64 prime_power = (u64)prime;
            int new_value_factor = (prime + 1) / 2;
            while (new_value_factor <= max_new_value_factor) {
                const int new_value = st.count * new_value_factor;
                states.push_back(State908{new_integer, new_value, new_prime_index + 1});

                const u128 next_prime_power = (u128)prime_power * (u64)prime;
                if (next_prime_power > (u128)std::numeric_limits<u64>::max()) {
                    break;
                }
                prime_power = (u64)next_prime_power;
                new_value_factor = (int)(((u64)prime * prime_power) / (2ULL * (u64)prime + 2ULL) + 1ULL);

                const u128 next_integer = (u128)new_integer * (u64)prime;
                if (next_integer > (u128)std::numeric_limits<u64>::max()) {
                    break;
                }
                new_integer = (u64)next_integer;
            }
        }
    }

    return states;
}

i64 count_clock_sequences(int n_limit) {
    const auto states = enumerate_states(n_limit);

    std::vector<i64> inv(n_limit + 2, 0);
    inv[1] = 1;
    for (int i = 2; i <= n_limit + 1; ++i) {
        inv[i] = (kMod - (kMod / i) * inv[kMod % i] % kMod) % kMod;
    }

    std::vector<i64> rep(n_limit + 1, 0);

    struct Task {
        const std::vector<State908>* states = nullptr;
        const std::vector<i64>* inv = nullptr;
        int n_limit = 0;
        std::size_t begin = 0;
        std::size_t end = 0;
        std::vector<i64> local_rep;
    };

    auto worker = [](void* arg) -> void* {
        auto* task = static_cast<Task*>(arg);
        task->local_rep.assign((std::size_t)task->n_limit + 1, 0);
        const auto& states_ref = *task->states;
        const auto& inv_ref = *task->inv;
        const int n_limit_local = task->n_limit;

        for (std::size_t idx = task->begin; idx < task->end; ++idx) {
            const auto& st = states_ref[idx];
            if (st.count > n_limit_local) {
                continue;
            }
            const u64 free = st.integer - (u64)st.count;
            const int max_k = (int)std::min<u64>(free, (u64)(n_limit_local - st.count));

            i64 comb = 1;
            int p = st.count;
            for (int k = 0; k <= max_k; ++k) {
                i64& cell = task->local_rep[p];
                cell += comb;
                if (cell >= kMod) {
                    cell -= kMod;
                }
                if (k == max_k) {
                    break;
                }
                comb = static_cast<i64>((__int128)comb * (i64)(free - (u64)k) % kMod);
                comb = static_cast<i64>((__int128)comb * inv_ref[k + 1] % kMod);
                ++p;
            }
        }
        return nullptr;
    };

    long cpu_count = ::sysconf(_SC_NPROCESSORS_ONLN);
    std::size_t thread_count = 1;
    if (cpu_count > 1) {
        thread_count = static_cast<std::size_t>(cpu_count);
        if (thread_count > 16) {
            thread_count = 16;
        }
        if (thread_count > states.size()) {
            thread_count = states.size();
        }
        if (thread_count * 256 > states.size()) {
            thread_count = std::max<std::size_t>(1, states.size() / 256);
        }
        if (thread_count == 0) {
            thread_count = 1;
        }
    }

    if (thread_count == 1) {
        Task task;
        task.states = &states;
        task.inv = &inv;
        task.n_limit = n_limit;
        task.begin = 0;
        task.end = states.size();
        worker(&task);
        rep = std::move(task.local_rep);
    } else {
        std::vector<pthread_t> threads(thread_count);
        std::vector<Task> tasks(thread_count);
        const std::size_t base = states.size() / thread_count;
        const std::size_t rem = states.size() % thread_count;
        std::size_t cur = 0;
        for (std::size_t i = 0; i < thread_count; ++i) {
            const std::size_t len = base + (i < rem ? 1 : 0);
            tasks[i].states = &states;
            tasks[i].inv = &inv;
            tasks[i].n_limit = n_limit;
            tasks[i].begin = cur;
            tasks[i].end = cur + len;
            cur += len;
            const int rc = ::pthread_create(&threads[i], nullptr, worker, &tasks[i]);
            assert(rc == 0);
        }
        for (std::size_t i = 0; i < thread_count; ++i) {
            const int rc = ::pthread_join(threads[i], nullptr);
            assert(rc == 0);
            const auto& local = tasks[i].local_rep;
            for (int p = 1; p <= n_limit; ++p) {
                rep[p] += local[p];
                if (rep[p] >= kMod) {
                    rep[p] -= kMod;
                }
            }
        }
    }

    for (int i = 2; i <= n_limit; ++i) {
        rep[i] += rep[i - 1];
        if (rep[i] >= kMod) {
            rep[i] -= kMod;
        }
    }

    const auto mu = mobius_array(n_limit);
    i64 result = 0;
    for (int i = 1; i <= n_limit; ++i) {
        result += (i64)mu[i] * rep[n_limit / i];
        result %= kMod;
    }
    if (result < 0) {
        result += kMod;
    }
    return result;
}

void validate() {
    const auto primes = primes_up_to(500);
    for (int s = 1; s <= 200; ++s) {
        assert(u_from_factorization(s, primes) == triangular_residue_count_bruteforce(s));
    }

    assert(count_clock_sequences(3) == 3);
    assert(count_clock_sequences(4) == 7);
    assert(count_clock_sequences(10) == 561);
}

}  // namespace

int main() {
    validate();
    std::cout << count_clock_sequences(10'000) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1111211113; N = 10000

    def mp(a, e):
        r = 1; a %= MOD
        while e > 0:
            if e&1: r=r*a%MOD
            a=a*a%MOD; e >>= 1
        return r

    # Sieve primes up to 2*N
    limit = 2*N; is_prime = [True]*(limit+1); is_prime[0]=is_prime[1]=False
    for p in range(2, int(limit**0.5)+1):
        if is_prime[p]:
            for q in range(p*p, limit+1, p): is_prime[q] = False
    primes = [x for x in range(2, limit+1) if is_prime[x]]

    def u_pp(p, e):
        if p == 2: return 1 << e
        u = (p+1)//2
        for k in range(2, e+1):
            d = (p-1) if k%2==0 else (p-1)//2
            u = p*u - d
        return u

    def u_from_fact(s):
        x = s; u = 1
        for p in primes:
            if p*p > x: break
            if x%p != 0: continue
            e = 0
            while x%p == 0: x //= p; e += 1
            u *= u_pp(p, e)
        if x > 1: u *= u_pp(x, 1)
        return u

    # Mobius sieve
    spf = [0]*(N+1); mu = [0]*(N+1); ps = []; mu[1] = 1
    for i in range(2, N+1):
        if spf[i]==0: spf[i]=i; ps.append(i); mu[i]=-1
        for p in ps:
            if p*i>N: break
            spf[p*i]=p
            if i%p==0: mu[p*i]=0; break
            mu[p*i]=-mu[i]

    # Enumerate states
    states = []
    for x_power in range(64):
        x = 1 << x_power
        if x > 2*N: break
        states.append((x, x, 1))

    head = 0
    while head < len(states):
        integer, count, pi = states[head]; head += 1
        if count > N: continue
        max_f = N // count
        if max_f < 2: continue
        max_p = 2*max_f - 1
        for idx in range(pi, len(primes)):
            p = primes[idx]
            if p > max_p: break
            ni = integer * p; pp = p; nvf = (p+1)//2
            while nvf <= max_f:
                nv = count * nvf
                states.append((ni, nv, idx+1))
                next_pp = pp * p
                if next_pp > 10**18: break
                pp = next_pp
                nvf = int(p*pp / (2*p + 2) + 1)
                ni_next = ni * p
                if ni_next > 10**18: break
                ni = ni_next

    # Compute inv
    inv = [0]*(N+2); inv[1] = 1
    for i in range(2, N+2): inv[i] = (MOD - MOD//i*inv[MOD%i]%MOD) % MOD

    rep = [0]*(N+1)
    for integer, count, pi in states:
        if count > N: continue
        free = integer - count
        max_k = min(free, N-count)
        comb = 1; p = count
        for k in range(int(max_k)+1):
            rep[p] = (rep[p]+comb) % MOD
            if k == max_k: break
            comb = comb*(free-k) % MOD * inv[k+1] % MOD
            p += 1

    for i in range(2, N+1): rep[i] = (rep[i]+rep[i-1]) % MOD

    result = 0
    for i in range(1, N+1):
        result += mu[i]*rep[N//i]
        result %= MOD
    if result < 0: result += MOD
    return str(result)

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

Java

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

public class Euler908 {

    static final long kMod = 1111211113L;

    static List<Integer> primesUpTo(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int p = 2; (long) p * p <= n; ++p) {
            if (isPrime[p]) {
                for (int q = p * p; q <= n; q += p) {
                    isPrime[q] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int x = 2; x <= n; ++x) {
            if (isPrime[x]) {
                primes.add(x);
            }
        }
        return primes;
    }

    static int[] mobiusArray(int n) {
        int[] mu = new int[n + 1];
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>();
        mu[1] = 1;
        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                long v = 1L * i * p;
                if (v > n || p > spf[i]) {
                    break;
                }
                spf[(int) v] = p;
                if (i % p == 0) {
                    mu[(int) v] = 0;
                    break;
                }
                mu[(int) v] = -mu[i];
            }
        }
        return mu;
    }

    static class State908 {
        long integer;
        int count;
        int primeIndex;

        State908(long integer, int count, int primeIndex) {
            this.integer = integer;
            this.count = count;
            this.primeIndex = primeIndex;
        }
    }

    static List<State908> enumerateStates(int nLimit) {
        List<Integer> primes = primesUpTo(2 * nLimit);
        List<State908> states = new ArrayList<>();

        for (long x = 1; x <= 2L * nLimit; x <<= 1) {
            states.add(new State908(x, (int) x, 1));
        }

        for (int head = 0; head < states.size(); ++head) {
            State908 st = states.get(head);
            if (st.count > nLimit) {
                continue;
            }
            int maxNewValueFactor = nLimit / st.count;
            if (maxNewValueFactor < 2) {
                continue;
            }
            int maxPrime = 2 * maxNewValueFactor - 1;
            for (int newPrimeIndex = st.primeIndex; newPrimeIndex < primes.size(); ++newPrimeIndex) {
                int prime = primes.get(newPrimeIndex);
                if (prime > maxPrime) {
                    break;
                }
                long newInteger = st.integer * prime;
                long primePower = prime;
                int newValueFactor = (prime + 1) / 2;

                while (newValueFactor <= maxNewValueFactor) {
                    int newValue = st.count * newValueFactor;
                    states.add(new State908(newInteger, newValue, newPrimeIndex + 1));

                    java.math.BigInteger nextPrimePower = java.math.BigInteger.valueOf(primePower)
                            .multiply(java.math.BigInteger.valueOf(prime));
                    if (nextPrimePower.compareTo(java.math.BigInteger.valueOf(Long.MAX_VALUE)) > 0) {
                        break;
                    }
                    primePower = nextPrimePower.longValue();

                    java.math.BigInteger numerator = java.math.BigInteger.valueOf(prime)
                            .multiply(java.math.BigInteger.valueOf(primePower));
                    long denom = 2L * prime + 2L;
                    newValueFactor = (int) (numerator.divide(java.math.BigInteger.valueOf(denom)).longValue() + 1L);

                    java.math.BigInteger nextInteger = java.math.BigInteger.valueOf(newInteger)
                            .multiply(java.math.BigInteger.valueOf(prime));
                    if (nextInteger.compareTo(java.math.BigInteger.valueOf(Long.MAX_VALUE)) > 0) {
                        break;
                    }
                    newInteger = nextInteger.longValue();
                }
            }
        }

        return states;
    }

    static long countClockSequences(int nLimit) {
        List<State908> states = enumerateStates(nLimit);

        long[] inv = new long[nLimit + 2];
        inv[1] = 1;
        for (int i = 2; i <= nLimit + 1; ++i) {
            inv[i] = (kMod - (kMod / i) * inv[(int) (kMod % i)] % kMod) % kMod;
        }

        long[] rep = new long[nLimit + 1];

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        if (threads > 16)
            threads = 16;
        if (threads > states.size())
            threads = states.size();

        long[][] localRep = new long[threads][nLimit + 1];
        Thread[] pool = new Thread[threads];
        int base = states.size() / threads;
        int rem = states.size() % threads;

        int cur = 0;
        for (int i = 0; i < threads; ++i) {
            final int tIdx = i;
            final int start = cur;
            final int len = base + (i < rem ? 1 : 0);
            final int end = start + len;
            cur += len;

            pool[i] = new Thread(() -> {
                for (int idx = start; idx < end; ++idx) {
                    State908 st = states.get(idx);
                    if (st.count > nLimit)
                        continue;

                    long free = st.integer - st.count;
                    int maxK = (int) Math.min(free, (long) (nLimit - st.count));

                    long comb = 1;
                    int p = st.count;
                    for (int k = 0; k <= maxK; ++k) {
                        localRep[tIdx][p] += comb;
                        if (localRep[tIdx][p] >= kMod) {
                            localRep[tIdx][p] -= kMod;
                        }
                        if (k == maxK)
                            break;
                        comb = (comb * ((free - k) % kMod)) % kMod;
                        comb = (comb * inv[k + 1]) % kMod;
                        p++;
                    }
                }
            });
            pool[i].start();
        }

        try {
            for (int i = 0; i < threads; i++) {
                if (pool[i] != null)
                    pool[i].join();
            }
        } catch (InterruptedException e) {
            e.printStackTrace();
        }

        for (int i = 0; i < threads; ++i) {
            for (int p = 1; p <= nLimit; ++p) {
                rep[p] += localRep[i][p];
                if (rep[p] >= kMod) {
                    rep[p] -= kMod;
                }
            }
        }

        for (int i = 2; i <= nLimit; ++i) {
            rep[i] += rep[i - 1];
            if (rep[i] >= kMod) {
                rep[i] -= kMod;
            }
        }

        int[] mu = mobiusArray(nLimit);
        long result = 0;
        for (int i = 1; i <= nLimit; ++i) {
            result += (long) mu[i] * rep[nLimit / i];
            result %= kMod;
        }
        if (result < 0) {
            result += kMod;
        }
        return result;
    }

    public static String solve() {
        return Long.toString(countClockSequences(10000));
    }

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