Problem 873: Words with Gaps

Project Euler Problem 873 Solution

EulerSolve provides an optimized solution for Project Euler Problem 873, Words with Gaps, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We are given $p$ copies of $A$, $q$ copies of $B$, and $r$ copies of $C$. A word is valid if every $A$ and every $B$ have at least two $C$'s between them. The task is to count all valid words and return the answer modulo $10^9+7$. The brute-force space has size $\binom{p+q+r}{p,q,r}$, so direct enumeration is impossible for the real parameters. The key simplification is that the exact positions of the $C$'s matter only through how many times the word switches between $A$ and $B$ after all $C$'s are removed. Mathematical Approach Let $W(p,q,r)$ denote the number of valid words built from $p$ letters $A$, $q$ letters $B$, and $r$ letters $C$. Step 1: Remove All $C$'s and Study the $A/B$ Skeleton Set $$m=p+q.$$ If we delete every $C$ from a word, the remaining letters form a binary skeleton $$s_1s_2\dots s_m,\qquad s_i\in\{A,B\},$$ with exactly $p$ letters $A$ and $q$ letters $B$. For such a skeleton define its transition count by $$t(s)=\#\{\,i\in\{1,\dots,m-1\}: s_i\ne s_{i+1}\,\}.$$ This number $t(s)$ counts how many adjacent boundaries in the skeleton change from $A$ to $B$ or from $B$ to $A$. Step 2: Each Transition Forces Two $C$'s Suppose two consecutive letters of the skeleton are different....

Detailed mathematical approach

Problem Summary

We are given $p$ copies of $A$, $q$ copies of $B$, and $r$ copies of $C$. A word is valid if every $A$ and every $B$ have at least two $C$'s between them. The task is to count all valid words and return the answer modulo $10^9+7$.

The brute-force space has size $\binom{p+q+r}{p,q,r}$, so direct enumeration is impossible for the real parameters. The key simplification is that the exact positions of the $C$'s matter only through how many times the word switches between $A$ and $B$ after all $C$'s are removed.

Mathematical Approach

Let $W(p,q,r)$ denote the number of valid words built from $p$ letters $A$, $q$ letters $B$, and $r$ letters $C$.

Step 1: Remove All $C$'s and Study the $A/B$ Skeleton

Set

$$m=p+q.$$

If we delete every $C$ from a word, the remaining letters form a binary skeleton

$$s_1s_2\dots s_m,\qquad s_i\in\{A,B\},$$

with exactly $p$ letters $A$ and $q$ letters $B$. For such a skeleton define its transition count by

$$t(s)=\#\{\,i\in\{1,\dots,m-1\}: s_i\ne s_{i+1}\,\}.$$

This number $t(s)$ counts how many adjacent boundaries in the skeleton change from $A$ to $B$ or from $B$ to $A$.

Step 2: Each Transition Forces Two $C$'s

Suppose two consecutive letters of the skeleton are different. In the full word, those two letters become an adjacent $A/B$ pair once the $C$'s between them are ignored, so that boundary must contain at least two $C$'s. Otherwise those two letters would violate the rule immediately.

The converse is also true. If every transition boundary of the skeleton receives at least two $C$'s, then any $A$ and any $B$ in the whole word are separated by at least one such boundary, hence by at least two $C$'s in total.

Therefore a skeleton with $t$ transitions consumes exactly

$$2t$$

mandatory letters $C$. The whole constraint is reduced to a transition cost.

Step 3: Count Skeletons with a Fixed Number of Transitions

Let $N_{p,q}(t)$ be the number of $A/B$ skeletons with $p$ letters $A$, $q$ letters $B$, and exactly $t$ transitions. A skeleton with $t$ transitions has $t+1$ runs.

If $t=2u+1$ is odd, then the word starts and ends with different letters, so it has $u+1$ runs of $A$ and $u+1$ runs of $B$. Splitting $p$ into $u+1$ positive run lengths gives $\binom{p-1}{u}$ choices, and similarly for $q$. Since the skeleton may start with either $A$ or $B$,

$$N_{p,q}(2u+1)=2\binom{p-1}{u}\binom{q-1}{u}.$$

If $t=2u$ is even, then the skeleton starts and ends with the same letter. Either $A$ uses $u+1$ runs and $B$ uses $u$ runs, or the roles are reversed, so

$$N_{p,q}(2u)=\binom{p-1}{u}\binom{q-1}{u-1}+\binom{p-1}{u-1}\binom{q-1}{u}.$$

As usual, we interpret $\binom{n}{k}=0$ when $k<0$ or $k>n$.

Step 4: Distribute the Remaining $C$'s

Fix a skeleton with $t$ transitions. After reserving $2t$ mandatory $C$'s, there remain

$$r-2t$$

letters $C$ to place freely. There are $m+1$ slots: before the first skeleton letter, after the last one, and one slot between each consecutive pair of skeleton letters.

By the stars-and-bars formula, the number of ways to distribute the remaining $C$'s is

$$\binom{(r-2t)+(m+1)-1}{(m+1)-1}=\binom{r-2t+m}{m},$$

provided $2t\le r$. So every skeleton with the same transition count contributes the same amount.

Step 5: Final Formula and the Natural Bound for $t$

If one of the letters $A$ or $B$ is absent, there is no restriction at all, so

$$W(p,q,r)=\binom{r+m}{m}\qquad\text{when }p=0\text{ or }q=0.$$

Otherwise

$$W(p,q,r)=\sum_{t\ge 1} N_{p,q}(t)\binom{r-2t+m}{m},$$

where terms with $2t>r$ vanish automatically. Because the count is symmetric in $p$ and $q$, the implementation first swaps them so that $p\le q$. Then the maximum possible transition count is

$$t_{\max}=\min\left(m-1,\left\lfloor\frac{r}{2}\right\rfloor,\begin{cases} 2p-1,& p=q,\\ 2p,& p<q. \end{cases}\right).$$

The last bound is purely combinatorial: the scarcer letter cannot support more alternations than this.

Worked Example: $(p,q,r)=(2,2,4)$

Here $m=4$. Since both $A$ and $B$ are present, we sum over transition counts.

For $t=1$, the valid skeletons are $AABB$ and $BBAA$. Each needs $2$ mandatory $C$'s, leaving $2$ free $C$'s. The number of insertions is

$$\binom{4-2+4}{4}=\binom{6}{4}=15,$$

so the contribution is $2\cdot 15=30$.

For $t=2$, the valid skeletons are $ABBA$ and $BAAB$. Now all $4$ letters $C$ are forced, so the insertion factor is

$$\binom{4-4+4}{4}=\binom{4}{4}=1,$$

and the contribution is $2\cdot 1=2$.

For $t=3$, we would need $6$ letters $C$, but only $4$ are available, so this case contributes nothing. Therefore

$$W(2,2,4)=30+2=32,$$

which matches the exact checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same formula. They first swap the counts so that $p\le q$, compute $m=p+q$, and handle the easy case $p=0$ or $q=0$ by returning $\binom{r+m}{m}$ modulo $10^9+7$.

Instead of building factorial tables up to $r+m$, the implementation computes

$$\binom{r+m}{m}=\prod_{i=1}^{m}\frac{r+i}{i}\pmod{10^9+7}$$

using modular inverses of $1,2,\dots,m$. It also builds the short binomial rows $\binom{p-1}{u}$ and $\binom{q-1}{u}$ incrementally, which is enough to evaluate the odd and even formulas for $N_{p,q}(t)$.

The second important recurrence updates the insertion factor without recomputing a fresh binomial coefficient each time:

$$\binom{n-2}{m}=\binom{n}{m}\cdot\frac{(n-m)(n-m-1)}{n(n-1)}.$$

Starting from $n=r+m$, one loop step moves from transition count $t-1$ to $t$. To make the divisions legal modulo $10^9+7$, the implementation precomputes inverses for the consecutive denominators that appear in this recurrence.

Finally, it loops over $t=1,\dots,t_{\max}$, evaluates the appropriate run-count formula for $N_{p,q}(t)$, multiplies by $\binom{r-2t+m}{m}$, and accumulates the result modulo $10^9+7$.

Complexity Analysis

After the swap, let $m=p+q$ and $t_{\max}$ be the transition limit above. Precomputing inverses up to $m$ costs $O(m)$ time and memory. The short binomial tables cost $O(t_{\max})$, and the final summation over all transition counts also costs $O(t_{\max})$.

Hence the total complexity is

$$O(m+t_{\max})=O(p+q)$$

time and

$$O(m+t_{\max})=O(p+q)$$

memory. The large value of $r$ does not create a table of size $r$; it appears only inside modular products.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=873
Stars and bars: Wikipedia — Stars and bars
Compositions of integers: Wikipedia — Composition (combinatorics)
Binomial coefficient: Wikipedia — Binomial coefficient
Modular arithmetic: Wikipedia — Modular arithmetic

Problem 873 source code

C++

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

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

static constexpr u64 MOD = 1'000'000'007ULL;

static u64 mod_pow(u64 a, u64 e) {
    u64 r = 1;
    while (e > 0) {
        if (e & 1ULL) r = (r * a) % MOD;
        a = (a * a) % MOD;
        e >>= 1ULL;
    }
    return r;
}

static u64 choose_u64(int n, int k) {
    if (k < 0 || k > n) return 0;
    k = std::min(k, n - k);
    unsigned __int128 res = 1;
    for (int i = 1; i <= k; ++i) {
        res = (res * static_cast<unsigned __int128>(n - k + i)) / static_cast<unsigned __int128>(i);
    }
    return static_cast<u64>(res);
}

static bool valid_word(const std::vector<char>& w) {
    int n = static_cast<int>(w.size());
    std::vector<int> pref(n + 1, 0);
    std::vector<int> pos_a;
    std::vector<int> pos_b;
    pos_a.reserve(n);
    pos_b.reserve(n);

    for (int i = 0; i < n; ++i) {
        pref[i + 1] = pref[i] + (w[i] == 'C' ? 1 : 0);
        if (w[i] == 'A') pos_a.push_back(i);
        if (w[i] == 'B') pos_b.push_back(i);
    }

    for (int a : pos_a) {
        for (int b : pos_b) {
            int l = std::min(a, b);
            int r = std::max(a, b);
            int cs = pref[r] - pref[l + 1];
            if (cs < 2) return false;
        }
    }
    return true;
}

static u64 brute_count(int p, int q, int r) {
    std::vector<char> w;
    w.reserve(p + q + r);
    u64 cnt = 0;

    std::function<void(int, int, int)> dfs = [&](int a, int b, int c) {
        if (a == 0 && b == 0 && c == 0) {
            if (valid_word(w)) ++cnt;
            return;
        }
        if (a > 0) {
            w.push_back('A');
            dfs(a - 1, b, c);
            w.pop_back();
        }
        if (b > 0) {
            w.push_back('B');
            dfs(a, b - 1, c);
            w.pop_back();
        }
        if (c > 0) {
            w.push_back('C');
            dfs(a, b, c - 1);
            w.pop_back();
        }
    };

    dfs(p, q, r);
    return cnt;
}

static u64 words_exact(int p, int q, int r) {
    int m = p + q;
    if (p == 0 || q == 0) return choose_u64(m + r, m);

    u64 ans = 0;
    for (int t = 1; t <= m - 1 && 2 * t <= r; ++t) {
        u64 ways_ab = 0;
        if (t & 1) {
            int u = t / 2;
            ways_ab = 2ULL * choose_u64(p - 1, u) * choose_u64(q - 1, u);
        } else {
            int u = t / 2;
            ways_ab = choose_u64(p - 1, u) * choose_u64(q - 1, u - 1)
                    + choose_u64(p - 1, u - 1) * choose_u64(q - 1, u);
        }
        if (ways_ab == 0) continue;
        ans += ways_ab * choose_u64(r - 2 * t + m, m);
    }
    return ans;
}

static u64 words_mod(i64 p, i64 q, i64 r) {
    if (p > q) std::swap(p, q);
    i64 m = p + q;

    std::vector<u64> inv(static_cast<std::size_t>(m + 1), 0);
    if (m >= 1) inv[1] = 1;
    for (i64 i = 2; i <= m; ++i) {
        inv[static_cast<std::size_t>(i)] = MOD - ((MOD / static_cast<u64>(i)) * inv[static_cast<std::size_t>(MOD % static_cast<u64>(i))]) % MOD;
    }

    u64 comb0 = 1;
    for (i64 i = 1; i <= m; ++i) {
        comb0 = (comb0 * ((r + i) % MOD)) % MOD;
        comb0 = (comb0 * inv[static_cast<std::size_t>(i)]) % MOD;
    }

    if (p == 0 || q == 0) return comb0;

    i64 max_trans = (p == q) ? (2 * p - 1) : (2 * p);
    i64 tmax = std::min({m - 1, r / 2, max_trans});
    if (tmax <= 0) return 0;

    i64 umax = tmax / 2;
    std::vector<u64> cp(static_cast<std::size_t>(umax + 1), 0);
    std::vector<u64> cq(static_cast<std::size_t>(umax + 1), 0);
    cp[0] = 1;
    cq[0] = 1;

    for (i64 u = 1; u <= umax; ++u) {
        if (u <= p - 1) {
            cp[static_cast<std::size_t>(u)] = cp[static_cast<std::size_t>(u - 1)] * ((p - u) % MOD) % MOD * inv[static_cast<std::size_t>(u)] % MOD;
        }
        if (u <= q - 1) {
            cq[static_cast<std::size_t>(u)] = cq[static_cast<std::size_t>(u - 1)] * ((q - u) % MOD) % MOD * inv[static_cast<std::size_t>(u)] % MOD;
        }
    }

    i64 n0 = r + m;
    i64 n_end = n0 - 2 * tmax;
    i64 low = n_end + 1;
    i64 high = n0;
    i64 len = high - low + 1;

    std::vector<u64> pref(static_cast<std::size_t>(len + 1), 1);
    for (i64 i = 0; i < len; ++i) {
        pref[static_cast<std::size_t>(i + 1)] = pref[static_cast<std::size_t>(i)] * ((low + i) % MOD) % MOD;
    }

    u64 inv_total = mod_pow(pref[static_cast<std::size_t>(len)], MOD - 2);
    std::vector<u64> inv_range(static_cast<std::size_t>(len), 0);
    for (i64 i = len - 1; i >= 0; --i) {
        u64 x = static_cast<u64>(low + i);
        inv_range[static_cast<std::size_t>(i)] = pref[static_cast<std::size_t>(i)] * inv_total % MOD;
        inv_total = inv_total * (x % MOD) % MOD;
    }

    auto inv_large = [&](i64 x) -> u64 {
        return inv_range[static_cast<std::size_t>(x - low)];
    };

    u64 ans = 0;
    u64 comb_t = comb0;
    i64 n_cur = n0;

    for (i64 t = 1; t <= tmax; ++t) {
        i64 n_prev = n_cur;
        comb_t = comb_t * ((n_prev - m) % MOD) % MOD;
        comb_t = comb_t * ((n_prev - m - 1) % MOD) % MOD;
        comb_t = comb_t * inv_large(n_prev) % MOD;
        comb_t = comb_t * inv_large(n_prev - 1) % MOD;
        n_cur = n_prev - 2;

        u64 ways_ab = 0;
        if (t & 1) {
            i64 u = t / 2;
            ways_ab = 2ULL * cp[static_cast<std::size_t>(u)] % MOD * cq[static_cast<std::size_t>(u)] % MOD;
        } else {
            i64 u = t / 2;
            ways_ab = (cp[static_cast<std::size_t>(u)] * cq[static_cast<std::size_t>(u - 1)]
                     + cp[static_cast<std::size_t>(u - 1)] * cq[static_cast<std::size_t>(u)]) % MOD;
        }

        ans += ways_ab * comb_t % MOD;
        if (ans >= MOD) ans -= MOD;
    }

    return ans;
}

int main() {
    assert(words_exact(2, 2, 4) == 32ULL);
    assert(words_exact(4, 4, 44) == 13'908'607'644ULL);

    for (int p = 0; p <= 3; ++p) {
        for (int q = 0; q <= 3; ++q) {
            for (int r = 0; r <= 5; ++r) {
                if (p + q + r > 9) continue;
                u64 brute = brute_count(p, q, r);
                u64 exact = words_exact(p, q, r);
                assert(brute == exact);
                assert(words_mod(p, q, r) == (exact % MOD));
            }
        }
    }

    std::cout << words_mod(1'000'000, 10'000'000, 100'000'000) << '\n';
    return 0;
}

Python

MOD = 1000000007

def mod_pow(a, e):
    return pow(a, e, MOD)

def words_mod(p, q, r):
    if p > q:
        p, q = q, p
    m = p + q

    inv = [0] * (m + 1)
    if m >= 1:
        inv[1] = 1
    for i in range(2, m + 1):
        inv[i] = MOD - ((MOD // i) * inv[MOD % i]) % MOD

    comb0 = 1
    for i in range(1, m + 1):
        comb0 = (comb0 * ((r + i) % MOD)) % MOD
        comb0 = (comb0 * inv[i]) % MOD

    if p == 0 or q == 0:
        return comb0

    max_trans = (2 * p - 1) if p == q else (2 * p)
    tmax = min(m - 1, r // 2, max_trans)
    if tmax <= 0:
        return 0

    umax = tmax // 2
    cp = [0] * (umax + 1)
    cq = [0] * (umax + 1)
    cp[0] = 1
    cq[0] = 1

    for u in range(1, umax + 1):
        if u <= p - 1:
            cp[u] = cp[u - 1] * ((p - u) % MOD) % MOD * inv[u] % MOD
        if u <= q - 1:
            cq[u] = cq[u - 1] * ((q - u) % MOD) % MOD * inv[u] % MOD

    n0 = r + m
    n_end = n0 - 2 * tmax
    low = n_end + 1
    high = n0
    length = high - low + 1

    pref = [1] * (length + 1)
    for i in range(length):
        pref[i + 1] = pref[i] * ((low + i) % MOD) % MOD

    inv_total = mod_pow(pref[length], MOD - 2)
    inv_range = [0] * length
    for i in range(length - 1, -1, -1):
        x = low + i
        inv_range[i] = pref[i] * inv_total % MOD
        inv_total = inv_total * (x % MOD) % MOD

    def inv_large(x):
        return inv_range[x - low]

    ans = 0
    comb_t = comb0
    n_cur = n0

    for t in range(1, tmax + 1):
        n_prev = n_cur
        comb_t = comb_t * ((n_prev - m) % MOD) % MOD
        comb_t = comb_t * ((n_prev - m - 1) % MOD) % MOD
        comb_t = comb_t * inv_large(n_prev) % MOD
        comb_t = comb_t * inv_large(n_prev - 1) % MOD
        n_cur = n_prev - 2

        ways_ab = 0
        if (t & 1) == 1:
            u = t // 2
            ways_ab = (2 * cp[u] * cq[u]) % MOD
        else:
            u = t // 2
            ways_ab = (cp[u] * cq[u - 1] + cp[u - 1] * cq[u]) % MOD

        ans = (ans + ways_ab * comb_t) % MOD

    return ans

def solve():
    return str(words_mod(1000000, 10000000, 100000000))

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

Java

public class Euler873 {

    static final long MOD = 1000000007L;

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

    static long wordsMod(long p, long q, long r) {
        if (p > q) {
            long temp = p;
            p = q;
            q = temp;
        }
        long m = p + q;

        long[] inv = new long[(int) (m + 1)];
        if (m >= 1)
            inv[1] = 1;
        for (int i = 2; i <= (int) m; ++i) {
            inv[i] = MOD - ((MOD / i) * inv[(int) (MOD % i)]) % MOD;
        }

        long comb0 = 1;
        for (int i = 1; i <= (int) m; ++i) {
            comb0 = (comb0 * ((r + i) % MOD)) % MOD;
            comb0 = (comb0 * inv[i]) % MOD;
        }

        if (p == 0 || q == 0)
            return comb0;

        long maxTrans = (p == q) ? (2 * p - 1) : (2 * p);
        long tmax = Math.min(m - 1, Math.min(r / 2, maxTrans));
        if (tmax <= 0)
            return 0;

        int umax = (int) (tmax / 2);
        long[] cp = new long[umax + 1];
        long[] cq = new long[umax + 1];
        cp[0] = 1;
        cq[0] = 1;

        for (int u = 1; u <= umax; ++u) {
            if (u <= p - 1) {
                cp[u] = cp[u - 1] * ((p - u) % MOD) % MOD * inv[u] % MOD;
            }
            if (u <= q - 1) {
                cq[u] = cq[u - 1] * ((q - u) % MOD) % MOD * inv[u] % MOD;
            }
        }

        long n0 = r + m;
        long nEnd = n0 - 2 * tmax;
        long low = nEnd + 1;
        long high = n0;
        int len = (int) (high - low + 1);

        long[] pref = new long[len + 1];
        pref[0] = 1;
        for (int i = 0; i < len; ++i) {
            pref[i + 1] = pref[i] * ((low + i) % MOD) % MOD;
        }

        long invTotal = modPow(pref[len], MOD - 2);
        long[] invRange = new long[len];
        for (int i = len - 1; i >= 0; --i) {
            long x = low + i;
            invRange[i] = pref[i] * invTotal % MOD;
            invTotal = invTotal * (x % MOD) % MOD;
        }

        long ans = 0;
        long combT = comb0;
        long nCur = n0;

        for (long t = 1; t <= tmax; ++t) {
            long nPrev = nCur;
            combT = combT * ((nPrev - m) % MOD) % MOD;
            combT = combT * ((nPrev - m - 1) % MOD) % MOD;
            combT = combT * invRange[(int) (nPrev - low)] % MOD;
            combT = combT * invRange[(int) (nPrev - 1 - low)] % MOD;
            nCur = nPrev - 2;

            long waysAb = 0;
            if ((t & 1) == 1) {
                int u = (int) (t / 2);
                waysAb = 2L * cp[u] % MOD * cq[u] % MOD;
            } else {
                int u = (int) (t / 2);
                waysAb = (cp[u] * cq[u - 1] + cp[u - 1] * cq[u]) % MOD;
            }

            ans = (ans + waysAb * combT) % MOD;
        }

        return ans;
    }

    public static String solve() {
        return Long.toString(wordsMod(1000000L, 10000000L, 100000000L));
    }

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