Problem 844: $k$-Markov Numbers

View on Project Euler

Project Euler Problem 844 Solution

EulerSolve provides an optimized solution for Project Euler Problem 844, $k$-Markov Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem studies the \(k\)-variable Markov-type equation $$x_1^2+x_2^2+\cdots+x_k^2=kx_1x_2\cdots x_k.$$ A Vieta-style move chooses one coordinate and replaces it by $$x_i'=k\prod_{j\ne i}x_j-x_i,$$ which preserves the equation. Starting from \((1,1,\dots,1)\), let \(\mathcal{V}_k(N)\) be the set of distinct integers \(x\le N\) that appear in at least one reachable positive tuple, and define $$M_k(N)=\sum_{x\in \mathcal{V}_k(N)}x.$$ The required quantity is $$S(K,N)=\sum_{k=3}^{K} M_k(N)\pmod{1405695061}.$$ The key difficulty is that small \(k\) produce a branching state graph, while large \(k\) collapse into a short deterministic chain. Mathematical Approach The implementation uses a hybrid argument: exact breadth-first exploration while genuine branching is still possible, and a closed-form tail once the bound \(N\) rules out tuples with three nontrivial entries. Step 1: The Vieta move preserves the Diophantine surface Fix all coordinates except \(x_i\). Then the defining equation becomes a quadratic in one variable: $$X^2-k\left(\prod_{j\ne i}x_j\right)X+\sum_{j\ne i}x_j^2=0.$$ If one root is \(X=x_i\), Vieta's formulas show that the other root is $$X'=k\prod_{j\ne i}x_j-x_i.$$ So replacing \(x_i\) by \(x_i'\) keeps the tuple on the same equation. The start tuple \((1,\dots,1)\) is valid because both sides equal \(k\)....

Detailed mathematical approach

Problem Summary

The problem studies the \(k\)-variable Markov-type equation

$$x_1^2+x_2^2+\cdots+x_k^2=kx_1x_2\cdots x_k.$$

A Vieta-style move chooses one coordinate and replaces it by

$$x_i'=k\prod_{j\ne i}x_j-x_i,$$

which preserves the equation. Starting from \((1,1,\dots,1)\), let \(\mathcal{V}_k(N)\) be the set of distinct integers \(x\le N\) that appear in at least one reachable positive tuple, and define

$$M_k(N)=\sum_{x\in \mathcal{V}_k(N)}x.$$

The required quantity is

$$S(K,N)=\sum_{k=3}^{K} M_k(N)\pmod{1405695061}.$$

The key difficulty is that small \(k\) produce a branching state graph, while large \(k\) collapse into a short deterministic chain.

Mathematical Approach

The implementation uses a hybrid argument: exact breadth-first exploration while genuine branching is still possible, and a closed-form tail once the bound \(N\) rules out tuples with three nontrivial entries.

Step 1: The Vieta move preserves the Diophantine surface

Fix all coordinates except \(x_i\). Then the defining equation becomes a quadratic in one variable:

$$X^2-k\left(\prod_{j\ne i}x_j\right)X+\sum_{j\ne i}x_j^2=0.$$

If one root is \(X=x_i\), Vieta's formulas show that the other root is

$$X'=k\prod_{j\ne i}x_j-x_i.$$

So replacing \(x_i\) by \(x_i'\) keeps the tuple on the same equation. The start tuple \((1,\dots,1)\) is valid because both sides equal \(k\). Every tuple examined by the exact search is generated from this seed by repeated applications of the same invariant move.

Step 2: Compress states by hiding all coordinates equal to \(1\)

The equation and the move are symmetric in the coordinates, so order does not matter. The implementation therefore stores only the entries greater than \(1\), sorted increasingly. A stored state

$$a_1\le a_2\le \cdots \le a_t,\qquad a_r>1,$$

represents the full \(k\)-tuple with \(k-t\) hidden ones. This compression is exact because the empty product is \(1\), and the move formulas become:

If a hidden \(1\) is replaced, the new value is

$$y=k\prod_{r=1}^{t}a_r-1.$$

If the stored value \(a_i\) is replaced, the new value is

$$y=k\prod_{r\ne i}a_r-a_i.$$

If \(y=1\), that coordinate simply disappears from the stored state. If \(y>1\), it is inserted back and the state is sorted again. This turns a \(k\)-tuple problem into a search over small sorted multisets.

Step 3: Exact exploration for the genuinely branching range

For fixed \(k\), a breadth-first search starts from the empty stored state, which represents \((1,\dots,1)\). Each reachable state with all entries at most \(N\) is visited once, and every distinct value appearing in any visited state is inserted into \(\mathcal{V}_k(N)\). Because states are stored canonically as sorted multisets, permuting coordinates does not create duplicates.

This phase is exact: it does not rely on heuristics, only on the move formulas above and on the bound \(x\le N\). The sum \(M_k(N)\) is therefore the sum of distinct reachable values, not the sum over states.

Step 4: Find the threshold where three nontrivial coordinates become impossible

Let

$$a_0=1,\qquad a_1=k-1.$$

The smallest possible state with one stored entry is therefore \((a_1)\). Replacing one more hidden \(1\) gives the smallest possible two-entry state, whose new value is

$$a_2=ka_1-a_0=k^2-k-1.$$

Now ask for the smallest possible state with three stored entries. The minimal way to create it is to start from the minimal two-entry state \((a_1,a_2)\) and replace yet another hidden \(1\). That produces

$$b_3=ka_1a_2-1=k(k-1)(k^2-k-1)-1.$$

Hence if \(b_3>N\), then no reachable tuple contributing to \(M_k(N)\) can contain three coordinates greater than \(1\). This gives the exact cutoff

$$\kappa_4=\max\left\{k\ge 3: k(k-1)(k^2-k-1)-1\le N\right\}.$$

The implementation explores \(k\le \kappa_4\) exactly and handles \(k>\kappa_4\) by formulas.

Step 5: Above the cutoff, the state graph collapses to a linear chain

Once three stored entries are impossible, every relevant state has at most two coordinates greater than \(1\). In that regime the reachable values follow the recurrence

$$a_{n+1}=ka_n-a_{n-1},\qquad a_0=1,\qquad a_1=k-1.$$

The first terms are

$$a_2=k^2-k-1,\qquad a_3=k^3-k^2-2k+1.$$

For \(k\ge 3\), the sequence is strictly increasing: if \(a_n>a_{n-1}\), then

$$a_{n+1}=ka_n-a_{n-1}>2a_n-a_{n-1}>a_n.$$

The next term is

$$a_4=ka_3-a_2=k^4-k^3-3k^2+2k+1,$$

and it satisfies

$$a_4-b_3=(k-2)(k^2-k-1)>0.$$

So if \(k>\kappa_4\), then \(b_3>N\), hence \(a_4>N\) and every later chain term is also \(>N\). Therefore only \(1\), \(a_1\), \(a_2\), and \(a_3\) can contribute in the large-\(k\) range.

Define the two remaining thresholds

$$\kappa_2=\max\left\{k\ge 3: k^2-k-1\le N\right\},$$

$$\kappa_3=\max\left\{k\ge 3: k^3-k^2-2k+1\le N\right\}.$$

Then for every \(k>\kappa_4\),

$$M_k(N)=k+\begin{cases} k^2-k-1,&k\le \kappa_2,\\ 0,&k>\kappa_2, \end{cases} +\begin{cases} k^3-k^2-2k+1,&k\le \kappa_3,\\ 0,&k>\kappa_3. \end{cases}$$

The term \(k\) is simply \(1+(k-1)\), the guaranteed contribution of the seed value and the first nontrivial value.

Worked Example: \(k=4\) and \(N=100\)

For \(k=4\), the chain begins with

$$a_0=1,\qquad a_1=3,\qquad a_2=11,\qquad a_3=41.$$

The smallest possible third stored entry is

$$b_3=4\cdot 3\cdot 11-1=131>100.$$

So no tuple with three entries greater than \(1\) can contribute. The only reachable values at or below \(100\) are therefore

$$1,\ 3,\ 11,\ 41,$$

which gives

$$M_4(100)=1+3+11+41=56.$$

Together with the \(k=3\) contribution, the full check value is

$$S(4,100)=229,$$

matching the program's built-in sanity test.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First, they compute the quartic cutoff \(\kappa_4\) and run an exact breadth-first search separately for each \(k\) from \(3\) up to \(\min(K,\kappa_4)\). A queue stores canonical states, a visited set prevents repeats, and a second set records the distinct integers that have appeared so they are counted only once in \(M_k(N)\).

During each transition, the implementation multiplies only up to the cap needed to decide whether the next value could still be at most \(N\). That prevents unnecessary big intermediate products and prunes branches as soon as they are certainly too large. When several equal entries appear in a state, only one representative replacement needs to be tried, because replacing identical coordinates yields the same canonical successor.

After the exact range, the implementation switches to the chain formulas. Three binary searches locate \(\kappa_2\), \(\kappa_3\), and \(\kappa_4\). The remaining contribution is then written as a sum of piecewise polynomials in \(k\), so it can be accumulated from closed forms for

$$\sum k,\qquad \sum k^2,\qquad \sum k^3$$

modulo \(1405695061\). Modular inverses of \(2\) and \(6\) are used to evaluate the usual prefix-sum formulas inside the modulus.

Complexity Analysis

Let \(T_k(N)\) be the number of canonical states visited in the exact search for a fixed \(k\). The exact phase costs

$$O\left(\sum_{k=3}^{\min(K,\kappa_4)} T_k(N)\cdot c_k\right),$$

where \(c_k\) is the average work per state for capped products, successor construction, and set operations. The memory usage of that phase is \(O(T_k(N))\) per active \(k\), dominated by the visited states and the set of distinct numbers.

The important structural fact is that \(\kappa_4\) is defined by a quartic inequality, so the exact outer loop covers only \(O(N^{1/4})\) values of \(k\). The large-\(k\) tail needs only three binary searches, hence \(O(\log N)\) time, plus \(O(1)\) arithmetic for the closed-form interval sums. In practice the runtime is dominated by the exact searches near the cutoff, and the tail is negligible.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=844
  2. Markov numbers and the classical three-variable case: Wikipedia — Markov number
  3. Vieta jumping and root replacement arguments: Wikipedia — Vieta jumping
  4. Breadth-first search: Wikipedia — Breadth-first search
  5. Linear recurrences: Wikipedia — Recurrence relation

Problem 844 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <deque>
#include <iostream>
#include <unordered_set>
#include <vector>

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

static constexpr u64 kMod = 1'405'695'061ULL;

static inline u64 mod_add(u64 a, u64 b) {
    a += b;
    if (a >= kMod) a -= kMod;
    return a;
}

static inline u64 mod_sub(u64 a, u64 b) {
    return (a >= b) ? (a - b) : (a + kMod - b);
}

static inline u64 mod_mul(u64 a, u64 b) {
    return static_cast<u64>((static_cast<u128>(a) * b) % kMod);
}

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

static u64 sum1_prefix(u64 n) {
    static const u64 inv2 = mod_pow(2, kMod - 2);
    const u64 a = n % kMod;
    const u64 b = (n + 1) % kMod;
    return mod_mul(mod_mul(a, b), inv2);
}

static u64 sum2_prefix(u64 n) {
    static const u64 inv6 = mod_pow(6, kMod - 2);
    const u64 a = n % kMod;
    const u64 b = (n + 1) % kMod;
    const u64 c = (2 * (n % kMod) + 1) % kMod;
    return mod_mul(mod_mul(mod_mul(a, b), c), inv6);
}

static u64 sum3_prefix(u64 n) {
    const u64 s1 = sum1_prefix(n);
    return mod_mul(s1, s1);
}

static u64 range_sum1(u64 l, u64 r) {
    if (l > r) return 0;
    return mod_sub(sum1_prefix(r), sum1_prefix(l - 1));
}

static u64 range_sum2(u64 l, u64 r) {
    if (l > r) return 0;
    return mod_sub(sum2_prefix(r), sum2_prefix(l - 1));
}

static u64 range_sum3(u64 l, u64 r) {
    if (l > r) return 0;
    return mod_sub(sum3_prefix(r), sum3_prefix(l - 1));
}

static u64 u2_value(u64 k) {
    return k * k - k - 1;
}

static u128 u3_value_u128(u64 k) {
    return static_cast<u128>(k) * k * k - static_cast<u128>(k) * k - 2 * static_cast<u128>(k) + 1;
}

static u128 p4_seed_u128(u64 k) {
    const u128 a = static_cast<u128>(k);
    return a * (k - 1) * (static_cast<u128>(k) * k - k - 1) - 1;
}

static u64 max_k_u2(u64 n) {
    u64 lo = 0, hi = 2'000'000'000ULL;
    while (lo < hi) {
        u64 mid = lo + (hi - lo + 1) / 2;
        if (u2_value(mid) <= n) lo = mid;
        else hi = mid - 1;
    }
    return lo;
}

static u64 max_k_u3(u64 n) {
    u64 lo = 0, hi = 2'000'000ULL;
    while (lo < hi) {
        u64 mid = lo + (hi - lo + 1) / 2;
        if (u3_value_u128(mid) <= static_cast<u128>(n)) lo = mid;
        else hi = mid - 1;
    }
    return lo;
}

static u64 max_k_seed4(u64 n) {
    u64 lo = 0, hi = 2'000'000ULL;
    while (lo < hi) {
        u64 mid = lo + (hi - lo + 1) / 2;
        if (p4_seed_u128(mid) <= static_cast<u128>(n)) lo = mid;
        else hi = mid - 1;
    }
    return lo;
}

struct VecHash {
    std::size_t operator()(const std::vector<u64>& v) const noexcept {
        std::size_t h = v.size();
        for (u64 x : v) {
            std::size_t y = std::hash<u64>{}(x);
            h ^= y + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        }
        return h;
    }
};

static bool capped_product_except(const std::vector<u64>& v, int skip, u64 cap, u64& out) {
    u128 prod = 1;
    for (int i = 0; i < static_cast<int>(v.size()); ++i) {
        if (i == skip) continue;
        if (prod > static_cast<u128>(cap) / v[i]) return false;
        prod *= v[i];
    }
    out = static_cast<u64>(prod);
    return true;
}

static u64 exact_M_k(u64 k, u64 n) {
    std::deque<std::vector<u64>> q;
    std::unordered_set<std::vector<u64>, VecHash> seen;
    std::unordered_set<u64> numbers;

    q.push_back({});
    seen.insert({});
    numbers.insert(1);

    while (!q.empty()) {
        std::vector<u64> cur = q.front();
        q.pop_front();

        for (u64 x : cur) {
            if (x <= n) numbers.insert(x);
        }

        const u64 ones = k - static_cast<u64>(cur.size());

        if (ones > 0) {
            const u64 cap = (n + 1) / k;
            u64 prod = 1;
            if (capped_product_except(cur, -1, cap, prod)) {
                u64 y = k * prod - 1;
                if (y <= n) {
                    std::vector<u64> nxt = cur;
                    nxt.push_back(y);
                    auto it = std::upper_bound(nxt.begin(), nxt.end() - 1, y);
                    std::rotate(it, nxt.end() - 1, nxt.end());
                    if (seen.insert(nxt).second) q.push_back(std::move(nxt));
                }
            }
        }

        for (int i = 0; i < static_cast<int>(cur.size()); ++i) {
            if (i > 0 && cur[i] == cur[i - 1]) continue;
            u64 x = cur[i];
            const u64 cap = (n + x) / k;
            u64 prod_others = 1;
            if (!capped_product_except(cur, i, cap, prod_others)) continue;

            u64 y = k * prod_others - x;
            if (y == 0) continue;

            std::vector<u64> nxt;
            nxt.reserve(cur.size());
            for (int j = 0; j < static_cast<int>(cur.size()); ++j) {
                if (j != i) nxt.push_back(cur[j]);
            }
            if (y > 1) {
                auto it = std::upper_bound(nxt.begin(), nxt.end(), y);
                nxt.insert(it, y);
            }
            if (!nxt.empty() && nxt.back() > n) continue;
            if (seen.insert(nxt).second) q.push_back(std::move(nxt));
        }
    }

    u64 out = 0;
    for (u64 x : numbers) out = (out + (x % kMod)) % kMod;
    return out;
}

static u64 chain_sum_range(u64 l, u64 r, u64 n, u64 k2max, u64 k3max) {
    if (l > r) return 0;

    u64 ans = range_sum1(l, r);

    if (l <= k2max) {
        u64 rr = std::min(r, k2max);
        u64 cnt = (rr - l + 1) % kMod;
        u64 s2 = range_sum2(l, rr);
        u64 s1 = range_sum1(l, rr);
        u64 add = mod_sub(mod_sub(s2, s1), cnt);
        ans = mod_add(ans, add);
    }

    if (l <= k3max) {
        u64 rr = std::min(r, k3max);
        u64 cnt = (rr - l + 1) % kMod;
        u64 s3 = range_sum3(l, rr);
        u64 s2 = range_sum2(l, rr);
        u64 s1 = range_sum1(l, rr);
        u64 add = mod_sub(mod_sub(s3, s2), mod_mul(2, s1));
        add = mod_add(add, cnt);
        ans = mod_add(ans, add);
    }

    (void)n;
    return ans;
}

static u64 solve(u64 K, u64 N) {
    if (K < 3) return 0;

    const u64 k2max = max_k_u2(N);
    const u64 k3max = max_k_u3(N);
    const u64 k0 = max_k_seed4(N);

    const u64 exact_to = std::min(K, k0);
    u64 ans = 0;
    for (u64 k = 3; k <= exact_to; ++k) {
        ans = mod_add(ans, exact_M_k(k, N));
    }

    if (K > exact_to) {
        const u64 l = std::max<u64>(3, exact_to + 1);
        ans = mod_add(ans, chain_sum_range(l, K, N, k2max, k3max));
    }

    return ans;
}

int main() {
    assert(solve(4, 100) == 229U);
    assert(solve(10, 100'000'000ULL) == 2'383'369'980ULL % kMod);
    assert(exact_M_k(8, 100'000'000ULL) == 131'493'335ULL % kMod);

    std::cout << solve(1'000'000'000'000'000'000ULL, 1'000'000'000'000'000'000ULL) << '\n';
    return 0;
}

Python

from collections import deque

kMod = 1405695061

def mod_add(a, b):
    a += b
    if a >= kMod:
        a -= kMod
    return a

def mod_sub(a, b):
    return (a - b) if a >= b else (a + kMod - b)

def mod_mul(a, b):
    return (a * b) % kMod

def mod_pow(a, e):
    r = 1
    while e > 0:
        if e & 1:
            r = mod_mul(r, a)
        a = mod_mul(a, a)
        e >>= 1
    return r

inv2 = mod_pow(2, kMod - 2)
inv6 = mod_pow(6, kMod - 2)

def sum1_prefix(n):
    a = n % kMod
    b = (n + 1) % kMod
    return mod_mul(mod_mul(a, b), inv2)

def sum2_prefix(n):
    a = n % kMod
    b = (n + 1) % kMod
    c = (2 * (n % kMod) + 1) % kMod
    return mod_mul(mod_mul(mod_mul(a, b), c), inv6)

def sum3_prefix(n):
    s1 = sum1_prefix(n)
    return mod_mul(s1, s1)

def range_sum1(l, r):
    if l > r: return 0
    return mod_sub(sum1_prefix(r), sum1_prefix(l - 1))

def range_sum2(l, r):
    if l > r: return 0
    return mod_sub(sum2_prefix(r), sum2_prefix(l - 1))

def range_sum3(l, r):
    if l > r: return 0
    return mod_sub(sum3_prefix(r), sum3_prefix(l - 1))

def max_k_u2(n):
    lo, hi = 0, 2000000000
    while lo < hi:
        mid = lo + (hi - lo + 1) // 2
        val = mid * mid - mid - 1
        if val <= n:
            lo = mid
        else:
            hi = mid - 1
    return lo

def max_k_u3(n):
    lo, hi = 0, 2000000
    while lo < hi:
        mid = lo + (hi - lo + 1) // 2
        val = mid * mid * mid - mid * mid - 2 * mid + 1
        if val <= n:
            lo = mid
        else:
            hi = mid - 1
    return lo

def max_k_seed4(n):
    lo, hi = 0, 2000000
    while lo < hi:
        mid = lo + (hi - lo + 1) // 2
        val = mid * (mid - 1) * (mid * mid - mid - 1) - 1
        if val <= n:
            lo = mid
        else:
            hi = mid - 1
    return lo

def capped_product_except(v, skip, cap):
    prod = 1
    for i in range(len(v)):
        if i == skip: continue
        if prod > cap // v[i]: return False, 0
        prod *= v[i]
    return True, prod

def exact_M_k(k, n):
    q = deque([tuple()])
    seen = {tuple()}
    numbers = {1}
    
    while q:
        cur = q.popleft()
        for x in cur:
            if x <= n: numbers.add(x)
            
        ones = k - len(cur)
        if ones > 0:
            cap = (n + 1) // k
            ok, prod = capped_product_except(cur, -1, cap)
            if ok:
                y = k * prod - 1
                if y <= n:
                    nxt = list(cur)
                    nxt.append(y)
                    nxt.sort()
                    nxt_t = tuple(nxt)
                    if nxt_t not in seen:
                        seen.add(nxt_t)
                        q.append(nxt_t)
                        
        for i in range(len(cur)):
            if i > 0 and cur[i] == cur[i - 1]: continue
            x = cur[i]
            cap = (n + x) // k
            ok, prod_others = capped_product_except(cur, i, cap)
            if not ok: continue
            
            y = k * prod_others - x
            if y == 0: continue
            
            nxt = list(cur[:i] + cur[i + 1:])
            if y > 1:
                import bisect
                bisect.insort(nxt, y)
            if nxt and nxt[-1] > n: continue
            
            nxt_t = tuple(nxt)
            if nxt_t not in seen:
                seen.add(nxt_t)
                q.append(nxt_t)
                
    out = 0
    for x in numbers:
        out = (out + (x % kMod)) % kMod
    return out

def chain_sum_range(l, r, n, k2max, k3max):
    if l > r: return 0
    ans = range_sum1(l, r)
    
    if l <= k2max:
        rr = min(r, k2max)
        cnt = (rr - l + 1) % kMod
        s2 = range_sum2(l, rr)
        s1 = range_sum1(l, rr)
        add = mod_sub(mod_sub(s2, s1), cnt)
        ans = mod_add(ans, add)
        
    if l <= k3max:
        rr = min(r, k3max)
        cnt = (rr - l + 1) % kMod
        s3 = range_sum3(l, rr)
        s2 = range_sum2(l, rr)
        s1 = range_sum1(l, rr)
        add = mod_sub(mod_sub(s3, s2), mod_mul(2, s1))
        add = mod_add(add, cnt)
        ans = mod_add(ans, add)
        
    return ans

def solve():
    K = 10**18
    N = 10**18
    
    k2max = max_k_u2(N)
    k3max = max_k_u3(N)
    k0 = max_k_seed4(N)
    
    exact_to = min(K, k0)
    ans = 0
    for k in range(3, exact_to + 1):
        ans = mod_add(ans, exact_M_k(k, N))
        
    if K > exact_to:
        l = max(3, exact_to + 1)
        ans = mod_add(ans, chain_sum_range(l, K, N, k2max, k3max))
        
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.Queue;

public class Euler844 {
    static final long kMod = 1405695061L;

    static long modAdd(long a, long b) {
        a += b;
        if (a >= kMod)
            a -= kMod;
        return a;
    }

    static long modSub(long a, long b) {
        return (a >= b) ? (a - b) : (a + kMod - b);
    }

    static long modMul(long a, long b) {
        return (a * b) % kMod;
    }

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

    static final long inv2 = modPow(2, kMod - 2);
    static final long inv6 = modPow(6, kMod - 2);

    static long sum1Prefix(long n) {
        long a = n % kMod;
        long b = (n + 1) % kMod;
        return modMul(modMul(a, b), inv2);
    }

    static long sum2Prefix(long n) {
        long a = n % kMod;
        long b = (n + 1) % kMod;
        long c = (2 * (n % kMod) + 1) % kMod;
        return modMul(modMul(modMul(a, b), c), inv6);
    }

    static long sum3Prefix(long n) {
        long s1 = sum1Prefix(n);
        return modMul(s1, s1);
    }

    static long rangeSum1(long l, long r) {
        if (l > r)
            return 0;
        return modSub(sum1Prefix(r), sum1Prefix(l - 1));
    }

    static long rangeSum2(long l, long r) {
        if (l > r)
            return 0;
        return modSub(sum2Prefix(r), sum2Prefix(l - 1));
    }

    static long rangeSum3(long l, long r) {
        if (l > r)
            return 0;
        return modSub(sum3Prefix(r), sum3Prefix(l - 1));
    }

    static long maxKU2(long n) {
        long lo = 0, hi = 2000000000L;
        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            long val = mid * mid - mid - 1;
            if (val <= n)
                lo = mid;
            else
                hi = mid - 1;
        }
        return lo;
    }

    static java.math.BigInteger toBI(long v) {
        return java.math.BigInteger.valueOf(v);
    }

    static long maxKU3(long n) {
        long lo = 0, hi = 2000000L;
        java.math.BigInteger nBi = toBI(n);
        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            java.math.BigInteger mBi = toBI(mid);
            java.math.BigInteger val = mBi.pow(3).subtract(mBi.pow(2)).subtract(toBI(2).multiply(mBi))
                    .add(java.math.BigInteger.ONE);
            if (val.compareTo(nBi) <= 0)
                lo = mid;
            else
                hi = mid - 1;
        }
        return lo;
    }

    static long maxKSeed4(long n) {
        long lo = 0, hi = 2000000L;
        java.math.BigInteger nBi = toBI(n);
        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            java.math.BigInteger mBi = toBI(mid);
            java.math.BigInteger val = mBi.multiply(mBi.subtract(java.math.BigInteger.ONE))
                    .multiply(mBi.pow(2).subtract(mBi).subtract(java.math.BigInteger.ONE))
                    .subtract(java.math.BigInteger.ONE);
            if (val.compareTo(nBi) <= 0)
                lo = mid;
            else
                hi = mid - 1;
        }
        return lo;
    }

    static boolean cappedProductExcept(ArrayList<Long> v, int skip, long cap, long[] out) {
        java.math.BigInteger prod = java.math.BigInteger.ONE;
        java.math.BigInteger capBi = toBI(cap);
        for (int i = 0; i < v.size(); ++i) {
            if (i == skip)
                continue;
            if (prod.compareTo(capBi.divide(toBI(v.get(i)))) > 0)
                return false;
            prod = prod.multiply(toBI(v.get(i)));
        }
        out[0] = prod.longValue();
        return true;
    }

    static long exactMK(long k, long n) {
        Queue<ArrayList<Long>> q = new LinkedList<>();
        HashSet<ArrayList<Long>> seen = new HashSet<>();
        HashSet<Long> numbers = new HashSet<>();

        ArrayList<Long> start = new ArrayList<>();
        q.add(start);
        seen.add(start);
        numbers.add(1L);

        while (!q.isEmpty()) {
            ArrayList<Long> cur = q.poll();

            for (long x : cur) {
                if (x <= n)
                    numbers.add(x);
            }

            long ones = k - cur.size();

            if (ones > 0) {
                long cap = (n + 1) / k;
                long[] prodArr = new long[1];
                if (cappedProductExcept(cur, -1, cap, prodArr)) {
                    long prod = prodArr[0];
                    long y = k * prod - 1;
                    if (y <= n) {
                        ArrayList<Long> nxt = new ArrayList<>(cur);
                        int pos = Collections.binarySearch(nxt, y);
                        if (pos < 0)
                            pos = -(pos + 1);
                        nxt.add(pos, y);
                        if (seen.add(nxt))
                            q.add(nxt);
                    }
                }
            }

            for (int i = 0; i < cur.size(); ++i) {
                if (i > 0 && cur.get(i).equals(cur.get(i - 1)))
                    continue;
                long x = cur.get(i);
                long cap = (n + x) / k;
                long[] prodArr = new long[1];
                if (!cappedProductExcept(cur, i, cap, prodArr))
                    continue;
                long prodOthers = prodArr[0];

                long y = k * prodOthers - x;
                if (y == 0)
                    continue;

                ArrayList<Long> nxt = new ArrayList<>(cur);
                nxt.remove(i);
                if (y > 1) {
                    int pos = Collections.binarySearch(nxt, y);
                    if (pos < 0)
                        pos = -(pos + 1);
                    nxt.add(pos, y);
                }
                if (!nxt.isEmpty() && nxt.get(nxt.size() - 1) > n)
                    continue;
                if (seen.add(nxt))
                    q.add(nxt);
            }
        }

        long out = 0;
        for (long x : numbers)
            out = modAdd(out, x % kMod);
        return out;
    }

    static long chainSumRange(long l, long r, long n, long k2max, long k3max) {
        if (l > r)
            return 0;

        long ans = rangeSum1(l, r);

        if (l <= k2max) {
            long rr = Math.min(r, k2max);
            long cnt = (rr - l + 1) % kMod;
            long s2 = rangeSum2(l, rr);
            long s1 = rangeSum1(l, rr);
            long add = modSub(modSub(s2, s1), cnt);
            ans = modAdd(ans, add);
        }

        if (l <= k3max) {
            long rr = Math.min(r, k3max);
            long cnt = (rr - l + 1) % kMod;
            long s3 = rangeSum3(l, rr);
            long s2 = rangeSum2(l, rr);
            long s1 = rangeSum1(l, rr);
            long add = modSub(modSub(s3, s2), modMul(2, s1));
            add = modAdd(add, cnt);
            ans = modAdd(ans, add);
        }

        return ans;
    }

    static long solveBig(long K, long N) {
        if (K < 3)
            return 0;

        long k2max = maxKU2(N);
        long k3max = maxKU3(N);
        long k0 = maxKSeed4(N);

        long exactTo = Math.min(K, k0);
        long ans = 0;
        for (long k = 3; k <= exactTo; ++k) {
            ans = modAdd(ans, exactMK(k, N));
        }

        if (K > exactTo) {
            long l = Math.max(3L, exactTo + 1);
            ans = modAdd(ans, chainSumRange(l, K, N, k2max, k3max));
        }

        return ans;
    }

    public static String solve() {
        return Long.toString(solveBig(1000000000000000000L, 1000000000000000000L));
    }

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