Problem 929: Odd-Run Compositions

Project Euler Problem 929 Solution

EulerSolve provides an optimized solution for Project Euler Problem 929, Odd-Run Compositions, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary An odd-run composition of $n$ is a composition $n=p_1+p_2+\cdots+p_m$ in which every maximal block of equal adjacent parts has odd length. For example, $2+1+1+1$ is valid because the run of 1's has length 3, while $2+2+1$ is invalid because the run of 2's has length 2. The problem asks for the number $u_{100000}$ of such compositions of $100000$, reduced modulo $1{,}111{,}124{,}111$. The solution does not enumerate compositions. Instead it builds a generating function for allowed runs, converts that to a divisor-sum coefficient sequence, and then extracts the required coefficient by inverting a formal power series. Mathematical Approach Let $u_n$ be the number of odd-run compositions of $n$, and let $$U(x)=\sum_{n\ge 0}u_nx^n,$$ with $u_0=1$ for the empty composition. The derivation becomes clean once compositions are viewed as sequences of runs instead of sequences of individual parts. Runs of a Fixed Part Size Fix a part size $j\ge 1$. A maximal run consisting of $j$'s may have length $1,3,5,\dots$, so its contribution to the ordinary generating function is $$R_j(x)=x^j+x^{3j}+x^{5j}+\cdots=\frac{x^j}{1-x^{2j}}.$$ An odd-run composition is therefore a sequence of such runs where two consecutive runs are not allowed to have the same part size, because otherwise they would merge into a longer run....

Detailed mathematical approach

Problem Summary

An odd-run composition of $n$ is a composition $n=p_1+p_2+\cdots+p_m$ in which every maximal block of equal adjacent parts has odd length. For example, $2+1+1+1$ is valid because the run of 1's has length 3, while $2+2+1$ is invalid because the run of 2's has length 2.

The problem asks for the number $u_{100000}$ of such compositions of $100000$, reduced modulo $1{,}111{,}124{,}111$. The solution does not enumerate compositions. Instead it builds a generating function for allowed runs, converts that to a divisor-sum coefficient sequence, and then extracts the required coefficient by inverting a formal power series.

Mathematical Approach

Let $u_n$ be the number of odd-run compositions of $n$, and let

$$U(x)=\sum_{n\ge 0}u_nx^n,$$

with $u_0=1$ for the empty composition. The derivation becomes clean once compositions are viewed as sequences of runs instead of sequences of individual parts.

Runs of a Fixed Part Size

Fix a part size $j\ge 1$. A maximal run consisting of $j$'s may have length $1,3,5,\dots$, so its contribution to the ordinary generating function is

$$R_j(x)=x^j+x^{3j}+x^{5j}+\cdots=\frac{x^j}{1-x^{2j}}.$$

An odd-run composition is therefore a sequence of such runs where two consecutive runs are not allowed to have the same part size, because otherwise they would merge into a longer run.

Smirnov Words and the Global Generating Function

A composition can be regarded as a word over the infinite alphabet $\{1,2,3,\dots\}$, where the letter $j$ carries weight $x^j$. After compressing each maximal run into one symbol, we get a word with no equal adjacent letters. For such words, the standard Smirnov-word generating function with letter weights $y_j$ is

$$S(\{y_j\})=\frac{1}{1-\sum_{j\ge 1}\frac{y_j}{1+y_j}}.$$

One way to see this is to start from unrestricted words, whose generating function is $1/(1-\sum z_j)$, and replace a letter $j$ by a positive run of $j$'s. Then $y_j=z_j+z_j^2+z_j^3+\cdots=z_j/(1-z_j)$, so $z_j=y_j/(1+y_j)$.

For odd-run compositions we substitute $y_j=R_j(x)$, which gives

$$U(x)=\frac{1}{1-\sum_{j\ge 1}\frac{R_j(x)}{1+R_j(x)}}=\frac{1}{1-\sum_{j\ge 1}\frac{x^j}{1+x^j-x^{2j}}}.$$

So the whole problem is reduced to understanding the series

$$A(x)=\sum_{j\ge 1}\frac{x^j}{1+x^j-x^{2j}}.$$

Why Fibonacci Numbers Appear

The Fibonacci generating function is

$$\sum_{m\ge 1}F_m t^m=\frac{t}{1-t-t^2},\qquad F_1=F_2=1.$$

Substituting $t=-x^j$ yields

$$\frac{x^j}{1+x^j-x^{2j}}=\sum_{m\ge 1}(-1)^{m-1}F_m x^{jm}.$$

Now collect equal powers of $x$. Writing

$$A(x)=\sum_{n\ge 1}a_nx^n,$$

the coefficient of $x^n$ receives one contribution from each factorization $n=jm$. Equivalently,

$$a_n=\sum_{d\mid n}(-1)^{d-1}F_d.$$

This divisor sum is exactly the sequence constructed by the implementations before any fast polynomial work begins.

The Coefficient Recurrence

Since $U(x)=1/(1-A(x))$, we have

$$\bigl(1-A(x)\bigr)U(x)=1.$$

Comparing coefficients gives the fundamental recurrence

$$u_0=1,\qquad u_n=\sum_{k=1}^{n}a_k\,u_{n-k}\quad(n\ge 1).$$

This already solves the problem conceptually: once the divisor-sum coefficients $a_k$ are known, the answer sequence is determined uniquely.

Worked Example: $n=5$

The first divisor-sum coefficients are

$$a_1=1,\quad a_2=0,\quad a_3=3,\quad a_4=-3,\quad a_5=6.$$

Using the recurrence,

$$u_1=1,\qquad u_2=1,\qquad u_3=4,\qquad u_4=4,$$

and then

$$u_5=a_1u_4+a_2u_3+a_3u_2+a_4u_1+a_5u_0=4+0+3-3+6=10.$$

This matches direct inspection. Valid examples include $5$, $1+3+1$, $2+1+2$, $2+1+1+1$, and $1+1+1+1+1$. Invalid examples such as $2+2+1$ or $1+1+3$ fail because they contain an even-length run.

How the Code Works

Building the Coefficients of $A(x)$

The C++, Python, and Java implementations first generate Fibonacci numbers modulo $1{,}111{,}124{,}111$. Then, for each $d$, they add the signed Fibonacci value $(-1)^{d-1}F_d$ to every multiple of $d$. After this divisor sweep, the array holds $a_1,a_2,\dots,a_n$, the coefficients of $A(x)$.

Recovering $U(x)=1/(1-A(x))$

A direct use of the coefficient recurrence is quadratic, so it is only suitable for small verification cases. For the real target $n=100000$, the implementations form the series $B(x)=1-A(x)$ and compute its inverse modulo $x^{n+1}$ by Newton iteration. If $G(x)$ is correct up to degree $m-1$, the next approximation is

$$G_{\text{new}}(x)=G(x)\bigl(2-B(x)G(x)\bigr)\pmod{x^m}.$$

Because the constant term of $B(x)$ is 1, this inverse exists as a formal power series. Each Newton step roughly doubles the number of correct coefficients, and the coefficient of $x^n$ in the final inverse is exactly $u_n$.

Fast Polynomial Multiplication

The expensive part of Newton iteration is polynomial multiplication. The C++ and Java implementations evaluate these convolutions with a number-theoretic transform under three auxiliary NTT-friendly prime moduli, then recombine the results with the Chinese remainder theorem to recover coefficients modulo $1{,}111{,}124{,}111$. The Python implementation is a thin wrapper around the same C++ computation, so all three language entries follow the same mathematical pipeline.

Complexity Analysis

Constructing the divisor-sum coefficients costs

$$\sum_{d=1}^{n}\left\lfloor\frac{n}{d}\right\rfloor=O(n\log n).$$

The fast phase is the formal power-series inversion. With FFT-style multiplication cost $M(n)=O(n\log n)$, Newton iteration runs in $O(M(n)\log n)$, which here is $O(n\log^2 n)$. Memory usage is $O(n)$ for the coefficient arrays and transform buffers.

The slower recurrence $u_n=\sum_{k=1}^{n}a_ku_{n-k}$ is $O(n^2)$ and appears only as a correctness check on small inputs. The actual computation for $n=100000$ relies on the quasi-linear series inversion.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=929
Compositions of integers: Wikipedia - Composition (combinatorics)
Fibonacci numbers and their generating function: Wikipedia - Fibonacci number
Formal power series: Wikipedia - Formal power series
Generating functions: Wikipedia - Generating function
Number-theoretic transform: Wikipedia - Number-theoretic transform
Chinese remainder theorem: Wikipedia - Chinese remainder theorem

Problem 929 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

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

constexpr i64 kMod = 1'111'124'111LL;

struct NttMod {
    i64 mod;
    i64 primitive_root;
};

constexpr std::array<NttMod, 3> kNttMods = {
    NttMod{998'244'353LL, 3LL},
    NttMod{1'004'535'809LL, 3LL},
    NttMod{469'762'049LL, 3LL},
};

i64 mod_pow(i64 a, i64 e, i64 mod) {
    i64 r = 1 % mod;
    a %= mod;
    while (e > 0) {
        if (e & 1LL) {
            r = static_cast<i64>((static_cast<u128>(r) * a) % mod);
        }
        a = static_cast<i64>((static_cast<u128>(a) * a) % mod);
        e >>= 1LL;
    }
    return r;
}

i64 mod_inv_prime(i64 a, i64 mod) {
    return mod_pow(a, mod - 2, mod);
}

void ntt(std::vector<i64>& a, i64 mod, i64 primitive_root, bool invert) {
    const int n = static_cast<int>(a.size());

    for (int i = 1, j = 0; i < n; ++i) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) {
            j ^= bit;
        }
        j ^= bit;
        if (i < j) {
            std::swap(a[i], a[j]);
        }
    }

    for (int len = 2; len <= n; len <<= 1) {
        const i64 wlen = mod_pow(primitive_root, (mod - 1) / len, mod);
        const i64 wlen_inv = mod_inv_prime(wlen, mod);
        const i64 step = invert ? wlen_inv : wlen;

        for (int i = 0; i < n; i += len) {
            i64 w = 1;
            for (int j = 0; j < len / 2; ++j) {
                const i64 u = a[i + j];
                const i64 v = static_cast<i64>((static_cast<u128>(a[i + j + len / 2]) * w) % mod);
                i64 x = u + v;
                if (x >= mod) {
                    x -= mod;
                }
                i64 y = u - v;
                if (y < 0) {
                    y += mod;
                }
                a[i + j] = x;
                a[i + j + len / 2] = y;
                w = static_cast<i64>((static_cast<u128>(w) * step) % mod);
            }
        }
    }

    if (invert) {
        const i64 n_inv = mod_inv_prime(n, mod);
        for (i64& x : a) {
            x = static_cast<i64>((static_cast<u128>(x) * n_inv) % mod);
        }
    }
}

std::vector<i64> convolution_single_mod(const std::vector<i64>& a,
                                        const std::vector<i64>& b,
                                        i64 mod,
                                        i64 primitive_root) {
    if (a.empty() || b.empty()) {
        return {};
    }

    if (std::min(a.size(), b.size()) <= 64) {
        std::vector<i64> out(a.size() + b.size() - 1, 0);
        for (std::size_t i = 0; i < a.size(); ++i) {
            for (std::size_t j = 0; j < b.size(); ++j) {
                out[i + j] = (out[i + j] + static_cast<i64>((static_cast<u128>(a[i]) * b[j]) % mod)) % mod;
            }
        }
        return out;
    }

    int n = 1;
    while (n < static_cast<int>(a.size() + b.size() - 1)) {
        n <<= 1;
    }

    std::vector<i64> fa(n, 0), fb(n, 0);
    for (std::size_t i = 0; i < a.size(); ++i) {
        fa[i] = a[i] % mod;
    }
    for (std::size_t i = 0; i < b.size(); ++i) {
        fb[i] = b[i] % mod;
    }

    ntt(fa, mod, primitive_root, false);
    ntt(fb, mod, primitive_root, false);

    for (int i = 0; i < n; ++i) {
        fa[i] = static_cast<i64>((static_cast<u128>(fa[i]) * fb[i]) % mod);
    }

    ntt(fa, mod, primitive_root, true);
    fa.resize(a.size() + b.size() - 1);
    return fa;
}

std::vector<i64> convolution_mod(const std::vector<i64>& a, const std::vector<i64>& b, i64 target_mod) {
    if (a.empty() || b.empty()) {
        return {};
    }

    std::array<std::vector<i64>, 3> convs;
    for (int t = 0; t < 3; ++t) {
        convs[t] = convolution_single_mod(a, b, kNttMods[t].mod, kNttMods[t].primitive_root);
    }

    const i64 p1 = kNttMods[0].mod;
    const i64 p2 = kNttMods[1].mod;
    const i64 p3 = kNttMods[2].mod;

    const i64 inv_p1_mod_p2 = mod_inv_prime(p1 % p2, p2);
    const i64 p1_mod_p3 = p1 % p3;
    const i64 p1p2_mod_p3 = static_cast<i64>((static_cast<u128>(p1_mod_p3) * (p2 % p3)) % p3);
    const i64 inv_p1p2_mod_p3 = mod_inv_prime(p1p2_mod_p3, p3);

    const i64 p1_mod_m = p1 % target_mod;
    const i64 p1p2_mod_m = static_cast<i64>((static_cast<u128>(p1_mod_m) * (p2 % target_mod)) % target_mod);

    std::vector<i64> out(convs[0].size(), 0);
    for (std::size_t i = 0; i < out.size(); ++i) {
        const i64 r1 = convs[0][i];
        const i64 r2 = convs[1][i];
        const i64 r3 = convs[2][i];

        i64 t1 = r2 - r1;
        if (t1 < 0) {
            t1 += p2;
        }
        t1 = static_cast<i64>((static_cast<u128>(t1) * inv_p1_mod_p2) % p2);

        const i64 x12_mod_p3 = (r1 + static_cast<i64>((static_cast<u128>(p1_mod_p3) * t1) % p3)) % p3;

        i64 t2 = r3 - x12_mod_p3;
        if (t2 < 0) {
            t2 += p3;
        }
        t2 = static_cast<i64>((static_cast<u128>(t2) * inv_p1p2_mod_p3) % p3);

        i64 x = r1 % target_mod;
        x += static_cast<i64>((static_cast<u128>(p1_mod_m) * (t1 % target_mod)) % target_mod);
        if (x >= target_mod) {
            x -= target_mod;
        }
        x += static_cast<i64>((static_cast<u128>(p1p2_mod_m) * (t2 % target_mod)) % target_mod);
        x %= target_mod;
        out[i] = x;
    }

    return out;
}

std::vector<i64> poly_inv_one_minus(const std::vector<i64>& b, int n) {
    std::vector<i64> inv(1, 1);

    int len = 1;
    while (len < n) {
        const int m = std::min(2 * len, n);

        std::vector<i64> b_cut(m, 0);
        for (int i = 0; i < m; ++i) {
            b_cut[i] = b[i] % kMod;
        }

        std::vector<i64> prod = convolution_mod(inv, b_cut, kMod);
        prod.resize(m, 0);

        for (int i = 0; i < m; ++i) {
            prod[i] = (kMod - prod[i]) % kMod;
        }
        prod[0] = (prod[0] + 2) % kMod;

        std::vector<i64> next = convolution_mod(inv, prod, kMod);
        next.resize(m, 0);
        inv.swap(next);

        len = m;
    }

    inv.resize(n, 0);
    return inv;
}

std::vector<i64> build_A(int n) {
    std::vector<i64> fib(n + 1, 0);
    if (n >= 1) {
        fib[1] = 1;
    }
    for (int i = 2; i <= n; ++i) {
        fib[i] = fib[i - 1] + fib[i - 2];
        if (fib[i] >= kMod) {
            fib[i] -= kMod;
        }
    }

    std::vector<i64> a(n + 1, 0);

    for (int d = 1; d <= n; ++d) {
        i64 b = fib[d];
        if ((d & 1) == 0) {
            b = (kMod - b) % kMod;
        }

        for (int m = d; m <= n; m += d) {
            a[m] += b;
            if (a[m] >= kMod) {
                a[m] -= kMod;
            }
        }
    }

    return a;
}

i64 solve_fast(int n) {
    std::vector<i64> a = build_A(n);

    std::vector<i64> one_minus_a(n + 1, 0);
    one_minus_a[0] = 1;
    for (int i = 1; i <= n; ++i) {
        one_minus_a[i] = (kMod - a[i]) % kMod;
    }

    std::vector<i64> inv = poly_inv_one_minus(one_minus_a, n + 1);
    return inv[n] % kMod;
}

i64 solve_slow(int n) {
    std::vector<i64> a = build_A(n);
    std::vector<i64> u(n + 1, 0);
    u[0] = 1;

    for (int i = 1; i <= n; ++i) {
        i64 cur = 0;
        for (int k = 1; k <= i; ++k) {
            cur += static_cast<i64>((static_cast<u128>(a[k]) * u[i - k]) % kMod);
            if (cur >= kMod) {
                cur -= kMod;
            }
        }
        u[i] = cur;
    }

    return u[n];
}

void run_validations() {
    assert(solve_slow(5) == 10);

    for (int n = 1; n <= 80; ++n) {
        assert(solve_fast(n) == solve_slow(n));
    }
}

}  // namespace

int main() {
    run_validations();
    constexpr int kN = 100'000;
    std::cout << solve_fast(kN) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.*;

public class Euler929 {

    static final long kMod = 1111124111L;

    static class NttMod {
        long mod;
        long primitiveRoot;

        NttMod(long m, long r) {
            mod = m;
            primitiveRoot = r;
        }
    }

    static final NttMod[] kNttMods = {
            new NttMod(998244353L, 3L),
            new NttMod(1004535809L, 3L),
            new NttMod(469762049L, 3L)
    };

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

    static long modInvPrime(long a, long mod) {
        return modPow(a, mod - 2, mod);
    }

    static void ntt(long[] a, long mod, long primitiveRoot, boolean invert) {
        int n = a.length;
        int j = 0;
        for (int i = 1; i < n; i++) {
            int bit = n >> 1;
            for (; (j & bit) != 0; bit >>= 1)
                j ^= bit;
            j ^= bit;
            if (i < j) {
                long temp = a[i];
                a[i] = a[j];
                a[j] = temp;
            }
        }

        for (int len = 2; len <= n; len <<= 1) {
            long wlen = modPow(primitiveRoot, (mod - 1) / len, mod);
            if (invert)
                wlen = modInvPrime(wlen, mod);

            long step = wlen;
            for (int i = 0; i < n; i += len) {
                long w = 1;
                for (int k = 0; k < len / 2; k++) {
                    long u = a[i + k];
                    long v = (a[i + k + len / 2] * w) % mod;
                    long x = u + v;
                    if (x >= mod)
                        x -= mod;
                    long y = u - v;
                    if (y < 0)
                        y += mod;
                    a[i + k] = x;
                    a[i + k + len / 2] = y;
                    w = (w * step) % mod;
                }
            }
        }

        if (invert) {
            long nInv = modInvPrime(n, mod);
            for (int i = 0; i < n; i++) {
                a[i] = (a[i] * nInv) % mod;
            }
        }
    }

    static long[] convolutionSingleMod(long[] a, long[] b, long mod, long primitiveRoot) {
        if (a.length == 0 || b.length == 0)
            return new long[0];

        if (Math.min(a.length, b.length) <= 64) {
            long[] out = new long[a.length + b.length - 1];
            for (int i = 0; i < a.length; i++) {
                for (int j = 0; j < b.length; j++) {
                    out[i + j] = (out[i + j] + (a[i] * b[j]) % mod) % mod;
                }
            }
            return out;
        }

        int n = 1;
        int targetLen = a.length + b.length - 1;
        while (n < targetLen)
            n <<= 1;

        long[] fa = new long[n];
        long[] fb = new long[n];
        for (int i = 0; i < a.length; i++)
            fa[i] = a[i] % mod;
        for (int i = 0; i < b.length; i++)
            fb[i] = b[i] % mod;

        ntt(fa, mod, primitiveRoot, false);
        ntt(fb, mod, primitiveRoot, false);

        for (int i = 0; i < n; i++) {
            fa[i] = (fa[i] * fb[i]) % mod;
        }

        ntt(fa, mod, primitiveRoot, true);
        return Arrays.copyOf(fa, targetLen);
    }

    static long[] convolutionMod(long[] a, long[] b, long targetMod) {
        if (a.length == 0 || b.length == 0)
            return new long[0];

        long[][] convs = new long[3][];
        for (int t = 0; t < 3; t++) {
            convs[t] = convolutionSingleMod(a, b, kNttMods[t].mod, kNttMods[t].primitiveRoot);
        }

        long p1 = kNttMods[0].mod;
        long p2 = kNttMods[1].mod;
        long p3 = kNttMods[2].mod;

        long invP1ModP2 = modInvPrime(p1 % p2, p2);
        long p1ModP3 = p1 % p3;
        long p1P2ModP3 = (p1ModP3 * (p2 % p3)) % p3;
        long invP1P2ModP3 = modInvPrime(p1P2ModP3, p3);

        long p1ModM = p1 % targetMod;
        long p1P2ModM = (p1ModM * (p2 % targetMod)) % targetMod;

        long[] out = new long[convs[0].length];
        for (int i = 0; i < out.length; i++) {
            long r1 = convs[0][i];
            long r2 = convs[1][i];
            long r3 = convs[2][i];

            long t1 = r2 - r1;
            if (t1 < 0)
                t1 += p2;
            t1 = (t1 * invP1ModP2) % p2;

            long x12ModP3 = (r1 + p1ModP3 * t1) % p3;

            long t2 = r3 - x12ModP3;
            if (t2 < 0)
                t2 += p3;
            t2 = (t2 * invP1P2ModP3) % p3;

            long x = r1 % targetMod;
            x = (x + p1ModM * (t1 % targetMod)) % targetMod;
            x = (x + p1P2ModM * (t2 % targetMod)) % targetMod;
            out[i] = x;
        }

        return out;
    }

    static long[] polyInvOneMinus(long[] b, int n) {
        long[] inv = { 1 };
        int len = 1;
        while (len < n) {
            int m = Math.min(2 * len, n);
            long[] bCut = new long[m];
            for (int i = 0; i < m && i < b.length; i++) {
                bCut[i] = b[i] % kMod;
            }

            long[] prod = convolutionMod(inv, bCut, kMod);
            long[] nextProd = new long[m];
            for (int i = 0; i < m && i < prod.length; i++) {
                nextProd[i] = (kMod - prod[i]) % kMod;
            }
            nextProd[0] = (nextProd[0] + 2) % kMod;

            long[] nextInv = convolutionMod(inv, nextProd, kMod);
            inv = Arrays.copyOf(nextInv, m);
            len = m;
        }
        return Arrays.copyOf(inv, n);
    }

    static long[] buildA(int n) {
        long[] fib = new long[n + 1];
        if (n >= 1)
            fib[1] = 1;
        for (int i = 2; i <= n; i++) {
            fib[i] = fib[i - 1] + fib[i - 2];
            if (fib[i] >= kMod)
                fib[i] -= kMod;
        }

        long[] a = new long[n + 1];
        for (int d = 1; d <= n; d++) {
            long b = fib[d];
            if ((d % 2) == 0) {
                b = (kMod - b) % kMod;
            }
            for (int m = d; m <= n; m += d) {
                a[m] += b;
                if (a[m] >= kMod)
                    a[m] -= kMod;
            }
        }
        return a;
    }

    static long solveFast(int n) {
        long[] a = buildA(n);
        long[] oneMinusA = new long[n + 1];
        oneMinusA[0] = 1;
        for (int i = 1; i <= n; i++) {
            oneMinusA[i] = (kMod - a[i]) % kMod;
        }

        long[] inv = polyInvOneMinus(oneMinusA, n + 1);
        return inv[n] % kMod;
    }

    public static String solve() {
        return Long.toString(solveFast(100000));
    }

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