Problem 767: Window into a Matrix II

Project Euler Problem 767 Solution

EulerSolve provides an optimized solution for Project Euler Problem 767, Window into a Matrix II, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The quantity to compute is $B(k,n)$ modulo $M=10^9+7$. The C++, Python, and Java implementations do not enumerate matrix states directly. Instead, they rewrite the target value as a coefficient extraction problem built from factorial-based series. The difficult part is a family of sums involving sixteenth powers of binomial coefficients, and the implementations evaluate all of those sums at once by fast polynomial convolution. For the target instance, $k$ is as large as $100000$, so a direct quadratic summation would be far too slow. The entire strategy is therefore organized around turning the problem into two convolutions of length about $k$, followed by a final modular reconstruction. Mathematical Approach Write $$q=\left\lfloor\frac{n}{k}\right\rfloor,\qquad M=10^9+7,\qquad \mu\equiv 2^q-2 \pmod{M}.$$ The implemented formula can be derived cleanly in the following steps. Step 1: Reduce the exponential parameter All dependence on $n$ enters through the single scalar $\mu$. Because $M$ is prime and $2\not\equiv 0\pmod{M}$, Fermat's little theorem gives $$2^q \equiv 2^{\,q\bmod (M-1)} \pmod{M}.$$ So even when $q$ is enormous, the implementation only needs the reduced exponent $q\bmod(M-1)$. After this reduction, the rest of the computation depends only on $k$ and on the already-computed value of $\mu$....

Detailed mathematical approach

Problem Summary

The quantity to compute is $B(k,n)$ modulo $M=10^9+7$. The C++, Python, and Java implementations do not enumerate matrix states directly. Instead, they rewrite the target value as a coefficient extraction problem built from factorial-based series. The difficult part is a family of sums involving sixteenth powers of binomial coefficients, and the implementations evaluate all of those sums at once by fast polynomial convolution.

For the target instance, $k$ is as large as $100000$, so a direct quadratic summation would be far too slow. The entire strategy is therefore organized around turning the problem into two convolutions of length about $k$, followed by a final modular reconstruction.

Mathematical Approach

Write

$$q=\left\lfloor\frac{n}{k}\right\rfloor,\qquad M=10^9+7,\qquad \mu\equiv 2^q-2 \pmod{M}.$$

The implemented formula can be derived cleanly in the following steps.

Step 1: Reduce the exponential parameter

All dependence on $n$ enters through the single scalar $\mu$. Because $M$ is prime and $2\not\equiv 0\pmod{M}$, Fermat's little theorem gives

$$2^q \equiv 2^{\,q\bmod (M-1)} \pmod{M}.$$

So even when $q$ is enormous, the implementation only needs the reduced exponent $q\bmod(M-1)$. After this reduction, the rest of the computation depends only on $k$ and on the already-computed value of $\mu$.

Step 2: Convert the inner binomial-power sum into a convolution

Define the sequence

$$c_i=(i!)^{-16}\pmod{M}\qquad (0\le i\le k).$$

Now form its self-convolution:

$$d_s=\sum_{i=0}^{s} c_i\,c_{s-i} \qquad (0\le s\le k).$$

Multiplying by $(s!)^{16}$ transforms this into

$$a_s=(s!)^{16}d_s=(s!)^{16}\sum_{i=0}^{s}\frac{1}{(i!)^{16}((s-i)!)^{16}}=\sum_{i=0}^{s}\left(\frac{s!}{i!(s-i)!}\right)^{16}.$$

Therefore

$$\boxed{a_s=\sum_{i=0}^{s}\binom{s}{i}^{16}}.$$

This identity is the key simplification: instead of evaluating each $a_s$ separately, one self-convolution computes every value $a_0,a_1,\dots,a_k$ in one shot.

Step 3: Rewrite $B(k,n)$ as one more coefficient extraction

After the quantities $a_s$ are known, define

$$g_s=\frac{a_s}{s!},\qquad e_j=\frac{\mu^j}{j!}.$$

Convolving these two sequences gives

$$f_t=\sum_{s=0}^{t} g_s\,e_{t-s}=\sum_{s=0}^{t}\frac{a_s}{s!}\frac{\mu^{\,t-s}}{(t-s)!}.$$

The final answer is the $t=k$ case multiplied by $k!$:

$$B(k,n)=k!\,f_k.$$

Expanding the factorials yields the explicit summation formula

$$\boxed{B(k,n)=\sum_{s=0}^{k}\binom{k}{s}a_s\,\mu^{\,k-s}\pmod{M}.}$$

Equivalently, if we package the coefficients into an exponential generating function, then

$$\frac{B(k,n)}{k!}=\left[x^k\right]\left(\sum_{s\ge 0}a_s\frac{x^s}{s!}\right)e^{\mu x}.$$

So the whole problem becomes a coefficient lookup after a second convolution.

Step 4: Use NTT and CRT for both convolutions

A naive evaluation of all $a_s$ from the formula $\sum_{i=0}^{s}\binom{s}{i}^{16}$ would cost

$$\sum_{s=0}^{k} O(s)=O(k^2),$$

which is too slow when $k=100000$. Both required convolutions therefore use the Number Theoretic Transform. The target modulus $M=10^9+7$ is not suitable for an NTT of the needed lengths, so the implementation performs each convolution under three NTT-friendly primes:

$$998244353,\qquad 1004535809,\qquad 469762049.$$

Each coefficient is first computed modulo these three primes and then reconstructed by the Chinese Remainder Theorem. The product of the three auxiliary moduli is far larger than the raw coefficient sizes that occur here, so the reconstruction is exact before the final reduction modulo $M$.

Step 5: Worked example

The small checkpoint $B(2,4)$ illustrates every layer of the formula. Here

$$q=\left\lfloor\frac{4}{2}\right\rfloor=2,\qquad \mu=2^2-2=2.$$

Next compute the first three values of $a_s$:

$$a_0=\binom{0}{0}^{16}=1,$$

$$a_1=\binom{1}{0}^{16}+\binom{1}{1}^{16}=1+1=2,$$

$$a_2=\binom{2}{0}^{16}+\binom{2}{1}^{16}+\binom{2}{2}^{16}=1+2^{16}+1=65538.$$

Substituting into the closed form gives

$$\begin{aligned} B(2,4) &=\binom{2}{0}a_0\mu^2+\binom{2}{1}a_1\mu+\binom{2}{2}a_2\\ &=1\cdot 1\cdot 4+2\cdot 2\cdot 2+1\cdot 65538\\ &=4+8+65538\\ &=65550. \end{aligned}$$

This matches the small-value check used by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they precompute factorials and inverse factorials modulo $M$ up to $k$. That makes it cheap to evaluate $(i!)^{-16}$, $(s!)^{16}$, and the factorial denominators that appear in the two generating series.

Next they compute $\mu=2^{\lfloor n/k\rfloor}-2\pmod{M}$ with fast modular exponentiation, reducing the exponent modulo $M-1$ before the power is taken. Then they build the sequence $c_i=(i!)^{-16}$, convolve it with itself, and rescale the result to recover every value $a_s=\sum_{i=0}^{s}\binom{s}{i}^{16}$.

After that, the implementation forms one series from $a_s/s!$ and a second series from $\mu^j/j!$, performs a second convolution, and multiplies the coefficient of degree $k$ by $k!$. The C++ version additionally parallelizes the three auxiliary-modulus convolutions and the CRT recombination, while the Python and Java versions execute the same mathematical steps sequentially.

Complexity Analysis

Precomputing factorials, inverse factorials, and the powers needed for the exponential series costs $O(k)$ time and $O(k)$ memory. Each NTT-based convolution costs $O(k\log k)$ time. There are two logical convolutions, each evaluated under a fixed set of three auxiliary primes, so the asymptotic running time remains $O(k\log k)$ and the memory usage remains $O(k)$. Parallel execution in the C++ implementation improves wall-clock time but does not change the asymptotic order.

Footnotes and References

Problem page: Project Euler 767
Binomial coefficients: Wikipedia - Binomial coefficient
Generating functions: Wikipedia - Generating function
Number Theoretic Transform: cp-algorithms - Fast Fourier transform and NTT
Chinese Remainder Theorem: Wikipedia - Chinese remainder theorem

Problem 767 source code

C++

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <thread>
#include <future>

using namespace std;

const long long TARGET_MOD = 1000000007;

long long power(long long base, long long exp, long long mod) {
    long long res = 1;
    base %= mod;
    while (exp > 0) {
        if (exp % 2 == 1) res = (__int128)res * base % mod;
        base = (__int128)base * base % mod;
        exp /= 2;
    }
    return res;
}

long long modInverse(long long n, long long mod) {
    return power(n, mod - 2, mod);
}

struct NTT_Mod {
    long long mod;
    long long root;
    long long root_pw;
};

const NTT_Mod P1 = {998244353, 3, 1 << 23};
const NTT_Mod P2 = {1004535809, 3, 1 << 21};
const NTT_Mod P3 = {469762049, 3, 1 << 26};

void ntt(vector<long long>& a, bool invert, const NTT_Mod& p) {
    int n = a.size();
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        long long wlen = power(p.root, (p.mod - 1) / len, p.mod);
        if (invert) wlen = modInverse(wlen, p.mod);
        for (int i = 0; i < n; i += len) {
            long long w = 1;
            for (int j = 0; j < len / 2; j++) {
                long long u = a[i + j];
                long long v = (__int128)a[i + j + len / 2] * w % p.mod;
                a[i + j] = (u + v < p.mod) ? u + v : u + v - p.mod;
                a[i + j + len / 2] = (u - v >= 0) ? u - v : u - v + p.mod;
                w = (__int128)w * wlen % p.mod;
            }
        }
    }

    if (invert) {
        long long n_inv = modInverse(n, p.mod);
        for (long long& x : a) x = (__int128)x * n_inv % p.mod;
    }
}

vector<long long> convolve(const vector<long long>& a, const vector<long long>& b, const NTT_Mod& p) {
    vector<long long> fa(a.begin(), a.end());
    vector<long long> fb(b.begin(), b.end());
    int n = 1;
    while (n < a.size() + b.size()) n <<= 1;
    fa.resize(n);
    fb.resize(n);

    ntt(fa, false, p);
    ntt(fb, false, p);
    for (int i = 0; i < n; i++) fa[i] = (__int128)fa[i] * fb[i] % p.mod;
    ntt(fa, true, p);

    vector<long long> result = fa;
    return result;
}

long long crt_combine(long long a1, long long a2, long long a3) {
    long long m1 = P1.mod, m2 = P2.mod, m3 = P3.mod;
    long long M = m1 * m2;
    long long inv_m1_m2 = modInverse(m1, m2);
    long long x = (a1 + (__int128)(a2 - a1 + m2) % m2 * inv_m1_m2 % m2 * m1);
    long long M_mod_m3 = M % m3;
    long long inv_M_m3 = modInverse(M_mod_m3, m3);
    long long k = (__int128)(a3 - (x % m3) + m3) % m3 * inv_M_m3 % m3;
    __int128 res = x + (__int128)k * M;
    return (long long)(res % TARGET_MOD);
}

vector<long long> threaded_convolution_crt(const vector<long long>& A, const vector<long long>& B) {
    auto f1 = async(launch::async, convolve, cref(A), cref(B), cref(P1));
    auto f2 = async(launch::async, convolve, cref(A), cref(B), cref(P2));
    auto f3 = async(launch::async, convolve, cref(A), cref(B), cref(P3));

    vector<long long> r1 = f1.get();
    vector<long long> r2 = f2.get();
    vector<long long> r3 = f3.get();

    int n = r1.size();
    vector<long long> res(n);

    int num_threads = thread::hardware_concurrency();
    vector<thread> threads;
    int chunk = n / num_threads + 1;
    
    for(int t=0; t<num_threads; ++t) {
        threads.emplace_back([&, t, chunk]() {
            int start = t * chunk;
            int end = min(start + chunk, n);
            for(int i=start; i<end; ++i) {
                res[i] = crt_combine(r1[i], r2[i], r3[i]);
            }
        });
    }
    for(auto& th : threads) th.join();
    
    return res;
}

class Solver {
    vector<long long> fact, invFact;
    int max_k;

public:
    Solver(int k) : max_k(k) {
        fact.resize(k + 1);
        invFact.resize(k + 1);
        fact[0] = 1;
        for (int i = 1; i <= k; i++) fact[i] = fact[i - 1] * i % TARGET_MOD;
        invFact[k] = modInverse(fact[k], TARGET_MOD);
        for (int i = k - 1; i >= 0; i--) invFact[i] = invFact[i + 1] * (i + 1) % TARGET_MOD;
    }

    long long nCr_pow16(int n, int r) {
        if (r < 0 || r > n) return 0;
        long long num = fact[n];
        long long den = (invFact[r] * invFact[n - r]) % TARGET_MOD;
        long long comb = (num * den) % TARGET_MOD;
        return power(comb, 16, TARGET_MOD); 
    }

    long long solve(int k, long long n) {
        long long phi = TARGET_MOD - 1;
        long long m_mod_phi = power(10, 11, phi); 
        long long m = n / k;
        long long m_exp = m % phi;
        
        long long lambda = (power(2, m_exp, TARGET_MOD) - 2 + TARGET_MOD) % TARGET_MOD;

        vector<long long> H(k + 1);
        for(int i=0; i<=k; ++i) {
            H[i] = power(invFact[i], 16, TARGET_MOD);
        }

        vector<long long> H2 = threaded_convolution_crt(H, H);

        vector<long long> U(k + 1);
        for(int s=0; s<=k; ++s) {
            long long A_s = (power(fact[s], 16, TARGET_MOD) * H2[s]) % TARGET_MOD;
            U[s] = (A_s * invFact[s]) % TARGET_MOD;
        }

        vector<long long> V(k + 1);
        long long lam_pow = 1;
        for(int j=0; j<=k; ++j) {
            V[j] = (lam_pow * invFact[j]) % TARGET_MOD;
            lam_pow = (lam_pow * lambda) % TARGET_MOD;
        }

        vector<long long> W = threaded_convolution_crt(U, V);

        long long ans = (fact[k] * W[k]) % TARGET_MOD;
        return ans;
    }
};

int main() {
    cout << "Running validation checks..." << endl;
    
    {
        Solver s(2);
        long long res = s.solve(2, 4);
        cout << "B(2, 4) = " << res << (res == 65550 ? " [PASS]" : " [FAIL]") << endl;
    }

    {
        Solver s(3);
        long long res = s.solve(3, 9);
        cout << "B(3, 9) = " << res << (res == 87273560 ? " [PASS]" : " [FAIL]") << endl;
    }

    cout << "-----------------------------------" << endl;
    
    int k_target = 100000;
    long long n_target = 10000000000000000LL;

    cout << "Solving for B(" << k_target << ", 10^16)..." << endl;
    
    Solver finalSolver(k_target); 
    
    auto start = chrono::high_resolution_clock::now();
    
    long long answer = finalSolver.solve(k_target, n_target);
    
    auto end = chrono::high_resolution_clock::now();
    chrono::duration<double> diff = end - start;

    cout << "Computation finished in " << diff.count() << " s" << endl;
    cout << "Answer: " << answer << endl;

    return 0;
}

Python

def solve():
    MOD = 1000000007
    k_target = 100000; n_target = 10000000000000000

    def mod_pow(base, exp, mod=MOD):
        r = 1; base %= mod
        while exp > 0:
            if exp & 1: r = r * base % mod
            base = base * base % mod; exp >>= 1
        return r

    def mod_inv(n, mod=MOD): return mod_pow(n, mod-2, mod)

    phi = MOD - 1
    m = n_target // k_target
    m_exp = m % phi
    lam = (mod_pow(2, m_exp) - 2 + MOD) % MOD

    fact = [1]*(k_target+1)
    for i in range(1, k_target+1): fact[i] = fact[i-1]*i % MOD
    inv_fact = [1]*(k_target+1)
    inv_fact[k_target] = mod_inv(fact[k_target])
    for i in range(k_target-1, -1, -1): inv_fact[i] = inv_fact[i+1]*(i+1) % MOD

    # NTT via 3 NTT primes + CRT
    P1, R1, RPW1 = 998244353, 3, 1<<23
    P2, R2, RPW2 = 1004535809, 3, 1<<21
    P3, R3, RPW3 = 469762049, 3, 1<<26

    def ntt(a, inv_flag, p, root):
        n = len(a)
        j = 0
        for i in range(1, n):
            bit = n >> 1
            while j & bit: j ^= bit; bit >>= 1
            j ^= bit
            if i < j: a[i], a[j] = a[j], a[i]
        ln = 2
        while ln <= n:
            wl = pow(root, (p-1)//ln, p)
            if inv_flag: wl = pow(wl, p-2, p)
            for i in range(0, n, ln):
                w = 1
                for jj in range(ln//2):
                    u = a[i+jj]; v = a[i+jj+ln//2]*w % p
                    a[i+jj] = (u+v) % p; a[i+jj+ln//2] = (u-v) % p
                    w = w*wl % p
            ln <<= 1
        if inv_flag:
            ni = pow(n, p-2, p)
            for i in range(n): a[i] = a[i]*ni % p

    def conv(a, b, p, root):
        n = 1
        while n < len(a)+len(b): n <<= 1
        fa = list(a) + [0]*(n-len(a)); fb = list(b) + [0]*(n-len(b))
        ntt(fa, False, p, root); ntt(fb, False, p, root)
        for i in range(n): fa[i] = fa[i]*fb[i] % p
        ntt(fa, True, p, root)
        return fa

    def crt3(a1, a2, a3):
        m1, m2, m3 = P1, P2, P3; M = m1*m2
        inv_m1_m2 = pow(m1, m2-2, m2)
        x = a1 + (a2-a1) % m2 * inv_m1_m2 % m2 * m1
        Mm3 = M % m3; inv_M = pow(Mm3, m3-2, m3)
        k = (a3 - x%m3) % m3 * inv_M % m3
        return (x + k*M) % MOD

    def conv_crt(a, b):
        r1 = conv(a, b, P1, R1); r2 = conv(a, b, P2, R2); r3 = conv(a, b, P3, R3)
        return [crt3(r1[i], r2[i], r3[i]) for i in range(len(r1))]

    # H[i] = invFact[i]^16
    H = [mod_pow(inv_fact[i], 16) for i in range(k_target+1)]
    H2 = conv_crt(H, H)

    U = [0]*(k_target+1)
    for s in range(k_target+1):
        A_s = mod_pow(fact[s], 16) * H2[s] % MOD
        U[s] = A_s * inv_fact[s] % MOD

    V = [0]*(k_target+1); lp = 1
    for j in range(k_target+1):
        V[j] = lp * inv_fact[j] % MOD; lp = lp * lam % MOD

    W = conv_crt(U, V)
    return str(fact[k_target] * W[k_target] % MOD)

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

Java

import java.util.Arrays;

public class Euler767 {
    static final long TARGET_MOD = 1000000007L;

    static long power(long base, long exp, long mod) {
        long res = 1;
        base %= mod;
        while (exp > 0) {
            if (exp % 2 == 1) {
                // To avoid BigInteger, we could use modular multiplication
                // But since mod is ~10^9, base * res can be ~10^18, which fits in signed 64-bit
                // long (up to 9 * 10^18)
                res = (res * base) % mod;
            }
            base = (base * base) % mod;
            exp /= 2;
        }
        return res;
    }

    static long modInverse(long n, long mod) {
        return power(n, mod - 2, mod);
    }

    static class NTTMod {
        long mod;
        long root;

        NTTMod(long mod, long root) {
            this.mod = mod;
            this.root = root;
        }
    }

    static final NTTMod P1 = new NTTMod(998244353L, 3L);
    static final NTTMod P2 = new NTTMod(1004535809L, 3L);
    static final NTTMod P3 = new NTTMod(469762049L, 3L);

    static void ntt(long[] a, boolean invert, NTTMod p) {
        int n = a.length;
        for (int i = 1, j = 0; i < n; i++) {
            int bit = n >> 1;
            for (; (j & bit) != 0; bit >>= 1)
                j ^= bit;
            j ^= bit;
            if (i < j) {
                long temp = a[i];
                a[i] = a[j];
                a[j] = temp;
            }
        }

        for (int len = 2; len <= n; len <<= 1) {
            long wlen = power(p.root, (p.mod - 1) / len, p.mod);
            if (invert)
                wlen = modInverse(wlen, p.mod);
            int half = len / 2;
            for (int i = 0; i < n; i += len) {
                long w = 1;
                for (int j = 0; j < half; j++) {
                    long u = a[i + j];
                    long v = (a[i + j + half] * w) % p.mod;
                    a[i + j] = (u + v < p.mod) ? (u + v) : (u + v - p.mod);
                    a[i + j + half] = (u - v >= 0) ? (u - v) : (u - v + p.mod);
                    w = (w * wlen) % p.mod;
                }
            }
        }

        if (invert) {
            long nInv = modInverse(n, p.mod);
            for (int i = 0; i < n; i++) {
                a[i] = (a[i] * nInv) % p.mod;
            }
        }
    }

    static long[] convolve(long[] a, long[] b, NTTMod p) {
        int n = 1;
        while (n < a.length + b.length)
            n <<= 1;
        long[] fa = Arrays.copyOf(a, n);
        long[] fb = Arrays.copyOf(b, n);

        ntt(fa, false, p);
        ntt(fb, false, p);
        for (int i = 0; i < n; i++) {
            fa[i] = (fa[i] * fb[i]) % p.mod;
        }
        ntt(fa, true, p);
        return fa;
    }

    static long crtCombine(long a1, long a2, long a3) {
        long m1 = P1.mod, m2 = P2.mod, m3 = P3.mod;
        long M = m1 * m2;
        long invM1M2 = modInverse(m1, m2);

        long diff = (a2 - a1) % m2;
        if (diff < 0)
            diff += m2;
        long term1 = (diff * invM1M2) % m2;
        long x = a1 + term1 * m1;

        long MModM3 = M % m3;
        long invMM3 = modInverse(MModM3, m3);

        long diff2 = (a3 - (x % m3)) % m3;
        if (diff2 < 0)
            diff2 += m3;
        long k = (diff2 * invMM3) % m3;

        // x + k * M % TARGET_MOD
        // Note: x < M, so x fits in long
        // k * M can be up to m3 * M ~ 4.7 * 10^8 * 10^18 ~ 4.7 * 10^26, exceeds long!
        // We must compute using BigInteger or split M.
        // Or simply (x % TARGET_MOD + k * (M % TARGET_MOD)) % TARGET_MOD
        long xMod = x % TARGET_MOD;
        long MMod = M % TARGET_MOD;
        long res = (xMod + (k % TARGET_MOD) * MMod) % TARGET_MOD;
        return res;
    }

    static long[] convolutionCrt(long[] A, long[] B) {
        long[] r1 = convolve(A, B, P1);
        long[] r2 = convolve(A, B, P2);
        long[] r3 = convolve(A, B, P3);

        int n = r1.length;
        long[] res = new long[n];
        for (int i = 0; i < n; i++) {
            res[i] = crtCombine(r1[i], r2[i], r3[i]);
        }
        return res;
    }

    static class Solver {
        long[] fact, invFact;
        int maxK;

        Solver(int k) {
            maxK = k;
            fact = new long[k + 1];
            invFact = new long[k + 1];
            fact[0] = 1;
            for (int i = 1; i <= k; i++)
                fact[i] = (fact[i - 1] * i) % TARGET_MOD;
            invFact[k] = modInverse(fact[k], TARGET_MOD);
            for (int i = k - 1; i >= 0; i--)
                invFact[i] = (invFact[i + 1] * (i + 1)) % TARGET_MOD;
        }

        long solve(long n) {
            int k = maxK;
            long phi = TARGET_MOD - 1;
            long mExp = (n / k) % phi;

            long lambda = (power(2, mExp, TARGET_MOD) - 2 + TARGET_MOD) % TARGET_MOD;

            long[] H = new long[k + 1];
            for (int i = 0; i <= k; ++i) {
                H[i] = power(invFact[i], 16, TARGET_MOD);
            }

            long[] H2 = convolutionCrt(H, H);

            long[] U = new long[k + 1];
            for (int s = 0; s <= k; ++s) {
                long as = (power(fact[s], 16, TARGET_MOD) * H2[s]) % TARGET_MOD;
                U[s] = (as * invFact[s]) % TARGET_MOD;
            }

            long[] V = new long[k + 1];
            long lamPow = 1;
            for (int j = 0; j <= k; ++j) {
                V[j] = (lamPow * invFact[j]) % TARGET_MOD;
                lamPow = (lamPow * lambda) % TARGET_MOD;
            }

            long[] W = convolutionCrt(U, V);

            return (fact[k] * W[k]) % TARGET_MOD;
        }
    }

    public static String solve() {
        Solver s = new Solver(100000);
        return Long.toString(s.solve(10000000000000000L));
    }

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