Problem 969: Kangaroo Hopping

Project Euler Problem 969 Solution

EulerSolve provides an optimized solution for Project Euler Problem 969, Kangaroo Hopping, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem defines an integer $S(n)$ by summing over all $t=1,2,\dots,n$. The term indexed by $t$ contributes only when the divisibility condition $(n-t)! \mid t^{\,n-t}$ holds, and then its value is $$(-1)^{n-t}\frac{t^{\,n-t}}{(n-t)!}.$$ The goal is not to compute one isolated $S(n)$, but the enormous prefix sum $$A(N)=\sum_{n=1}^{N} S(n)$$ for $N=10^{18}$, modulo $10^9+7$. A direct evaluation is hopeless, because even checking the divisibility condition term by term would already be far too slow. The successful approach turns that divisibility test into a simple description of which $t$-values are allowed, then reduces everything to a short sum of power sums. Mathematical Approach The implementations revolve around one re-indexing step, one valuation argument, and one fast method for evaluating $\sum_{k=1}^{x} k^m$ at huge $x$. Fix the gap $m=n-t$ Write $$m=n-t,\qquad t=n-m.$$ Then $m$ ranges from $0$ to $n-1$, and the defining sum becomes $$S(n)=\sum_{m=0}^{n-1}\mathbf{1}_{\,m!\mid (n-m)^m}\,(-1)^m\frac{(n-m)^m}{m!}.$$ For the prefix sum, it is better to regard $m$ as the outer variable and $t$ as the free value....

Detailed mathematical approach

Problem Summary

The problem defines an integer $S(n)$ by summing over all $t=1,2,\dots,n$. The term indexed by $t$ contributes only when the divisibility condition $(n-t)! \mid t^{\,n-t}$ holds, and then its value is

$$(-1)^{n-t}\frac{t^{\,n-t}}{(n-t)!}.$$

The goal is not to compute one isolated $S(n)$, but the enormous prefix sum

$$A(N)=\sum_{n=1}^{N} S(n)$$

for $N=10^{18}$, modulo $10^9+7$. A direct evaluation is hopeless, because even checking the divisibility condition term by term would already be far too slow. The successful approach turns that divisibility test into a simple description of which $t$-values are allowed, then reduces everything to a short sum of power sums.

Mathematical Approach

The implementations revolve around one re-indexing step, one valuation argument, and one fast method for evaluating $\sum_{k=1}^{x} k^m$ at huge $x$.

Fix the gap $m=n-t$

Write

$$m=n-t,\qquad t=n-m.$$

Then $m$ ranges from $0$ to $n-1$, and the defining sum becomes

$$S(n)=\sum_{m=0}^{n-1}\mathbf{1}_{\,m!\mid (n-m)^m}\,(-1)^m\frac{(n-m)^m}{m!}.$$

For the prefix sum, it is better to regard $m$ as the outer variable and $t$ as the free value. Since $n=t+m$, we obtain

$$A(N)=\sum_{m\ge 0}\frac{(-1)^m}{m!}\sum_{\substack{t\ge 1\\ t+m\le N\\ m!\mid t^m}} t^m.$$

Now the whole problem is: for each fixed $m$, characterize the integers $t$ for which $m!$ divides $t^m$.

The divisibility criterion is exactly a primorial criterion

Define

$$P_m=\prod_{\substack{p\le m\\ p\text{ prime}}} p,$$

with the empty product interpreted as $1$, so $P_0=P_1=1$. This is the primorial up to $m$.

The crucial fact used by all three implementations is

$$m!\mid t^m \quad\Longleftrightarrow\quad P_m\mid t.$$

Why is this true? For every prime $p\le m$, the factorial $m!$ contains at least one factor $p$. So if $p\nmid t$, then $v_p(t^m)=0\lt v_p(m!)$, and divisibility fails immediately. Therefore every prime $p\le m$ must divide $t$, which means $P_m\mid t$.

Conversely, suppose $P_m\mid t$. Then every prime $p\le m$ divides $t$, so $v_p(t)\ge 1$ and hence $v_p(t^m)\ge m$. On the other hand, Legendre's formula gives

$$v_p(m!)=\sum_{j\ge 1}\left\lfloor \frac{m}{p^j}\right\rfloor \le m,$$

so $v_p(t^m)\ge v_p(m!)$ for every prime $p\le m$. No larger prime matters, because it does not divide $m!$. Thus $m!\mid t^m$.

This is the decisive simplification: the complicated-looking factorial divisibility test collapses to the simple statement that $t$ is a multiple of one explicit primorial.

Turn admissible $t$ into a power sum

Once we know that the admissible values are exactly $t=kP_m$, the inner sum becomes static:

$$A(N)=\sum_{m\ge 0}\frac{(-1)^m}{m!}\sum_{kP_m+m\le N}(kP_m)^m.$$

Pull $P_m^m$ out of the inner sum:

$$A(N)=\sum_{m\ge 0}(-1)^m\frac{P_m^m}{m!}\sum_{k=1}^{\left\lfloor\frac{N-m}{P_m}\right\rfloor} k^m.$$

So for each $m$ we only need the power sum

$$F_m(x)=\sum_{k=1}^{x} k^m,\qquad x=\left\lfloor\frac{N-m}{P_m}\right\rfloor.$$

The outer loop is short because primorials grow very fast. As soon as $P_m>N-m$, the floor becomes $0$, and no later $m$ can contribute.

Worked example: evaluating $S(10)$

This small case shows how the primorial filter really selects the surviving terms. Here

$$S(10)=\sum_{t=1}^{10}\mathbf{1}_{(10-t)!\mid t^{\,10-t}}(-1)^{10-t}\frac{t^{\,10-t}}{(10-t)!}.$$

Write $m=10-t$.

For $m=4$, we have $P_4=2\cdot 3=6$. Then $t=6$ is a multiple of $6$, so the term survives and contributes

$$(-1)^4\frac{6^4}{4!}=\frac{1296}{24}=54.$$

For $m=3$, we have $P_3=6$, but now $t=7$ is not a multiple of $6$, so the term vanishes. For $m=2$, $P_2=2$ and $t=8$ is allowed, giving

$$(-1)^2\frac{8^2}{2!}=32.$$

The terms with $m=1$ and $m=0$ always survive, contributing $-9$ and $1$. All larger $m$ fail the divisibility test. Therefore

$$S(10)=54+32-9+1=78.$$

This example captures the full logic of the large computation: each $m$ contributes through multiples of $P_m$, and nothing else matters.

Evaluate $F_m(x)$ at astronomically large $x$

For fixed $m$, the power sum $F_m(x)$ is a polynomial in $x$ of degree $m+1$ by Faulhaber's theorem. Therefore it is completely determined by the $m+2$ values

$$F_m(0),F_m(1),\dots,F_m(m+1).$$

If we set $d=m+2$ and $y_i=F_m(i)$, then for any $x$

$$F_m(x)=\sum_{i=0}^{d-1} y_i \prod_{\substack{0\le j\le d-1\\ j\ne i}}\frac{x-j}{i-j}.$$

Because the interpolation nodes are consecutive integers, the denominator simplifies to

$$\prod_{\substack{0\le j\le d-1\\ j\ne i}}(i-j)=(-1)^{d-1-i} i!(d-1-i)!.$$

The implementations therefore precompute factorials and inverse factorials modulo $10^9+7$, then use prefix and suffix products of $(x-j)$ to evaluate the numerator of every Lagrange basis term in linear time.

This is especially efficient here because the contributing $m$-values are tiny compared with the modulus. For the target $N=10^{18}$, only $m=0,1,\dots,52$ contribute, so the interpolation tables never exceed $d=54$ entries.

How the Code Works

Exact checkpoint phase

Each implementation contains a small exact routine that evaluates $S(n)$ directly with arbitrary-precision integers. That exact routine is used only for validation: it confirms $S(1)=1$, $S(3)=-1$, the prefix $\sum_{n=1}^{10}S(n)=43$, and then compares several larger prefixes against the fast method. This guards the derivation before the large-$N$ computation is trusted.

Fast modular phase

The optimized solver walks upward through $m$, maintaining three quantities: the actual primorial $P_m$ for the stopping test, its residue modulo $10^9+7$, and $m!$ modulo $10^9+7$. For each $m$, it computes

$$(-1)^m P_m^m (m!)^{-1}\pmod{10^9+7}$$

and multiplies it by the modular value of $F_m\!\left(\left\lfloor\frac{N-m}{P_m}\right\rfloor\right)$.

The power sum is not expanded symbolically. Instead, the code first tabulates the prefix values $F_m(0),F_m(1),\dots,F_m(m+1)$, then performs one Lagrange evaluation modulo the prime modulus. The C++, Python, and Java implementations all follow this same mathematical pipeline; they only differ in language-specific integer types.

Complexity Analysis

The outer summation is extremely short. For $N=10^{18}$, the primorial up to $53$ already exceeds $N-53$, so only $53$ values of $m$ contribute: $m=0$ through $m=52$.

For each fixed $m$, building the sample table $F_m(0),\dots,F_m(m+1)$, the factorial tables, and the prefix/suffix products all costs $O(m)$ modular operations. Summing over all contributing $m$, the total work is therefore $O(m_{\max}^2)$ with $m_{\max}=52$, and the memory usage is $O(m_{\max})$. In concrete terms, the problem becomes small because the hard part is not the size of $N$, but the very slow growth of the set of relevant exponents $m$.

Footnotes and References

Problem page: https://projecteuler.net/problem=969
Primorial: Wikipedia - Primorial
p-adic valuation and Legendre's formula: Wikipedia - Legendre's formula
Power sums and Faulhaber's theorem: Wikipedia - Faulhaber's formula
Lagrange interpolation: Wikipedia - Lagrange polynomial

Problem 969 source code

C++

#include <cstdint>
#include <iostream>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using boost::multiprecision::cpp_int;

constexpr i64 MOD = 1'000'000'007LL;

i64 mod_pow(i64 base, i64 exp) {
    i64 result = 1;
    base %= MOD;
    while (exp > 0) {
        if (exp & 1LL) {
            result = static_cast<i64>((__int128)result * base % MOD);
        }
        base = static_cast<i64>((__int128)base * base % MOD);
        exp >>= 1LL;
    }
    return result;
}

bool is_prime_small(int n) {
    if (n < 2) {
        return false;
    }
    for (int d = 2; static_cast<i64>(d) * d <= n; ++d) {
        if (n % d == 0) {
            return false;
        }
    }
    return true;
}

i64 sum_powers_mod(u64 n, int k) {
    if (n == 0) {
        return 0;
    }

    const int d = k + 2;
    std::vector<i64> y(static_cast<std::size_t>(d), 0);
    for (int i = 1; i < d; ++i) {
        y[static_cast<std::size_t>(i)] =
            (y[static_cast<std::size_t>(i - 1)] + mod_pow(i, k)) % MOD;
    }

    const i64 x = static_cast<i64>(n % static_cast<u64>(MOD));
    if (x < d) {
        return y[static_cast<std::size_t>(x)];
    }

    std::vector<i64> fact(static_cast<std::size_t>(d), 1);
    std::vector<i64> inv_fact(static_cast<std::size_t>(d), 1);
    for (int i = 1; i < d; ++i) {
        fact[static_cast<std::size_t>(i)] =
            static_cast<i64>((__int128)fact[static_cast<std::size_t>(i - 1)] * i % MOD);
    }
    inv_fact[static_cast<std::size_t>(d - 1)] =
        mod_pow(fact[static_cast<std::size_t>(d - 1)], MOD - 2);
    for (int i = d - 1; i > 0; --i) {
        inv_fact[static_cast<std::size_t>(i - 1)] =
            static_cast<i64>((__int128)inv_fact[static_cast<std::size_t>(i)] * i % MOD);
    }

    std::vector<i64> prefix(static_cast<std::size_t>(d + 1), 1);
    std::vector<i64> suffix(static_cast<std::size_t>(d + 1), 1);

    for (int i = 0; i < d; ++i) {
        const i64 term = (x - i + MOD) % MOD;
        prefix[static_cast<std::size_t>(i + 1)] =
            static_cast<i64>((__int128)prefix[static_cast<std::size_t>(i)] * term % MOD);
    }
    for (int i = d - 1; i >= 0; --i) {
        const i64 term = (x - i + MOD) % MOD;
        suffix[static_cast<std::size_t>(i)] =
            static_cast<i64>((__int128)suffix[static_cast<std::size_t>(i + 1)] * term % MOD);
    }

    i64 answer = 0;
    for (int i = 0; i < d; ++i) {
        const i64 numerator =
            static_cast<i64>((__int128)prefix[static_cast<std::size_t>(i)] *
                                 suffix[static_cast<std::size_t>(i + 1)] %
                             MOD);

        i64 denominator = static_cast<i64>((__int128)inv_fact[static_cast<std::size_t>(i)] *
                                           inv_fact[static_cast<std::size_t>(d - 1 - i)] %
                                           MOD);
        if ((d - 1 - i) & 1) {
            denominator = (MOD - denominator) % MOD;
        }

        const i64 add = static_cast<i64>((__int128)y[static_cast<std::size_t>(i)] * numerator % MOD *
                                         denominator % MOD);
        answer += add;
        if (answer >= MOD) {
            answer -= MOD;
        }
    }

    return answer;
}

i64 solve_prefix(u64 n) {
    if (n == 0) {
        return 0;
    }

    i64 answer = 0;
    u128 primorial = 1;
    i64 primorial_mod = 1;
    i64 factorial_mod = 1;

    for (int m = 0;; ++m) {
        if (m >= 2 && is_prime_small(m)) {
            primorial *= static_cast<u128>(m);
            primorial_mod = static_cast<i64>((__int128)primorial_mod * m % MOD);
        }

        if (m > 0) {
            factorial_mod = static_cast<i64>((__int128)factorial_mod * m % MOD);
        }

        const u64 mu = static_cast<u64>(m);
        if (mu > n) {
            break;
        }
        const u64 remaining = n - mu;
        if (primorial > static_cast<u128>(remaining)) {
            break;
        }

        const u64 count = static_cast<u64>(static_cast<u128>(remaining) / primorial);

        i64 coefficient =
            static_cast<i64>((__int128)mod_pow(primorial_mod, m) *
                             mod_pow(factorial_mod, MOD - 2) %
                             MOD);
        if (m & 1) {
            coefficient = (MOD - coefficient) % MOD;
        }

        const i64 power_sum = sum_powers_mod(count, m);
        const i64 add = static_cast<i64>((__int128)coefficient * power_sum % MOD);

        answer += add;
        if (answer >= MOD) {
            answer -= MOD;
        }
    }

    return answer;
}

cpp_int direct_s(int n) {
    std::vector<cpp_int> factorial(static_cast<std::size_t>(n + 1), 1);
    for (int i = 1; i <= n; ++i) {
        factorial[static_cast<std::size_t>(i)] =
            factorial[static_cast<std::size_t>(i - 1)] * i;
    }

    cpp_int total = 0;
    for (int t = 1; t <= n; ++t) {
        const int m = n - t;
        cpp_int numerator = 1;
        for (int i = 0; i < m; ++i) {
            numerator *= t;
        }
        if (m & 1) {
            numerator = -numerator;
        }

        const cpp_int& denominator = factorial[static_cast<std::size_t>(m)];
        if (numerator % denominator == 0) {
            total += numerator / denominator;
        }
    }

    return total;
}

i64 cpp_int_mod(const cpp_int& value) {
    cpp_int residue = value % MOD;
    if (residue < 0) {
        residue += MOD;
    }
    return residue.convert_to<i64>();
}

bool run_checkpoints() {
    if (direct_s(1) != 1) {
        std::cerr << "Checkpoint failed: S(1) != 1" << '\n';
        return false;
    }

    if (direct_s(3) != -1) {
        std::cerr << "Checkpoint failed: S(3) != -1" << '\n';
        return false;
    }

    cpp_int prefix = 0;
    for (int n = 1; n <= 10; ++n) {
        prefix += direct_s(n);
    }
    if (prefix != 43) {
        std::cerr << "Checkpoint failed: sum_{n=1}^{10} S(n) != 43" << '\n';
        return false;
    }

    if (solve_prefix(10) != 43) {
        std::cerr << "Checkpoint failed: fast sum_{n=1}^{10} S(n) != 43" << '\n';
        return false;
    }

    for (int limit : {20, 30, 40}) {
        cpp_int direct_prefix = 0;
        for (int n = 1; n <= limit; ++n) {
            direct_prefix += direct_s(n);
        }

        const i64 fast = solve_prefix(static_cast<u64>(limit));
        const i64 slow = cpp_int_mod(direct_prefix);
        if (fast != slow) {
            std::cerr << "Checkpoint failed for N=" << limit << ": fast=" << fast
                      << ", slow=" << slow << '\n';
            return false;
        }
    }

    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 target = 1'000'000'000'000'000'000ULL;
    std::cout << solve_prefix(target) << '\n';
    return 0;
}

Python

import math

MOD = 1000000007

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

def is_prime_small(n):
    if n < 2:
        return False
    for d in range(2, int(n**0.5) + 1):
        if n % d == 0:
            return False
    return True

def sum_powers_mod(n, k):
    if n == 0:
        return 0

    d = k + 2
    y = [0] * d
    for i in range(1, d):
        y[i] = (y[i - 1] + mod_pow(i, k)) % MOD

    x = n % MOD
    if x < d:
        return y[x]

    fact = [1] * d
    inv_fact = [1] * d
    for i in range(1, d):
        fact[i] = (fact[i - 1] * i) % MOD
        
    inv_fact[d - 1] = pow(fact[d - 1], MOD - 2, MOD)
    for i in range(d - 1, 0, -1):
        inv_fact[i - 1] = (inv_fact[i] * i) % MOD

    prefix = [1] * (d + 1)
    suffix = [1] * (d + 1)

    for i in range(d):
        term = (x - i) % MOD
        prefix[i + 1] = (prefix[i] * term) % MOD
        
    for i in range(d - 1, -1, -1):
        term = (x - i) % MOD
        suffix[i] = (suffix[i + 1] * term) % MOD

    answer = 0
    for i in range(d):
        numerator = (prefix[i] * suffix[i + 1]) % MOD
        denominator = (inv_fact[i] * inv_fact[d - 1 - i]) % MOD
        if (d - 1 - i) % 2 == 1:
            denominator = (MOD - denominator) % MOD

        add = (y[i] * numerator % MOD) * denominator % MOD
        answer = (answer + add) % MOD

    return answer

def solve_prefix(n):
    if n == 0:
        return 0

    answer = 0
    primorial = 1
    primorial_mod = 1
    factorial_mod = 1

    m = 0
    while True:
        if m >= 2 and is_prime_small(m):
            primorial *= m
            primorial_mod = (primorial_mod * m) % MOD

        if m > 0:
            factorial_mod = (factorial_mod * m) % MOD

        if m > n:
            break
            
        remaining = n - m
        if primorial > remaining:
            break

        count = remaining // primorial

        coefficient = (mod_pow(primorial_mod, m) * pow(factorial_mod, MOD - 2, MOD)) % MOD
        if m % 2 == 1:
            coefficient = (MOD - coefficient) % MOD

        power_sum = sum_powers_mod(count, m)
        add = (coefficient * power_sum) % MOD

        answer = (answer + add) % MOD
        m += 1

    return answer

def direct_s(n):
    factorial = [1] * (n + 1)
    for i in range(1, n + 1):
        factorial[i] = factorial[i - 1] * i

    total = 0
    for t in range(1, n + 1):
        m = n - t
        numerator = 1
        for _ in range(m):
            numerator *= t
        if m % 2 == 1:
            numerator = -numerator

        denominator = factorial[m]
        if numerator % denominator == 0:
            total += numerator // denominator

    return total

def run_checkpoints():
    assert direct_s(1) == 1
    assert direct_s(3) == -1

    prefix = 0
    for n in range(1, 11):
        prefix += direct_s(n)
    assert prefix == 43
    assert solve_prefix(10) == 43

    for limit in [20, 30, 40]:
        direct_prefix = 0
        for n in range(1, limit + 1):
            direct_prefix += direct_s(n)

        fast = solve_prefix(limit)
        slow = direct_prefix % MOD
        assert fast == slow

def solve():
    target = 1000000000000000000
    return str(solve_prefix(target))

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

Java

import java.math.BigInteger;

public class Euler969 {

    static final long MOD = 1_000_000_007L;

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

    static boolean isPrimeSmall(int n) {
        if (n < 2)
            return false;
        for (int d = 2; (long) d * d <= n; ++d) {
            if (n % d == 0)
                return false;
        }
        return true;
    }

    static long sumPowersMod(long n, int k) {
        if (n == 0)
            return 0;

        int d = k + 2;
        long[] y = new long[d];
        for (int i = 1; i < d; ++i) {
            y[i] = (y[i - 1] + modPow(i, k)) % MOD;
        }

        long x = n % MOD;
        if (x < d)
            return y[(int) x];

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

        invFact[d - 1] = modPow(fact[d - 1], MOD - 2);
        for (int i = d - 1; i > 0; --i) {
            invFact[i - 1] = (invFact[i] * i) % MOD;
        }

        long[] prefix = new long[d + 1];
        long[] suffix = new long[d + 1];
        prefix[0] = 1;
        suffix[d] = 1;

        for (int i = 0; i < d; ++i) {
            long term = (x - i) % MOD;
            if (term < 0)
                term += MOD;
            prefix[i + 1] = (prefix[i] * term) % MOD;
        }
        for (int i = d - 1; i >= 0; --i) {
            long term = (x - i) % MOD;
            if (term < 0)
                term += MOD;
            suffix[i] = (suffix[i + 1] * term) % MOD;
        }

        long answer = 0;
        for (int i = 0; i < d; ++i) {
            long numerator = (prefix[i] * suffix[i + 1]) % MOD;
            long denominator = (invFact[i] * invFact[d - 1 - i]) % MOD;
            if (((d - 1 - i) & 1) != 0) {
                denominator = (MOD - denominator) % MOD;
            }

            long add = (y[i] * numerator) % MOD;
            add = (add * denominator) % MOD;
            answer = (answer + add) % MOD;
        }

        return answer;
    }

    static long solvePrefix(long n) {
        if (n == 0)
            return 0;

        long answer = 0;
        BigInteger primorial = BigInteger.ONE;
        long primorialMod = 1;
        long factorialMod = 1;

        for (int m = 0;; ++m) {
            if (m >= 2 && isPrimeSmall(m)) {
                primorial = primorial.multiply(BigInteger.valueOf(m));
                primorialMod = (primorialMod * m) % MOD;
            }

            if (m > 0) {
                factorialMod = (factorialMod * m) % MOD;
            }

            if (m > n)
                break;

            long remaining = n - m;
            if (primorial.compareTo(BigInteger.valueOf(remaining)) > 0) {
                break;
            }

            BigInteger[] divRem = BigInteger.valueOf(remaining).divideAndRemainder(primorial);
            long count = divRem[0].longValue();

            long coefficient = (modPow(primorialMod, m) * modPow(factorialMod, MOD - 2)) % MOD;
            if ((m & 1) != 0) {
                coefficient = (MOD - coefficient) % MOD;
            }

            long powerSum = sumPowersMod(count, m);
            long add = (coefficient * powerSum) % MOD;

            answer = (answer + add) % MOD;
        }

        return answer;
    }

    static BigInteger directS(int n) {
        BigInteger[] factorial = new BigInteger[n + 1];
        factorial[0] = BigInteger.ONE;
        for (int i = 1; i <= n; ++i) {
            factorial[i] = factorial[i - 1].multiply(BigInteger.valueOf(i));
        }

        BigInteger total = BigInteger.ZERO;
        for (int t = 1; t <= n; ++t) {
            int m = n - t;
            BigInteger numerator = BigInteger.ONE;
            for (int i = 0; i < m; ++i) {
                numerator = numerator.multiply(BigInteger.valueOf(t));
            }
            if ((m & 1) != 0) {
                numerator = numerator.negate();
            }

            BigInteger denominator = factorial[m];
            if (numerator.remainder(denominator).equals(BigInteger.ZERO)) {
                total = total.add(numerator.divide(denominator));
            }
        }

        return total;
    }

    static long bigIntMod(BigInteger value) {
        BigInteger residue = value.remainder(BigInteger.valueOf(MOD));
        if (residue.compareTo(BigInteger.ZERO) < 0) {
            residue = residue.add(BigInteger.valueOf(MOD));
        }
        return residue.longValue();
    }

    public static String solve() {
        long target = 1_000_000_000_000_000_000L;
        return Long.toString(solvePrefix(target));
    }

    public static void main(String[] args) {
        if (!directS(1).equals(BigInteger.ONE)) {
            System.out.println("Validation failed");
            return;
        }
        if (!directS(3).equals(BigInteger.valueOf(-1))) {
            System.out.println("Validation failed");
            return;
        }

        BigInteger prefix = BigInteger.ZERO;
        for (int n = 1; n <= 10; ++n) {
            prefix = prefix.add(directS(n));
        }
        if (!prefix.equals(BigInteger.valueOf(43))) {
            System.out.println("Validation failed");
            return;
        }
        if (solvePrefix(10) != 43) {
            System.out.println("Validation failed");
            return;
        }

        for (int limit : new int[] { 20, 30, 40 }) {
            BigInteger directPrefix = BigInteger.ZERO;
            for (int n = 1; n <= limit; ++n) {
                directPrefix = directPrefix.add(directS(n));
            }

            long fast = solvePrefix(limit);
            long slow = bigIntMod(directPrefix);
            if (fast != slow) {
                System.out.println("Validation failed for N=" + limit + ": fast=" + fast + ", slow=" + slow);
                return;
            }
        }

        System.out.println(solve());
    }
}