Problem 546: The Floor's Revenge

Project Euler Problem 546 Solution

EulerSolve provides an optimized solution for Project Euler Problem 546, The Floor's Revenge, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $k \ge 2$, define the sequence $$f_k(0)=1,\qquad f_k(n)=f_k(n-1)+f_k\left(\left\lfloor\frac{n}{k}\right\rfloor\right)\qquad (n\ge 1).$$ The required value is $$\sum_{k=2}^{10} f_k(N)\pmod{10^9+7},$$ with $N=10^{14}$ in the standard instance. A direct dynamic program up to $N$ is impossible, so the solution converts the recurrence into a digit-by-digit polynomial computation. Mathematical Approach Fix one base $k$, and write $f(n)=f_k(n)$. Let the base-$k$ expansion of $n$ be $$n=d_0+d_1k+\cdots+d_Lk^L,\qquad 0\le d_i\le k-1.$$ The aim is to evaluate $f(n)$ without iterating through every integer from $1$ to $n$. Step 1: Rewrite the Recurrence as a Cumulative Sum Since $f(0)=1$, the recurrence telescopes into $$f(n)=\sum_{m=0}^{n} f\left(\left\lfloor\frac{m}{k}\right\rfloor\right).$$ Now write $n=kx+s$ with $0\le s\le k-1$. Grouping the integers $m$ by the quotient $q=\lfloor m/k\rfloor$ gives $$f(kx+s)=k\sum_{q=0}^{x-1} f(q)+(s+1)f(x).$$ So one base-$k$ digit tells us how a prefix sum of previous values is sampled. Step 2: Expand Repeatedly Until Only Digits Remain Apply the same identity again to $f(x)$, then to the next quotient, and continue until the quotient becomes $0$....

Detailed mathematical approach

Problem Summary

For each integer $k \ge 2$, define the sequence

$$f_k(0)=1,\qquad f_k(n)=f_k(n-1)+f_k\left(\left\lfloor\frac{n}{k}\right\rfloor\right)\qquad (n\ge 1).$$

The required value is

$$\sum_{k=2}^{10} f_k(N)\pmod{10^9+7},$$

with $N=10^{14}$ in the standard instance. A direct dynamic program up to $N$ is impossible, so the solution converts the recurrence into a digit-by-digit polynomial computation.

Mathematical Approach

Fix one base $k$, and write $f(n)=f_k(n)$. Let the base-$k$ expansion of $n$ be

$$n=d_0+d_1k+\cdots+d_Lk^L,\qquad 0\le d_i\le k-1.$$

The aim is to evaluate $f(n)$ without iterating through every integer from $1$ to $n$.

Step 1: Rewrite the Recurrence as a Cumulative Sum

Since $f(0)=1$, the recurrence telescopes into

$$f(n)=\sum_{m=0}^{n} f\left(\left\lfloor\frac{m}{k}\right\rfloor\right).$$

Now write $n=kx+s$ with $0\le s\le k-1$. Grouping the integers $m$ by the quotient $q=\lfloor m/k\rfloor$ gives

$$f(kx+s)=k\sum_{q=0}^{x-1} f(q)+(s+1)f(x).$$

So one base-$k$ digit tells us how a prefix sum of previous values is sampled.

Step 2: Expand Repeatedly Until Only Digits Remain

Apply the same identity again to $f(x)$, then to the next quotient, and continue until the quotient becomes $0$. The result is an exact nested-sum formula:

$$f(n)=\sum_{x_L=0}^{d_L}\sum_{x_{L-1}=0}^{kx_L+d_{L-1}}\cdots\sum_{x_1=0}^{kx_2+d_1}\sum_{x_0=0}^{kx_1+d_0}1.$$

Each level introduces one variable whose upper bound is $k$ times the next variable plus the corresponding digit of $n$. The depth is only the number of base-$k$ digits, which is $O(\log_k n)$.

Step 3: Package One Layer as an Operator

For any function $P(x)$, define

$$T_s[P](x)=\sum_{t=0}^{kx+s} P(t),\qquad 0\le s\le k-1.$$

Then the nested sum can be written compactly as

$$f(n)=\left(T_{d_L}\circ T_{d_{L-1}}\circ\cdots\circ T_{d_0}\right)[1](0).$$

Thus the problem becomes: start from the constant polynomial $1$, and repeatedly apply the operator “take a prefix sum and evaluate it at $kx+s$” according to the digits of $n$.

Step 4: Why Polynomial States Are Enough

If $P(x)$ is a polynomial of degree $r$, then its prefix sum

$$S[P](x)=\sum_{t=0}^{x}P(t)$$

is a polynomial of degree $r+1$, and the substitution $x\mapsto kx+s$ preserves polynomial form. Therefore repeated applications of $T_s[P](x)=S[P](kx+s)$ never leave the polynomial world.

For a monomial $x^p$, the required prefix sum is

$$\sum_{t=0}^{x} t^p=\sum_{j=0}^{p} S(p,j)\,j!\binom{x+1}{j+1},$$

where $S(p,j)$ denotes a Stirling number of the second kind. This Faulhaber-type identity is exactly what lets the implementation precompute every power-sum polynomial once and reuse it.

Step 5: Handle the Upper Bound with Two Digit States

The implementation does not compose a single operator chain literally. Instead it scans the base-$k$ digits of $n$ from least significant to most significant and keeps two polynomial states:

one state for choices that still match the already processed part of $n$, and one state for choices that are already strictly smaller.

If the next digit is $d_j$, and $E_{j-1}$ and $B_{j-1}$ are the exact and below states from the previous step, then with $S[P](x)=\sum_{t=0}^{x}P(t)$ we have

$$B_j(x)=\sum_{s=0}^{k-1} S[B_{j-1}](kx+s),$$

$$E_j(x)=\sum_{s=0}^{d_j-1} S[B_{j-1}](kx+s)+S[E_{j-1}](kx+d_j).$$

The initial one-digit states are

$$B_0(x)=k,\qquad E_0(x)=d_0+1,$$

because a full last-digit block contributes $k$ choices, while the truncated block allowed by $d_0$ contributes $d_0+1$. After the final digit, the desired value is simply

$$f(n)=E_L(0).$$

Worked Example: $k=5,\ n=10$

Since $10=(20)_5$, the digits are $d_0=0$ and $d_1=2$. The nested-sum formula gives

$$f_5(10)=\sum_{x_1=0}^{2}\sum_{x_0=0}^{5x_1}1=1+6+11=18,$$

which matches the sample checkpoint used by the implementation.

The polynomial-state update gives the same answer. Initially, $B_0(x)=5$ and $E_0(x)=1$. Their prefix sums are $5x+5$ and $x+1$. Therefore

$$E_1(x)=\sum_{s=0}^{1}(25x+5s+5)+(5x+3)=55x+18,$$

$$f_5(10)=E_1(0)=18.$$

How the Code Works

The C++, Python, and Java implementations follow the same plan. They first determine a safe maximum polynomial degree from the binary length of $N$, because base $2$ yields the largest number of digits. Then they precompute binomial coefficients, Stirling numbers of the second kind, factorials, and the explicit polynomials for $\sum_{t=0}^{x} t^p$ for every degree that can appear.

For each $k=2,3,\dots,10$, the implementation converts $N$ to base $k$, initializes the exact and below states for the least significant digit, and processes the remaining digits one by one. Each digit step performs two generic algebraic operations:

$$P(x)\longmapsto \sum_{t=0}^{x}P(t),\qquad P(x)\longmapsto P(kx+s).$$

Because the states are stored by coefficients, both operations are done symbolically modulo $10^9+7$, never by iterating up to $N$. After the highest digit is processed, the exact-state polynomial is evaluated at $x=0$, giving $f_k(N)$. The final answer is the sum of these nine values modulo $10^9+7$.

Complexity Analysis

Let $D=\lfloor\log_2 N\rfloor+2$. This bounds the number of digits in every base from $2$ to $10$, so every polynomial degree is $O(D)$. The precomputation of binomial and power-sum data costs $O(D^3)$ time and $O(D^2)$ memory. For a fixed base $k$, processing all digits costs $O(kD^3)$ time in the worst case, because the degree grows linearly with the number of processed digits and each digit step performs a bounded number of $O(D^2)$ polynomial transforms. Here $k\le 10$ and $D$ is only about the number of binary digits of $10^{14}$, so the method is easily fast enough.

Footnotes and References

Problem page: https://projecteuler.net/problem=546
Positional notation: Wikipedia — Positional notation
Stirling numbers of the second kind: Wikipedia — Stirling numbers of the second kind
Faulhaber's formula: Wikipedia — Faulhaber's formula

Problem 546 source code

C++

#include <algorithm>
#include <cstdint>
#include <cstdlib>
#include <future>
#include <iostream>
#include <thread>
#include <vector>

namespace {

using i64 = long long;
using u64 = unsigned long long;

constexpr i64 kMod = 1000000007LL;
constexpr i64 kDefaultN = 100000000000000LL;

i64 mod_pow(i64 base, i64 exp) {
    i64 res = 1 % kMod;
    i64 cur = base % kMod;
    while (exp > 0) {
        if (exp & 1LL) res = (res * cur) % kMod;
        cur = (cur * cur) % kMod;
        exp >>= 1LL;
    }
    return res;
}

struct Precomputed {
    int max_deg = 0;
    std::vector<std::vector<i64>> binom;
    std::vector<std::vector<i64>> sum_pow;
};

std::vector<int> to_digits(i64 n, int base) {
    if (n == 0) return {0};
    std::vector<int> digits;
    while (n > 0) {
        digits.push_back(static_cast<int>(n % base));
        n /= base;
    }
    return digits;
}

Precomputed build_precomputed(int max_deg) {
    Precomputed pre;
    pre.max_deg = max_deg;
    const int size = max_deg + 2;

    std::vector<std::vector<i64>> binom(size, std::vector<i64>(size, 0));
    for (int n = 0; n < size; ++n) {
        binom[n][0] = 1;
        binom[n][n] = 1;
        for (int k = 1; k < n; ++k) {
            binom[n][k] = binom[n - 1][k - 1] + binom[n - 1][k];
            if (binom[n][k] >= kMod) binom[n][k] -= kMod;
        }
    }

    std::vector<std::vector<i64>> stirling(size, std::vector<i64>(size, 0));
    stirling[0][0] = 1;
    for (int n = 1; n < size; ++n) {
        for (int k = 1; k <= n; ++k) {
            stirling[n][k] = (stirling[n - 1][k - 1] + k * stirling[n - 1][k]) % kMod;
        }
    }

    std::vector<i64> fact(size + 1, 1);
    for (int i = 1; i < static_cast<int>(fact.size()); ++i) {
        fact[i] = (fact[i - 1] * i) % kMod;
    }

    std::vector<std::vector<i64>> poly_binom(size);
    poly_binom[0] = {1};
    for (int m = 1; m < size; ++m) {
        const auto& prev = poly_binom[m - 1];
        std::vector<i64> cur(prev.size() + 1, 0);
        for (int p = 0; p < static_cast<int>(prev.size()); ++p) {
            const i64 a = prev[p];
            const i64 dec = (a * (m - 1)) % kMod;
            cur[p] -= dec;
            if (cur[p] < 0) cur[p] += kMod;
            cur[p + 1] += a;
            if (cur[p + 1] >= kMod) cur[p + 1] -= kMod;
        }
        const i64 inv_m = mod_pow(m, kMod - 2);
        for (auto& v : cur) {
            v = (v * inv_m) % kMod;
        }
        poly_binom[m] = std::move(cur);
    }

    std::vector<std::vector<i64>> poly_binom_shift(size);
    for (int m = 0; m < size; ++m) {
        const auto& poly = poly_binom[m];
        std::vector<i64> shift(poly.size(), 0);
        for (int p = 0; p < static_cast<int>(poly.size()); ++p) {
            const i64 a = poly[p];
            if (a == 0) continue;
            for (int q = 0; q <= p; ++q) {
                i64 add = (a * binom[p][q]) % kMod;
                shift[q] += add;
                if (shift[q] >= kMod) shift[q] -= kMod;
            }
        }
        poly_binom_shift[m] = std::move(shift);
    }

    std::vector<std::vector<i64>> sum_pow(max_deg + 1);
    for (int p = 0; p <= max_deg; ++p) {
        std::vector<i64> poly(p + 2, 0);
        for (int j = 0; j <= p; ++j) {
            if (stirling[p][j] == 0) continue;
            const i64 coeff = (stirling[p][j] * fact[j]) % kMod;
            const auto& base = poly_binom_shift[j + 1];  // C(x+1, j+1)
            for (int idx = 0; idx < static_cast<int>(base.size()); ++idx) {
                i64 add = (coeff * base[idx]) % kMod;
                poly[idx] += add;
                if (poly[idx] >= kMod) poly[idx] -= kMod;
            }
        }
        sum_pow[p] = std::move(poly);
    }

    pre.binom = std::move(binom);
    pre.sum_pow = std::move(sum_pow);
    return pre;
}

std::vector<i64> prefix_sum_poly(const std::vector<i64>& poly, const Precomputed& pre) {
    std::vector<i64> res(poly.size() + 1, 0);
    for (int p = 0; p < static_cast<int>(poly.size()); ++p) {
        const i64 a = poly[p];
        if (a == 0) continue;
        const auto& base = pre.sum_pow[p];
        for (int i = 0; i < static_cast<int>(base.size()); ++i) {
            i64 add = (a * base[i]) % kMod;
            res[i] += add;
            if (res[i] >= kMod) res[i] -= kMod;
        }
    }
    return res;
}

std::vector<i64> substitute_poly(const std::vector<i64>& poly, int k, int s,
                                 const Precomputed& pre, const std::vector<i64>& pow_k) {
    const int d = static_cast<int>(poly.size()) - 1;
    std::vector<i64> pow_s(d + 1, 1);
    for (int i = 1; i <= d; ++i) {
        pow_s[i] = (pow_s[i - 1] * s) % kMod;
    }

    std::vector<i64> res(d + 1, 0);
    for (int p = 0; p <= d; ++p) {
        const i64 a = poly[p];
        if (a == 0) continue;
        for (int q = 0; q <= p; ++q) {
            i64 term = (a * pre.binom[p][q]) % kMod;
            term = (term * pow_k[q]) % kMod;
            term = (term * pow_s[p - q]) % kMod;
            res[q] += term;
            if (res[q] >= kMod) res[q] -= kMod;
        }
    }
    return res;
}

i64 compute_f(i64 n, int k, const Precomputed& pre) {
    if (n == 0) return 1;
    // Carry-based digit DP: polynomials encode counts as a function of the next carry.
    std::vector<int> digits = to_digits(n, k);
    const int L = static_cast<int>(digits.size()) - 1;

    std::vector<i64> pow_k(pre.max_deg + 2, 1);
    for (int i = 1; i < static_cast<int>(pow_k.size()); ++i) {
        pow_k[i] = (pow_k[i - 1] * k) % kMod;
    }

    std::vector<i64> free_poly(1, k % kMod);
    std::vector<i64> tight_poly(1, (digits[0] + 1) % kMod);
    if (L == 0) return tight_poly[0];

    for (int j = 1; j <= L; ++j) {
        std::vector<i64> sum_free = prefix_sum_poly(free_poly, pre);
        std::vector<i64> sum_tight = prefix_sum_poly(tight_poly, pre);

        std::vector<i64> new_free(sum_free.size(), 0);
        for (int s = 0; s < k; ++s) {
            std::vector<i64> sub = substitute_poly(sum_free, k, s, pre, pow_k);
            for (int i = 0; i < static_cast<int>(sub.size()); ++i) {
                new_free[i] += sub[i];
                if (new_free[i] >= kMod) new_free[i] -= kMod;
            }
        }

        std::vector<i64> new_tight(sum_free.size(), 0);
        const int dj = digits[j];
        for (int s = 0; s < dj; ++s) {
            std::vector<i64> sub = substitute_poly(sum_free, k, s, pre, pow_k);
            for (int i = 0; i < static_cast<int>(sub.size()); ++i) {
                new_tight[i] += sub[i];
                if (new_tight[i] >= kMod) new_tight[i] -= kMod;
            }
        }
        std::vector<i64> sub_tight = substitute_poly(sum_tight, k, dj, pre, pow_k);
        for (int i = 0; i < static_cast<int>(sub_tight.size()); ++i) {
            new_tight[i] += sub_tight[i];
            if (new_tight[i] >= kMod) new_tight[i] -= kMod;
        }

        free_poly = std::move(new_free);
        tight_poly = std::move(new_tight);
    }

    return tight_poly[0] % kMod;
}

u64 naive_f(int k, int n) {
    std::vector<u64> g(n + 1, 0);
    g[0] = 1;
    for (int i = 1; i <= n; ++i) {
        g[i] = g[i - 1] + g[i / k];
    }
    return g[n];
}

void run_validation(const Precomputed& pre) {
    struct Example { int k; int n; u64 expected; };
    const Example examples[] = {
        {5, 10, 18ULL},
        {7, 100, 1003ULL},
        {2, 1000, 264830889564ULL},
    };
    for (const auto& ex : examples) {
        const i64 got = compute_f(ex.n, ex.k, pre);
        const i64 want = static_cast<i64>(ex.expected % kMod);
        if (got != want) {
            std::cerr << "Validation failed for k=" << ex.k << " n=" << ex.n
                      << ": got " << got << " expected " << want << '\n';
            std::exit(1);
        }
    }

    for (int k = 2; k <= 10; ++k) {
        const int n = 200;
        const u64 exact = naive_f(k, n);
        const i64 got = compute_f(n, k, pre);
        if (got != static_cast<i64>(exact % kMod)) {
            std::cerr << "Validation failed for k=" << k << " n=" << n
                      << ": got " << got << " expected " << (exact % kMod) << '\n';
            std::exit(1);
        }
    }
}

}  // namespace

int main(int argc, char** argv) {
    i64 n = kDefaultN;
    if (argc > 1) {
        n = std::strtoll(argv[1], nullptr, 10);
    }

    const i64 max_n = std::max<i64>(n, 1000);
    const int max_deg = static_cast<int>(to_digits(max_n, 2).size()) + 2;
    const Precomputed pre = build_precomputed(max_deg);

    run_validation(pre);

    std::vector<std::future<i64>> futures;
    futures.reserve(9);
    for (int k = 2; k <= 10; ++k) {
        futures.emplace_back(std::async(std::launch::async, [&, k]() {
            return compute_f(n, k, pre);
        }));
    }

    i64 total = 0;
    for (auto& fut : futures) {
        i64 val = fut.get();
        total += val;
        if (total >= kMod) total -= kMod;
    }

    std::cout << total % kMod << '\n';
    return 0;
}

Python

MOD = 1000000007

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def to_digits(n, base):
    if n == 0: return [0]
    digits = []
    while n > 0:
        digits.append(n % base)
        n //= base
    return digits

class Precomputed:
    def __init__(self, max_deg):
        self.max_deg = max_deg
        size = max_deg + 2
        
        binom = [[0] * size for _ in range(size)]
        for n in range(size):
            binom[n][0] = 1
            binom[n][n] = 1
            for k in range(1, n):
                binom[n][k] = (binom[n - 1][k - 1] + binom[n - 1][k]) % MOD
        self.binom = binom
        
        stirling = [[0] * size for _ in range(size)]
        stirling[0][0] = 1
        for n in range(1, size):
            for k in range(1, n + 1):
                stirling[n][k] = (stirling[n - 1][k - 1] + k * stirling[n - 1][k]) % MOD
                
        fact = [1] * (size + 1)
        for i in range(1, size + 1):
            fact[i] = (fact[i - 1] * i) % MOD
            
        poly_binom = [[] for _ in range(size)]
        poly_binom[0] = [1]
        for m in range(1, size):
            prev = poly_binom[m - 1]
            cur = [0] * (len(prev) + 1)
            for p in range(len(prev)):
                a = prev[p]
                dec = (a * (m - 1)) % MOD
                cur[p] = (cur[p] - dec) % MOD
                cur[p + 1] = (cur[p + 1] + a) % MOD
                
            inv_m = mod_pow(m, MOD - 2)
            for i in range(len(cur)):
                cur[i] = (cur[i] * inv_m) % MOD
            poly_binom[m] = cur
            
        poly_binom_shift = [[] for _ in range(size)]
        for m in range(size):
            poly = poly_binom[m]
            shift = [0] * len(poly)
            for p in range(len(poly)):
                a = poly[p]
                if a == 0: continue
                for q in range(p + 1):
                    shift[q] = (shift[q] + a * binom[p][q]) % MOD
            poly_binom_shift[m] = shift
            
        sum_pow = [[] for _ in range(max_deg + 1)]
        for p in range(max_deg + 1):
            poly = [0] * (p + 2)
            for j in range(p + 1):
                if stirling[p][j] == 0: continue
                coeff = (stirling[p][j] * fact[j]) % MOD
                base = poly_binom_shift[j + 1]
                for idx in range(len(base)):
                    poly[idx] = (poly[idx] + coeff * base[idx]) % MOD
            sum_pow[p] = poly
            
        self.sum_pow = sum_pow

def prefix_sum_poly(poly, pre):
    res = [0] * (len(poly) + 1)
    for p in range(len(poly)):
        a = poly[p]
        if a == 0: continue
        base = pre.sum_pow[p]
        for i in range(len(base)):
            res[i] = (res[i] + a * base[i]) % MOD
    return res

def substitute_poly(poly, k, s, pre, pow_k):
    d = len(poly) - 1
    pow_s = [1] * (d + 1)
    pow_s[0] = 1
    for i in range(1, d + 1):
        pow_s[i] = (pow_s[i - 1] * s) % MOD
        
    res = [0] * (d + 1)
    for p in range(d + 1):
        a = poly[p]
        if a == 0: continue
        for q in range(p + 1):
            term = (a * pre.binom[p][q]) % MOD
            term = (term * pow_k[q]) % MOD
            term = (term * pow_s[p - q]) % MOD
            res[q] = (res[q] + term) % MOD
    return res

def compute_f(n, k, pre):
    if n == 0: return 1
    digits = to_digits(n, k)
    L = len(digits) - 1
    
    pow_k = [1] * (pre.max_deg + 2)
    for i in range(1, len(pow_k)):
        pow_k[i] = (pow_k[i - 1] * k) % MOD
        
    free_poly = [k % MOD]
    tight_poly = [(digits[0] + 1) % MOD]
    if L == 0: return tight_poly[0]
    
    for j in range(1, L + 1):
        sum_free = prefix_sum_poly(free_poly, pre)
        sum_tight = prefix_sum_poly(tight_poly, pre)
        
        new_free = [0] * len(sum_free)
        for s in range(k):
            sub = substitute_poly(sum_free, k, s, pre, pow_k)
            for i in range(len(sub)):
                new_free[i] = (new_free[i] + sub[i]) % MOD
                
        new_tight = [0] * len(sum_free)
        dj = digits[j]
        for s in range(dj):
            sub = substitute_poly(sum_free, k, s, pre, pow_k)
            for i in range(len(sub)):
                new_tight[i] = (new_tight[i] + sub[i]) % MOD
                
        sub_tight = substitute_poly(sum_tight, k, dj, pre, pow_k)
        for i in range(len(sub_tight)):
            new_tight[i] = (new_tight[i] + sub_tight[i]) % MOD
            
        free_poly = new_free
        tight_poly = new_tight
        
    return tight_poly[0] % MOD

def solve():
    n = 100000000000000
    max_n = max(n, 1000)
    max_deg = len(to_digits(max_n, 2)) + 2
    pre = Precomputed(max_deg)
    
    total = 0
    for k in range(2, 11):
        total = (total + compute_f(n, k, pre)) % MOD
            
    return str(total)

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

Java

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

public class Euler546 {
    static final long MOD = 1000000007L;

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

    static class Precomputed {
        int maxDeg;
        long[][] binom;
        long[][] sumPow;

        Precomputed(int maxDeg, long[][] binom, long[][] sumPow) {
            this.maxDeg = maxDeg;
            this.binom = binom;
            this.sumPow = sumPow;
        }
    }

    static List<Integer> toDigits(long n, int base) {
        if (n == 0) {
            List<Integer> list = new ArrayList<>();
            list.add(0);
            return list;
        }
        List<Integer> digits = new ArrayList<>();
        while (n > 0) {
            digits.add((int) (n % base));
            n /= base;
        }
        return digits;
    }

    static Precomputed buildPrecomputed(int maxDeg) {
        int size = maxDeg + 2;
        long[][] binom = new long[size][size];
        for (int n = 0; n < size; n++) {
            binom[n][0] = 1;
            binom[n][n] = 1;
            for (int k = 1; k < n; k++) {
                binom[n][k] = (binom[n - 1][k - 1] + binom[n - 1][k]) % MOD;
            }
        }

        long[][] stirling = new long[size][size];
        stirling[0][0] = 1;
        for (int n = 1; n < size; n++) {
            for (int k = 1; k <= n; k++) {
                stirling[n][k] = (stirling[n - 1][k - 1] + (long) k * stirling[n - 1][k]) % MOD;
            }
        }

        long[] fact = new long[size + 1];
        fact[0] = 1;
        for (int i = 1; i < fact.length; i++) {
            fact[i] = (fact[i - 1] * i) % MOD;
        }

        long[][] polyBinom = new long[size][];
        polyBinom[0] = new long[] { 1 };
        for (int m = 1; m < size; m++) {
            long[] prev = polyBinom[m - 1];
            long[] cur = new long[prev.length + 1];
            for (int p = 0; p < prev.length; p++) {
                long a = prev[p];
                long dec = (a * (m - 1)) % MOD;
                cur[p] = (cur[p] - dec + MOD) % MOD;
                cur[p + 1] = (cur[p + 1] + a) % MOD;
            }
            long invM = modPow(m, MOD - 2);
            for (int i = 0; i < cur.length; i++) {
                cur[i] = (cur[i] * invM) % MOD;
            }
            polyBinom[m] = cur;
        }

        long[][] polyBinomShift = new long[size][];
        for (int m = 0; m < size; m++) {
            long[] poly = polyBinom[m];
            long[] shift = new long[poly.length];
            for (int p = 0; p < poly.length; p++) {
                long a = poly[p];
                if (a == 0)
                    continue;
                for (int q = 0; q <= p; q++) {
                    long add = (a * binom[p][q]) % MOD;
                    shift[q] = (shift[q] + add) % MOD;
                }
            }
            polyBinomShift[m] = shift;
        }

        long[][] sumPow = new long[maxDeg + 1][];
        for (int p = 0; p <= maxDeg; p++) {
            long[] poly = new long[p + 2];
            for (int j = 0; j <= p; j++) {
                if (stirling[p][j] == 0)
                    continue;
                long coeff = (stirling[p][j] * fact[j]) % MOD;
                long[] baseArr = polyBinomShift[j + 1];
                for (int idx = 0; idx < baseArr.length; idx++) {
                    long add = (coeff * baseArr[idx]) % MOD;
                    poly[idx] = (poly[idx] + add) % MOD;
                }
            }
            sumPow[p] = poly;
        }

        return new Precomputed(maxDeg, binom, sumPow);
    }

    static long[] prefixSumPoly(long[] poly, Precomputed pre) {
        long[] res = new long[poly.length + 1];
        for (int p = 0; p < poly.length; p++) {
            long a = poly[p];
            if (a == 0)
                continue;
            long[] base = pre.sumPow[p];
            for (int i = 0; i < base.length; i++) {
                long add = (a * base[i]) % MOD;
                res[i] = (res[i] + add) % MOD;
            }
        }
        return res;
    }

    static long[] substitutePoly(long[] poly, int k, int s, Precomputed pre, long[] powK) {
        int d = poly.length - 1;
        long[] powS = new long[d + 1];
        powS[0] = 1;
        for (int i = 1; i <= d; i++) {
            powS[i] = (powS[i - 1] * s) % MOD;
        }

        long[] res = new long[d + 1];
        for (int p = 0; p <= d; p++) {
            long a = poly[p];
            if (a == 0)
                continue;
            for (int q = 0; q <= p; q++) {
                long term = (a * pre.binom[p][q]) % MOD;
                term = (term * powK[q]) % MOD;
                term = (term * powS[p - q]) % MOD;
                res[q] = (res[q] + term) % MOD;
            }
        }
        return res;
    }

    static long computeF(long n, int k, Precomputed pre) {
        if (n == 0)
            return 1;
        List<Integer> digits = toDigits(n, k);
        int L = digits.size() - 1;

        long[] powK = new long[pre.maxDeg + 2];
        powK[0] = 1;
        for (int i = 1; i < powK.length; i++) {
            powK[i] = (powK[i - 1] * k) % MOD;
        }

        long[] freePoly = { k % MOD };
        long[] tightPoly = { (digits.get(0) + 1) % MOD };
        if (L == 0)
            return tightPoly[0];

        for (int j = 1; j <= L; j++) {
            long[] sumFree = prefixSumPoly(freePoly, pre);
            long[] sumTight = prefixSumPoly(tightPoly, pre);

            long[] newFree = new long[sumFree.length];
            for (int s = 0; s < k; s++) {
                long[] sub = substitutePoly(sumFree, k, s, pre, powK);
                for (int i = 0; i < sub.length; i++) {
                    newFree[i] = (newFree[i] + sub[i]) % MOD;
                }
            }

            long[] newTight = new long[sumFree.length];
            int dj = digits.get(j);
            for (int s = 0; s < dj; s++) {
                long[] sub = substitutePoly(sumFree, k, s, pre, powK);
                for (int i = 0; i < sub.length; i++) {
                    newTight[i] = (newTight[i] + sub[i]) % MOD;
                }
            }
            long[] subTight = substitutePoly(sumTight, k, dj, pre, powK);
            for (int i = 0; i < subTight.length; i++) {
                newTight[i] = (newTight[i] + subTight[i]) % MOD;
            }

            freePoly = newFree;
            tightPoly = newTight;
        }

        return tightPoly[0] % MOD;
    }

    public static String solve() {
        long n = 100000000000000L;
        long maxN = Math.max(n, 1000L);
        int maxDeg = toDigits(maxN, 2).size() + 2;
        Precomputed pre = buildPrecomputed(maxDeg);

        long total = 0;
        for (int k = 2; k <= 10; k++) {
            long val = computeF(n, k, pre);
            total = (total + val) % MOD;
        }
        return Long.toString(total);
    }

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