Problem 798: Card Stacking Game

Project Euler Problem 798 Solution

EulerSolve provides an optimized solution for Project Euler Problem 798, Card Stacking Game, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The quantity $C(n,s)$ counts losing positions in the card-stacking game with $s$ independent suits and $n$ cards per suit. The real target is $C(10^7,10^7)$ modulo $M=10^9+7$, so brute force is impossible. The implementation therefore does not enumerate positions directly: it first computes the Grundy-value distribution for a single suit, then combines $s$ identical suit distributions by xor convolution. Mathematical Approach Let $a_g$ be the number of one-suit positions whose Grundy value is $g$. The implementation shows that only $0\le g\le n-1$ occur, and that these frequencies partition all one-suit states: $$\sum_{g=0}^{n-1} a_g = 2^n.$$ Once this one-suit distribution is known, the multi-suit count follows from standard impartial-game xor rules. Step 1: Reduce the Full Game to One-Suit Grundy Values If the $j$-th suit contributes Grundy value $x_j$, then the total position has value $$x_1\oplus x_2\oplus \cdots \oplus x_s.$$ A position is losing exactly when this xor is $0$. Therefore $$C(n,s)=\sum_{x_1\oplus\cdots\oplus x_s=0}\prod_{j=1}^{s} a_{x_j}.$$ This formula says that the whole problem is controlled by the frequency table $a_g$ for one suit. Step 2: Rewrite the Linear-Time Recurrence in Closed Form The one-suit distribution is generated by a recurrence indexed by $m=\lfloor g/2\rfloor$....

Detailed mathematical approach

Problem Summary

The quantity $C(n,s)$ counts losing positions in the card-stacking game with $s$ independent suits and $n$ cards per suit. The real target is $C(10^7,10^7)$ modulo $M=10^9+7$, so brute force is impossible. The implementation therefore does not enumerate positions directly: it first computes the Grundy-value distribution for a single suit, then combines $s$ identical suit distributions by xor convolution.

Mathematical Approach

Let $a_g$ be the number of one-suit positions whose Grundy value is $g$. The implementation shows that only $0\le g\le n-1$ occur, and that these frequencies partition all one-suit states:

$$\sum_{g=0}^{n-1} a_g = 2^n.$$

Once this one-suit distribution is known, the multi-suit count follows from standard impartial-game xor rules.

Step 1: Reduce the Full Game to One-Suit Grundy Values

If the $j$-th suit contributes Grundy value $x_j$, then the total position has value

$$x_1\oplus x_2\oplus \cdots \oplus x_s.$$

A position is losing exactly when this xor is $0$. Therefore

$$C(n,s)=\sum_{x_1\oplus\cdots\oplus x_s=0}\prod_{j=1}^{s} a_{x_j}.$$

This formula says that the whole problem is controlled by the frequency table $a_g$ for one suit.

Step 2: Rewrite the Linear-Time Recurrence in Closed Form

The one-suit distribution is generated by a recurrence indexed by $m=\lfloor g/2\rfloor$. After simplifying the ratio updates used by the implementation, the auxiliary quantities become ordinary binomial coefficients:

$$\beta_m=\binom{n-m-2}{m},\qquad \alpha_m=\begin{cases} 0,&m=0,\\ \binom{n-m-2}{m-1},&m\ge 1. \end{cases}$$

Applying Pascal's identity gives two more useful expressions:

$$\pi_m=\alpha_m+\beta_m=\binom{n-m-1}{m},$$

$$\chi_m=\pi_m\frac{n-2m-1}{m+1}=\binom{n-m-1}{m+1}.$$

The power-of-two baseline that appears in the code starts at $2^{n-2}$ and is halved at every step, so at stage $m$ it is exactly

$$2^{n-2-m}.$$

Step 3: Recover the One-Suit Distribution

The remaining correction term is another linear recurrence:

$$\tau_0=1,\qquad \tau_{m+1}=\frac{\tau_m+\alpha_{m+1}}{2}+\beta_{m+1}.$$

With this notation the one-suit frequencies are

$$a_0=2^{n-2}+2,$$

$$a_{2m+1}=2^{n-2-m}+\chi_m-\tau_m \qquad \left(0\le m\le \left\lfloor\frac{n-2}{2}\right\rfloor\right),$$

$$a_{2m}=2^{n-2-m}+\pi_m+\beta_m-\tau_m \qquad \left(1\le m\le \left\lfloor\frac{n-1}{2}\right\rfloor\right).$$

This produces every nonzero $a_g$ in $O(n)$ time, without traversing any of the $2^n$ one-suit states individually.

Step 4: Use Structural Identities as Sanity Checks

The implementation verifies two identities that are mathematically informative and computationally useful:

$$\sum_{g=0}^{n-1} a_g = 2^n,$$

$$a_{n-1}=1.$$

The first identity confirms that the recurrence accounts for every one-suit position exactly once. The second shows that the largest attainable Grundy value appears exactly once. Both checks are strong consistency tests before the xor-convolution stage starts.

Step 5: Combine the $s$ Suits with XOR Convolution

For arrays $f$ and $g$, define xor convolution by

$$ (f *_\oplus g)(t)=\sum_{x\oplus y=t} f(x)g(y). $$

Then $C(n,s)$ is the coefficient at $t=0$ of the $s$-fold xor convolution of the one-suit distribution with itself. Let

$$P=2^{\lceil \log_2 n\rceil}.$$

Pad $(a_0,\dots,a_{n-1})$ with zeros to length $P$, apply the Walsh-Hadamard transform, and write the transformed vector as $\widehat{a}$. The xor-convolution theorem yields

$$C(n,s)=\frac{1}{P}\sum_{i=0}^{P-1}\widehat{a}_i^s \pmod{M},\qquad M=10^9+7.$$

So the $s$-suit aggregation becomes a pointwise power in transform space.

Worked Example: $n=3,\ s=2$

For one suit, the base value is

$$a_0=2^{1}+2=4.$$

At $m=0$, we have $\beta_0=1$, $\alpha_0=0$, $\pi_0=1$, $\chi_0=2$, and $\tau_0=1$. Therefore

$$a_1=2^1+\chi_0-\tau_0=2+2-1=3.$$

Next, $\alpha_1=1$ and $\beta_1=0$, so

$$\tau_1=\frac{1+1}{2}+0=1,$$

and hence

$$a_2=2^0+\pi_1+\beta_1-\tau_1=1+1+0-1=1.$$

Thus the one-suit distribution is

$$a=(4,3,1).$$

Pad it to $P=4$:

$$a=(4,3,1,0).$$

Its Walsh-Hadamard transform is

$$\widehat{a}=(8,2,6,0).$$

For two suits, the losing count is

$$C(3,2)=\frac{8^2+2^2+6^2+0^2}{4}=26,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematics. They first precompute modular inverses up to $n$, so that every rational-looking factor in the recurrence can be evaluated modulo $10^9+7$. They then run one forward pass over $m$, updating the auxiliary sequences and filling the one-suit frequency table for all Grundy values below $n$.

After that, the distribution is padded to the next power of two and transformed in place with the Walsh-Hadamard transform. In the transformed domain, each spectral entry is raised to the $s$-th power, the results are summed, and the final multiplication by $P^{-1}$ extracts the zero-xor coefficient. The C++ and Java implementations also parallelize the large transform and exponentiation loops on big inputs, while the Python implementation delegates to the same compiled arithmetic core, so all three language versions evaluate the same recurrence and the same transform formula.

Complexity Analysis

Building the modular inverses and the one-suit distribution costs $O(n)$ time and $O(n)$ memory before padding. Let $P=2^{\lceil\log_2 n\rceil}$, so $n\le P<2n$. The Walsh-Hadamard transform costs $O(P\log P)$, and exponentiating all transformed entries to the $s$-th power costs $O(P\log s)$. Therefore the total running time is $O(n + P\log P + P\log s)$, which is effectively $O(n\log n + n\log s)$, and the working memory is $O(P)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=798
Sprague-Grundy theorem: Wikipedia - Sprague-Grundy theorem
Nim and xor sums: Wikipedia - Nim
Hadamard transform: Wikipedia - Hadamard transform
Walsh-Hadamard transform for xor convolution: cp-algorithms - Walsh-Hadamard transform

Problem 798 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <future>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <string>
#include <thread>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
#include <pthread.h>
#include <unistd.h>

using namespace std;

static constexpr int MOD = 1'000'000'007;

static inline int addmod(int a, int b) {
    int s = a + b;
    if (s >= MOD) s -= MOD;
    return s;
}

static inline int submod(int a, int b) {
    int s = a - b;
    if (s < 0) s += MOD;
    return s;
}

static long long mod_pow(long long a, long long e) {
    long long r = 1 % MOD;
    a %= MOD;
    while (e > 0) {
        if (e & 1) r = (r * a) % MOD;
        a = (a * a) % MOD;
        e >>= 1;
    }
    return r;
}

static int next_pow2(int x) {
    int p = 1;
    while (p < x) p <<= 1;
    return p;
}

static int get_thread_count() {
    const char* env = getenv("EULER_THREADS");
    if (env && *env) {
        long v = strtol(env, nullptr, 10);
        if (v > 0) return (int)v;
    }
    long n = sysconf(_SC_NPROCESSORS_ONLN);
    if (n < 1) n = 1;
    return (int)n;
}

static vector<int> compute_distribution_one_suit(int n, int L, bool do_checks) {
    vector<int> arr(L, 0);
    if (n <= 0) return arr;

    if (n == 1) {
        arr[0] = 2;
        return arr;
    }

    int maxInv = n;
    vector<int> inv(maxInv + 1);
    inv[1] = 1;
    for (int i = 2; i <= maxInv; i++) {
        inv[i] = int(MOD - (long long)(MOD / i) * inv[MOD % i] % MOD);
    }
    const int inv2 = (MOD + 1) / 2;

    long long pow2N = mod_pow(2, n - 2);
    arr[0] = int((pow2N + 2) % MOD);

    long long sumDist = 0;
    if (do_checks) sumDist = arr[0];

    long long B = 1;
    long long T = 1;

    int m_max = (n - 1) / 2;

    for (int m = 0; m <= m_max; m++) {
        long long denomA = (long long)n - 1 - 2LL * m;
        long long A = 0;
        if (m >= 1) {
            if (denomA == 0) {
                A = 1;
            } else {
                A = (B * m) % MOD * inv[(int)denomA] % MOD;
            }
        }

        long long P = (B + A) % MOD;
        long long Q = P * ((long long)n - 2LL * m - 1) % MOD * inv[m + 1] % MOD;

        int godd = 2 * m + 1;
        if (godd < n) {
            long long val = (pow2N + Q - T) % MOD;
            if (val < 0) val += MOD;
            arr[godd] = (int)val;
            if (do_checks) sumDist += arr[godd];
        }

        int geven = 2 * m;
        if (m >= 1 && geven < n) {
            long long val = (pow2N + P + B - T) % MOD;
            if (val < 0) val += MOD;
            arr[geven] = (int)val;
            if (do_checks) sumDist += arr[geven];
        }

        if (m == m_max) break;

        long long num1 = (long long)n - 2 - 2LL * m;
        long long num2 = (long long)n - 3 - 2LL * m;
        long long den1 = (long long)m + 1;
        long long den2 = (long long)n - 2 - m;
        B = B * (num1 % MOD) % MOD;
        B = B * (num2 % MOD) % MOD;
        B = B * inv[(int)den1] % MOD;
        B = B * inv[(int)den2] % MOD;

        pow2N = pow2N * inv2 % MOD;

        int m1 = m + 1;
        long long denomA_next = (long long)n - 1 - 2LL * m1;
        long long A_next = 0;
        if (denomA_next == 0) A_next = 1;
        else A_next = (B * m1) % MOD * inv[(int)denomA_next] % MOD;

        T = ((T + A_next) % MOD) * inv2 % MOD;
        T = (T + B) % MOD;
    }

    if (do_checks) {
        sumDist %= MOD;
        long long want = mod_pow(2, n);
        if (sumDist != want) {
            cerr << "[CHECK FAILED] sum(a[g]) != 2^n (mod MOD)\n";
            cerr << "sumDist=" << sumDist << " want=" << want << "\n";
            exit(1);
        }
        if (arr[n - 1] != 1) {
            cerr << "[CHECK FAILED] a[n-1] != 1\n";
            cerr << "a[n-1]=" << arr[n - 1] << "\n";
            exit(1);
        }
    }

    return arr;
}

struct WhtTask {
    vector<int>* a;
    int len;
    int start_block;
    int end_block;
};

static void* wht_worker(void* arg) {
    WhtTask* t = (WhtTask*)arg;
    vector<int>& a = *t->a;
    int len = t->len;
    int step = len << 1;
    for (int b = t->start_block; b < t->end_block; b++) {
        int i = b * step;
        for (int j = 0; j < len; j++) {
            int u = a[i + j];
            int v = a[i + j + len];
            a[i + j] = addmod(u, v);
            a[i + j + len] = submod(u, v);
        }
    }
    return nullptr;
}

static void walsh_hadamard_xor(vector<int>& a) {
    const int n = (int)a.size();
    const int threads = get_thread_count();
    for (int len = 1; len < n; len <<= 1) {
        const int step = len << 1;
        const int blocks = n / step;
        if (threads > 1 && blocks >= 2048) {
            int tcnt = min(threads, blocks);
            vector<pthread_t> th(tcnt);
            vector<WhtTask> tasks(tcnt);
            int per = (blocks + tcnt - 1) / tcnt;
            for (int i = 0; i < tcnt; i++) {
                int start = i * per;
                int end = min(blocks, start + per);
                tasks[i] = WhtTask{&a, len, start, end};
                pthread_create(&th[i], nullptr, wht_worker, &tasks[i]);
            }
            for (int i = 0; i < tcnt; i++) {
                pthread_join(th[i], nullptr);
            }
        } else {
            for (int i = 0; i < n; i += step) {
                for (int j = 0; j < len; j++) {
                    int u = a[i + j];
                    int v = a[i + j + len];
                    a[i + j] = addmod(u, v);
                    a[i + j + len] = submod(u, v);
                }
            }
        }
    }
}

struct PowTask {
    const vector<int>* a;
    int s;
    int start;
    int end;
    long long partial;
};

static void* pow_worker(void* arg) {
    PowTask* t = (PowTask*)arg;
    long long sum = 0;
    for (int i = t->start; i < t->end; i++) {
        sum += mod_pow((*t->a)[i], t->s);
    }
    t->partial = sum % MOD;
    return nullptr;
}

static int solve_C(int n, int s, bool do_checks) {
    int L = next_pow2(n);
    vector<int> A = compute_distribution_one_suit(n, L, do_checks);

    walsh_hadamard_xor(A);

    long long total = 0;
    int threads = get_thread_count();
    if (threads > 1 && L >= 4096) {
        int tcnt = min(threads, L);
        vector<pthread_t> th(tcnt);
        vector<PowTask> tasks(tcnt);
        int per = (L + tcnt - 1) / tcnt;
        for (int i = 0; i < tcnt; i++) {
            int start = i * per;
            int end = min(L, start + per);
            tasks[i] = PowTask{&A, s, start, end, 0};
            pthread_create(&th[i], nullptr, pow_worker, &tasks[i]);
        }
        for (int i = 0; i < tcnt; i++) {
            pthread_join(th[i], nullptr);
            total += tasks[i].partial;
        }
        total %= MOD;
    } else {
        for (int i = 0; i < L; i++) {
            total += mod_pow(A[i], s);
        }
        total %= MOD;
    }

    long long invL = mod_pow(L, MOD - 2);
    return int(total * invL % MOD);
}

static void run_selftest() {
    int a = solve_C(3, 2, true);
    int b = solve_C(13, 4, true);
    if (a != 26) {
        cerr << "[SELFTEST FAILED] C(3,2) expected 26 got " << a << "\n";
        exit(1);
    }
    if (b != 540318329) {
        cerr << "[SELFTEST FAILED] C(13,4) expected 540318329 got " << b << "\n";
        exit(1);
    }
    cerr << "[SELFTEST OK]\n";
}

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int n = 10'000'000;
    int s = 10'000'000;
    bool do_checks = true;

    if (argc >= 2) {
        string arg1 = argv[1];
        if (arg1 == "--selftest") {
            run_selftest();
            return 0;
        }
    }
    if (argc >= 3) {
        n = stoi(argv[1]);
        s = stoi(argv[2]);
    }

    int ans = solve_C(n, s, do_checks);
    cout << ans << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler798 {
    static final int MOD = 1000000007;

    static int addmod(int a, int b) {
        int s = a + b;
        if (s >= MOD)
            s -= MOD;
        return s;
    }

    static int submod(int a, int b) {
        int s = a - b;
        if (s < 0)
            s += MOD;
        return s;
    }

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

    static int nextPow2(int x) {
        int p = 1;
        while (p < x)
            p <<= 1;
        return p;
    }

    static int[] computeDistributionOneSuit(int n, int L, boolean doChecks) {
        int[] arr = new int[L];
        if (n <= 0)
            return arr;
        if (n == 1) {
            arr[0] = 2;
            return arr;
        }

        int maxInv = n;
        long[] inv = new long[maxInv + 1];
        inv[1] = 1;
        for (int i = 2; i <= maxInv; i++) {
            inv[i] = MOD - (MOD / i) * inv[MOD % i] % MOD;
        }

        long inv2 = (MOD + 1) / 2;
        long pow2N = modPow(2, n - 2);
        arr[0] = (int) ((pow2N + 2) % MOD);

        long sumDist = 0;
        if (doChecks)
            sumDist = arr[0];

        long B = 1;
        long T = 1;
        int mMax = (n - 1) / 2;

        for (int m = 0; m <= mMax; m++) {
            long denomA = (long) n - 1 - 2L * m;
            long A = 0;
            if (m >= 1) {
                if (denomA == 0)
                    A = 1;
                else
                    A = B * m % MOD * inv[(int) denomA] % MOD;
            }

            long P = (B + A) % MOD;
            long Q = P * (n - 2L * m - 1) % MOD * inv[m + 1] % MOD;

            int godd = 2 * m + 1;
            if (godd < n) {
                long val = (pow2N + Q - T) % MOD;
                if (val < 0)
                    val += MOD;
                arr[godd] = (int) val;
                if (doChecks)
                    sumDist += val;
            }

            int geven = 2 * m;
            if (m >= 1 && geven < n) {
                long val = (pow2N + P + B - T) % MOD;
                if (val < 0)
                    val += MOD;
                arr[geven] = (int) val;
                if (doChecks)
                    sumDist += val;
            }

            if (m == mMax)
                break;

            long num1 = (long) n - 2 - 2L * m;
            long num2 = (long) n - 3 - 2L * m;
            long den1 = m + 1;
            long den2 = n - 2 - m;

            B = B * (num1 % MOD) % MOD;
            B = B * (num2 % MOD) % MOD;
            B = B * inv[(int) den1] % MOD;
            B = B * inv[(int) den2] % MOD;

            pow2N = pow2N * inv2 % MOD;

            int m1 = m + 1;
            long denomANext = (long) n - 1 - 2L * m1;
            long ANext = 0;
            if (denomANext == 0)
                ANext = 1;
            else
                ANext = B * m1 % MOD * inv[(int) denomANext] % MOD;

            T = (T + ANext) % MOD * inv2 % MOD;
            T = (T + B) % MOD;
        }

        if (doChecks) {
            sumDist %= MOD;
            long want = modPow(2, n);
            if (sumDist != want) {
                throw new RuntimeException("CHECK FAILED sumDist");
            }
            if (arr[n - 1] != 1) {
                throw new RuntimeException("CHECK FAILED arr[start_n-1]");
            }
        }

        return arr;
    }

    static void walshHadamardXor(int[] a) {
        int n = a.length;
        int numThreads = Runtime.getRuntime().availableProcessors();
        if (numThreads <= 0)
            numThreads = 1;

        for (int len = 1; len < n; len <<= 1) {
            int step = len << 1;
            int blocks = n / step;
            if (numThreads > 1 && blocks >= 2048) {
                int tcnt = Math.min(numThreads, blocks);
                Thread[] threads = new Thread[tcnt];
                int per = (blocks + tcnt - 1) / tcnt;

                for (int t = 0; t < tcnt; t++) {
                    final int startBlock = t * per;
                    final int endBlock = Math.min(blocks, startBlock + per);
                    final int currentLen = len;
                    final int currentStep = step;

                    threads[t] = new Thread(() -> {
                        for (int b = startBlock; b < endBlock; b++) {
                            int i = b * currentStep;
                            for (int j = 0; j < currentLen; j++) {
                                int u = a[i + j];
                                int v = a[i + j + currentLen];
                                a[i + j] = addmod(u, v);
                                a[i + j + currentLen] = submod(u, v);
                            }
                        }
                    });
                    threads[t].start();
                }
                for (int t = 0; t < tcnt; t++) {
                    try {
                        threads[t].join();
                    } catch (Exception e) {
                    }
                }
            } else {
                for (int i = 0; i < n; i += step) {
                    for (int j = 0; j < len; j++) {
                        int u = a[i + j];
                        int v = a[i + j + len];
                        a[i + j] = addmod(u, v);
                        a[i + j + len] = submod(u, v);
                    }
                }
            }
        }
    }

    static int solveC(int n, int s, boolean doChecks) {
        int L = nextPow2(n);
        int[] A = computeDistributionOneSuit(n, L, doChecks);

        walshHadamardXor(A);

        int numThreads = Runtime.getRuntime().availableProcessors();
        if (numThreads <= 0)
            numThreads = 1;

        long[] partials = new long[numThreads];
        Thread[] threads = new Thread[numThreads];
        int per = (L + numThreads - 1) / numThreads;

        for (int t = 0; t < numThreads; t++) {
            final int tid = t;
            final int start = t * per;
            final int end = Math.min(L, start + per);
            if (start >= end)
                continue;

            threads[t] = new Thread(() -> {
                long sum = 0;
                for (int i = start; i < end; i++) {
                    if (A[i] != 0) {
                        sum += modPow(A[i], s);
                    }
                }
                partials[tid] = sum % MOD;
            });
            threads[t].start();
        }

        long total = 0;
        for (int t = 0; t < numThreads; t++) {
            if (threads[t] != null) {
                try {
                    threads[t].join();
                    total += partials[t];
                } catch (Exception e) {
                }
            }
        }
        total %= MOD;

        long invL = modPow(L, MOD - 2);
        return (int) ((total * invL) % MOD);
    }

    public static String solve() {
        return Integer.toString(solveC(10000000, 10000000, false));
    }

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