Problem 319: Bounded Sequences

Project Euler Problem 319 Solution

EulerSolve provides an optimized solution for Project Euler Problem 319, Bounded Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We count sequences $$x_1,x_2,\dots,x_n$$ such that $$x_1=2,$$ $$x_{i-1}\lt x_i\qquad (2\le i\le n),$$ and for all $1\le i,j\le n$, $$(x_i)^j \lt (x_j+1)^i.$$ Let $t(n)$ be the number of such sequences. We are given $$t(2)=5,\qquad t(5)=293,\qquad t(10)=86195,\qquad t(20)=5227991891,$$ and we must find $$t(10^{10})\pmod{10^9}.$$ Mathematical Approach 1) The defining inequalities mean "take floors of powers of one real number". Rewrite the condition $$(x_i)^j \lt (x_j+1)^i$$ as $$x_i^{1/i} \lt (x_j+1)^{1/j} \qquad \text{for all } i,j.$$ So the intervals $$I_i=\bigl(x_i^{1/i},\ (x_i+1)^{1/i}\bigr)$$ all overlap. Therefore there exists a real number $y$ belonging to every $I_i$, and for that $y$ we have $$x_i \lt y^i \lt x_i+1,$$ hence $$x_i=\lfloor y^i\rfloor.$$ Conversely, any $y\in(2,3)$ defines a valid sequence by setting $x_i=\lfloor y^i\rfloor$. The condition $x_1=2$ is exactly $2\le y \lt 3$, and because $y\gt1$, the sequence is strictly increasing. 2) So $t(n)$ is a boundary-counting problem on $y\in[2,3)$. As $y$ moves through $[2,3)$, the sequence $$\bigl(\lfloor y\rfloor,\lfloor y^2\rfloor,\dots,\lfloor y^n\rfloor\bigr)$$ changes only when some $y^k$ crosses an integer....

Detailed mathematical approach

Problem Summary

We count sequences

$$x_1,x_2,\dots,x_n$$

such that

$$x_1=2,$$

$$x_{i-1}\lt x_i\qquad (2\le i\le n),$$

and for all $1\le i,j\le n$,

$$(x_i)^j \lt (x_j+1)^i.$$

Let $t(n)$ be the number of such sequences. We are given

$$t(2)=5,\qquad t(5)=293,\qquad t(10)=86195,\qquad t(20)=5227991891,$$

and we must find

$$t(10^{10})\pmod{10^9}.$$

Mathematical Approach

1) The defining inequalities mean "take floors of powers of one real number".

Rewrite the condition

$$(x_i)^j \lt (x_j+1)^i$$

$$x_i^{1/i} \lt (x_j+1)^{1/j} \qquad \text{for all } i,j.$$

So the intervals

$$I_i=\bigl(x_i^{1/i},\ (x_i+1)^{1/i}\bigr)$$

all overlap. Therefore there exists a real number $y$ belonging to every $I_i$, and for that $y$ we have

$$x_i \lt y^i \lt x_i+1,$$

hence

$$x_i=\lfloor y^i\rfloor.$$

Conversely, any $y\in(2,3)$ defines a valid sequence by setting $x_i=\lfloor y^i\rfloor$. The condition $x_1=2$ is exactly $2\le y \lt 3$, and because $y\gt1$, the sequence is strictly increasing.

2) So $t(n)$ is a boundary-counting problem on $y\in[2,3)$.

As $y$ moves through $[2,3)$, the sequence

$$\bigl(\lfloor y\rfloor,\lfloor y^2\rfloor,\dots,\lfloor y^n\rfloor\bigr)$$

changes only when some $y^k$ crosses an integer. The change points are exactly the numbers

$$y=m^{1/k},\qquad 1\le k\le n,\qquad 2^k \lt m \lt 3^k.$$

Therefore $t(n)$ equals

$$1+\text{(number of distinct boundary points in }(2,3)\text{ visible up to exponent }n).$$

3) Distinct roots are the real difficulty.

If we simply count all pairs $(m,k)$ with $2^k \lt m \lt 3^k$, we overcount. For example, the same boundary might be representable both as $m^{1/k}$ and as $u^{1/d}$ if $m$ is a perfect power.

To remove duplicates, define:

$$f(k)=3^k-2^k-1,$$

the number of integers strictly between $2^k$ and $3^k$, and

$$g(k)=\text{number of boundary points whose minimal exponent is exactly }k.$$

Every integer $m$ counted by $f(k)$ corresponds to a boundary whose minimal exponent divides $k$. Hence

$$f(k)=\sum_{d\mid k} g(d).$$

4) Möbius inversion produces the primitive counts.

Applying Möbius inversion gives

$$g(k)=\sum_{d\mid k}\mu(d)\,f\!\left(\frac{k}{d}\right),$$

where $\mu$ is the Möbius function.

Since each primitive boundary of degree $k$ contributes exactly one new cut point once $k\le n$, we have

$$t(n)=1+\sum_{k=1}^{n} g(k).$$

Substituting the inversion formula and swapping the order of summation yields

$$t(n)=1+\sum_{k=1}^{n} f(k)\,M\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right),$$

where

$$M(x)=\sum_{m\le x}\mu(m)$$

is the Mertens function.

This is exactly the formula implemented by the code, with

$$f(k)=3^k-2^k-1.$$

5) Worked example: why $t(2)=5$.

For $n=2$, only $k=2$ contributes new boundaries, because

$$f(1)=3-2-1=0,\qquad f(2)=9-4-1=4.$$

The four boundary points are

$$\sqrt5,\ \sqrt6,\ \sqrt7,\ \sqrt8.$$

They split the interval $[2,3)$ into five parts, producing the five sequences

$$\{2,4\},\ \{2,5\},\ \{2,6\},\ \{2,7\},\ \{2,8\}.$$

$$t(2)=5,$$

exactly as stated.

6) Prefix sums of the coefficients.

The code needs many interval sums of

$$a_k=f(k)=3^k-2^k-1.$$

Define

$$A(m)=\sum_{k=1}^{m} a_k.$$

Using geometric-series formulas,

$$A(m)=\frac{3^{m+1}-3}{2}-(2^{m+1}-2)-m.$$

Therefore any block $[l,r]$ contributes

$$\sum_{k=l}^{r}a_k=A(r)-A(l-1).$$

Because the final modulus is $10^9$, division by $2$ is not invertible modulo $10^9$. The implementation avoids trouble by first computing $3^{m+1}$ modulo $2\cdot10^9$, so the numerator is still even and can be safely halved.

7) Harmonic partition of the outer sum.

The quotient

$$q=\left\lfloor\frac{n}{k}\right\rfloor$$

is constant on long ranges $[l,r]$. So instead of summing one $k$ at a time, we group all $k$ in the same block:

$$\sum_{k=1}^{n} a_k\,M\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right) =\sum_{\text{blocks }[l,r]} \bigl(A(r)-A(l-1)\bigr)\,M(q).$$

The number of such blocks is only about $2\sqrt n$, not $n$.

8) Fast computation of the Mertens function.

For small values, $\mu$ and $M$ are precomputed with a linear sieve. For large $x$, the code uses the classical identity

$$M(x)=1-\sum_{l=2}^{x}(r-l+1)\,M\!\left(\left\lfloor\frac{x}{l}\right\rfloor\right),$$

where each interval $[l,r]$ shares the same quotient $\lfloor x/l\rfloor$. Memoization makes each large $M(x)$ value get computed only once.

Algorithm

1) Use the combinatorial reformulation

$$t(n)=1+\sum_{k=1}^{n}(3^k-2^k-1)\,M\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right).$$

2) Precompute $\mu$ and $M$ up to about $n^{2/3}$ with a linear sieve.

3) Compute large $M(x)$ recursively with quotient grouping and memoization.

4) Group the outer sum into blocks of constant $\lfloor n/k\rfloor$.

5) Use the closed form for $A(m)$ to evaluate each block in $O(1)$.

Complexity Analysis

The sieve runs up to roughly $n^{2/3}$. The outer harmonic decomposition has only

$$O(\sqrt n)$$

blocks. Combined with memoized Mertens queries, this is dramatically faster than summing $10^{10}$ terms one by one.

Checks And Final Result

The source validates

$$t(2)=5,\qquad t(5)=293,\qquad t(10)=86195,\qquad t(20)=5227991891.$$

For

$$n=10^{10},$$

the program outputs

$$268457129$$

modulo $10^9$.

Problem 319 source code

C++

#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <unordered_map>
#include <vector>

using int64 = long long;
using i128 = __int128_t;

namespace {
constexpr int64 kMod = 1000000000LL;
constexpr int64 kMod2 = 2 * kMod;

int64 mod_norm(int64 x) {
    x %= kMod;
    if (x < 0) x += kMod;
    return x;
}

int64 mul_mod(int64 a, int64 b, int64 mod) {
    return static_cast<int64>((static_cast<__int128>(a) * b) % mod);
}

int64 pow_mod(int64 base, int64 exp, int64 mod) {
    int64 res = 1 % mod;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) res = mul_mod(res, base, mod);
        base = mul_mod(base, base, mod);
        exp >>= 1;
    }
    return res;
}

int64 sum_a_prefix(int64 m) {
    if (m <= 0) return 0;
    int64 pow3 = pow_mod(3, m + 1, kMod2);
    int64 num = pow3 - 3;
    if (num < 0) num += kMod2;
    int64 sum3 = num / 2;  // (3^{m+1}-3)/2 mod 1e9

    int64 pow2 = pow_mod(2, m + 1, kMod);
    int64 sum2 = pow2 - 2;
    if (sum2 < 0) sum2 += kMod;

    int64 res = sum3 - sum2 - (m % kMod);
    return mod_norm(res);
}

class Mertens {
public:
    explicit Mertens(int64 n) {
        limit_ = static_cast<int64>(std::pow(static_cast<long double>(n), 2.0L / 3.0L)) + 1;
        if (limit_ < 1) limit_ = 1;
        init_mu();
        cache_.reserve(1 << 20);
    }

    int64 get(int64 n) {
        if (n <= limit_) return prefix_[static_cast<size_t>(n)];
        auto it = cache_.find(n);
        if (it != cache_.end()) return it->second;
        int64 res = 1;
        int64 l = 2;
        while (l <= n) {
            int64 q = n / l;
            int64 r = n / q;
            res -= (r - l + 1) * get(q);
            l = r + 1;
        }
        cache_[n] = res;
        return res;
    }

private:
    int64 limit_ = 0;
    std::vector<int> mu_;
    std::vector<int> prefix_;
    std::unordered_map<int64, int64> cache_;

    void init_mu() {
        mu_.assign(static_cast<size_t>(limit_) + 1, 0);
        prefix_.assign(static_cast<size_t>(limit_) + 1, 0);
        std::vector<int> primes;
        std::vector<int> is_comp(static_cast<size_t>(limit_) + 1, 0);
        mu_[1] = 1;
        for (int i = 2; i <= limit_; ++i) {
            if (!is_comp[static_cast<size_t>(i)]) {
                primes.push_back(i);
                mu_[static_cast<size_t>(i)] = -1;
            }
            for (int p : primes) {
                int64 v = static_cast<int64>(i) * p;
                if (v > limit_) break;
                is_comp[static_cast<size_t>(v)] = 1;
                if (i % p == 0) {
                    mu_[static_cast<size_t>(v)] = 0;
                    break;
                } else {
                    mu_[static_cast<size_t>(v)] = -mu_[static_cast<size_t>(i)];
                }
            }
        }
        for (int i = 1; i <= limit_; ++i) {
            prefix_[static_cast<size_t>(i)] = prefix_[static_cast<size_t>(i - 1)] + mu_[static_cast<size_t>(i)];
        }
    }
};

struct Segment {
    int64 l;
    int64 r;
    int64 q;
    int64 mq;
};

int64 compute_t_mod(int64 n, int threads) {
    Mertens mertens(n);
    std::vector<Segment> segs;
    segs.reserve(static_cast<size_t>(2 * std::sqrt(static_cast<long double>(n)) + 10));

    for (int64 l = 1; l <= n; ) {
        int64 q = n / l;
        int64 r = n / q;
        segs.push_back({l, r, q, 0});
        l = r + 1;
    }

    for (auto& seg : segs) {
        seg.mq = mertens.get(seg.q);
    }

    if (threads <= 0) threads = 1;
    threads = std::min<int>(threads, static_cast<int>(segs.size()));

    auto worker = [&](size_t start, size_t step) -> int64 {
        int64 local = 0;
        for (size_t i = start; i < segs.size(); i += step) {
            const auto& seg = segs[i];
            int64 sum_a = sum_a_prefix(seg.r) - sum_a_prefix(seg.l - 1);
            sum_a = mod_norm(sum_a);
            int64 mq_mod = seg.mq % kMod;
            if (mq_mod < 0) mq_mod += kMod;
            int64 contrib = static_cast<int64>((static_cast<__int128>(sum_a) * mq_mod) % kMod);
            local += contrib;
            if (local >= kMod || local <= -kMod) local %= kMod;
        }
        return mod_norm(local);
    };

    int64 total = 1 % kMod;
    if (threads == 1 || segs.size() < 1024) {
        total = mod_norm(total + worker(0, 1));
        return total;
    }

    std::vector<int64> partial(static_cast<size_t>(threads), 0);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<size_t>(threads));
    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            partial[static_cast<size_t>(t)] = worker(static_cast<size_t>(t), static_cast<size_t>(threads));
        });
    }
    for (auto& th : pool) th.join();

    for (int64 val : partial) {
        total += val;
        if (total >= kMod || total <= -kMod) total %= kMod;
    }
    return mod_norm(total);
}

int64 pow_ll(int64 base, int64 exp) {
    int64 res = 1;
    while (exp > 0) {
        if (exp & 1) res *= base;
        base *= base;
        exp >>= 1;
    }
    return res;
}

int64 t_direct(int n) {
    std::vector<int> mu(n + 1, 0);
    std::vector<int> is_comp(n + 1, 0);
    std::vector<int> primes;
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (!is_comp[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            int v = i * p;
            if (v > n) break;
            is_comp[v] = 1;
            if (i % p == 0) {
                mu[v] = 0;
                break;
            } else {
                mu[v] = -mu[i];
            }
        }
    }

    std::vector<int> pref(n + 1, 0);
    for (int i = 1; i <= n; ++i) pref[i] = pref[i - 1] + mu[i];

    i128 total = 1;
    for (int k = 1; k <= n; ++k) {
        i128 a = static_cast<i128>(pow_ll(3, k)) - pow_ll(2, k) - 1;
        total += a * pref[n / k];
    }
    return static_cast<int64>(total);
}

void run_validation() {
    struct Test {
        int n;
        int64 expected;
    } tests[] = {
        {2, 5},
        {5, 293},
        {10, 86195},
        {20, 5227991891LL},
    };

    for (const auto& test : tests) {
        int64 val = t_direct(test.n);
        if (val != test.expected) {
            std::cerr << "Validation failed for n=" << test.n << ": got " << val
                      << ", expected " << test.expected << "\n";
            std::exit(1);
        }
        int64 mod_val = compute_t_mod(test.n, 1);
        int64 exp_mod = test.expected % kMod;
        if (mod_val != exp_mod) {
            std::cerr << "Mod validation failed for n=" << test.n << ": got " << mod_val
                      << ", expected " << exp_mod << "\n";
            std::exit(1);
        }
    }
}

}  // namespace

int main() {
    run_validation();

    const int64 n = 10000000000LL;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;
    int64 ans = compute_t_mod(n, threads);
    std::cout << ans << "\n";
    return 0;
}

Python

import math
import sys

sys.setrecursionlimit(20000)

MOD = 1000000000
MOD2 = 2 * MOD

def pow_mod(base, exp, mod):
    return pow(base, exp, mod)

def mod_norm(x):
    return x % MOD

def sum_a_prefix(m):
    if m <= 0:
        return 0
    pow3 = pow_mod(3, m + 1, MOD2)
    num = (pow3 - 3) % MOD2
    sum3 = num // 2
    
    pow2 = pow_mod(2, m + 1, MOD)
    sum2 = (pow2 - 2) % MOD
    
    res = (sum3 - sum2 - (m % MOD)) % MOD
    return res

class Mertens:
    def __init__(self, n):
        self.limit = int(math.pow(n, 2.0 / 3.0)) + 1
        if self.limit < 1:
            self.limit = 1
        self.mu = [0] * (self.limit + 1)
        self.prefix = [0] * (self.limit + 1)
        self.cache = {}
        self.init_mu()
        
    def init_mu(self):
        primes = []
        is_comp = bytearray(self.limit + 1)
        self.mu[1] = 1
        for i in range(2, self.limit + 1):
            if not is_comp[i]:
                primes.append(i)
                self.mu[i] = -1
            for p in primes:
                v = i * p
                if v > self.limit:
                    break
                is_comp[v] = 1
                if i % p == 0:
                    self.mu[v] = 0
                    break
                else:
                    self.mu[v] = -self.mu[i]
        
        for i in range(1, self.limit + 1):
            self.prefix[i] = self.prefix[i - 1] + self.mu[i]
            
    def get(self, n):
        if n <= self.limit:
            return self.prefix[n]
        if n in self.cache:
            return self.cache[n]
        res = 1
        l = 2
        while l <= n:
            q = n // l
            r = n // q
            res -= (r - l + 1) * self.get(q)
            l = r + 1
        self.cache[n] = res
        return res

def compute_t_mod(n):
    mertens = Mertens(n)
    segs = []
    l = 1
    while l <= n:
        q = n // l
        r = n // q
        segs.append((l, r, q))
        l = r + 1
        
    total = 1 % MOD
    for l_val, r_val, q_val in segs:
        sum_a = (sum_a_prefix(r_val) - sum_a_prefix(l_val - 1)) % MOD
        mq_mod = mertens.get(q_val) % MOD
        contrib = (sum_a * mq_mod) % MOD
        total = (total + contrib) % MOD
        
    return total

def solve():
    n = 10000000000
    ans = compute_t_mod(n)
    return str(ans)
    
if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler319 {
    static final long MOD = 1000000000L;
    static final long MOD2 = 2 * MOD;

    static long modNorm(long x) {
        x %= MOD;
        if (x < 0)
            x += MOD;
        return x;
    }

    static long powMod(long base, long exp, long mod) {
        long res = 1 % mod;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * base) % mod;
            base = (base * base) % mod;
            exp >>= 1;
        }
        return res;
    }

    static long sumAPrefix(long m) {
        if (m <= 0)
            return 0;
        long pow3 = powMod(3, m + 1, MOD2);
        long num = pow3 - 3;
        if (num < 0)
            num += MOD2;
        long sum3 = num / 2;

        long pow2 = powMod(2, m + 1, MOD);
        long sum2 = pow2 - 2;
        if (sum2 < 0)
            sum2 += MOD;

        long res = sum3 - sum2 - (m % MOD);
        return modNorm(res);
    }

    static class Mertens {
        int limit;
        int[] mu;
        int[] prefix;
        Map<Long, Long> cache;

        Mertens(long n) {
            limit = (int) Math.pow((double) n, 2.0 / 3.0) + 1;
            if (limit < 1)
                limit = 1;
            mu = new int[limit + 1];
            prefix = new int[limit + 1];
            cache = new HashMap<>();
            initMu();
        }

        void initMu() {
            List<Integer> primes = new ArrayList<>();
            byte[] isComp = new byte[limit + 1];
            mu[1] = 1;
            for (int i = 2; i <= limit; ++i) {
                if (isComp[i] == 0) {
                    primes.add(i);
                    mu[i] = -1;
                }
                for (int p : primes) {
                    long v = (long) i * p;
                    if (v > limit)
                        break;
                    isComp[(int) v] = 1;
                    if (i % p == 0) {
                        mu[(int) v] = 0;
                        break;
                    } else {
                        mu[(int) v] = -mu[i];
                    }
                }
            }
            for (int i = 1; i <= limit; ++i) {
                prefix[i] = prefix[i - 1] + mu[i];
            }
        }

        long get(long n) {
            if (n <= limit)
                return prefix[(int) n];
            if (cache.containsKey(n))
                return cache.get(n);
            long res = 1;
            long l = 2;
            while (l <= n) {
                long q = n / l;
                long r = n / q;
                res -= (r - l + 1) * get(q);
                l = r + 1;
            }
            cache.put(n, res);
            return res;
        }
    }

    static long computeTMod(long n) {
        Mertens mertens = new Mertens(n);
        long total = 1 % MOD;
        long l = 1;
        while (l <= n) {
            long q = n / l;
            long r = n / q;
            long sumA = modNorm(sumAPrefix(r) - sumAPrefix(l - 1));
            long mqMod = modNorm(mertens.get(q));
            long contrib = (sumA * mqMod) % MOD;
            total = modNorm(total + contrib);
            l = r + 1;
        }
        return total;
    }

    public static String solve() {
        long n = 10000000000L;
        long ans = computeTMod(n);
        return String.valueOf(ans);
    }

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