Problem 636: Restricted Factorisations

Project Euler Problem 636 Solution

EulerSolve provides an optimized solution for Project Euler Problem 636, Restricted Factorisations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The quantity $F(n)$ counts restricted factorisations of $n!$ of the shape $$n!=b_1\,b_2^2\,b_3^2\,b_4^3\,b_5^3\,b_6^3\,b_7^4\,b_8^4\,b_9^4\,b_{10}^4.$$ The computation is organized in two stages. First, the ten slots are treated as labeled so that equal-base collisions can be handled systematically by inclusion-exclusion. Second, the labels inside the equal-exponent groups are forgotten, because swapping the two square slots, the three cube slots, or the four fourth-power slots does not change the underlying factorisation. The goal is to evaluate $F(10^6)\bmod(10^9+7)$. Mathematical Approach Let the exponent pattern be $$\alpha=(\alpha_1,\dots,\alpha_{10})=(1,2,2,3,3,3,4,4,4,4).$$ If we temporarily label the bases as $b_1,\dots,b_{10}$, then the only real difficulty is that some of these bases may coincide. The implementations resolve that by working on the partition lattice of the ten labels and by treating each prime of $n!$ independently. Step 1: Compute the prime exponents of $n!$ For each prime $p\le n$, Legendre's formula gives $$e_p=v_p(n!)=\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor.$$ Prime powers are independent, so the global count factors over the multiset of values $\{e_p\}$....

Detailed mathematical approach

Problem Summary

The quantity $F(n)$ counts restricted factorisations of $n!$ of the shape

$$n!=b_1\,b_2^2\,b_3^2\,b_4^3\,b_5^3\,b_6^3\,b_7^4\,b_8^4\,b_9^4\,b_{10}^4.$$

The computation is organized in two stages. First, the ten slots are treated as labeled so that equal-base collisions can be handled systematically by inclusion-exclusion. Second, the labels inside the equal-exponent groups are forgotten, because swapping the two square slots, the three cube slots, or the four fourth-power slots does not change the underlying factorisation. The goal is to evaluate $F(10^6)\bmod(10^9+7)$.

Mathematical Approach

Let the exponent pattern be

$$\alpha=(\alpha_1,\dots,\alpha_{10})=(1,2,2,3,3,3,4,4,4,4).$$

If we temporarily label the bases as $b_1,\dots,b_{10}$, then the only real difficulty is that some of these bases may coincide. The implementations resolve that by working on the partition lattice of the ten labels and by treating each prime of $n!$ independently.

Step 1: Compute the prime exponents of $n!$

For each prime $p\le n$, Legendre's formula gives

$$e_p=v_p(n!)=\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor.$$

Prime powers are independent, so the global count factors over the multiset of values $\{e_p\}$. It is therefore enough to store the histogram

$$c_e=\#\{p\le n:\ v_p(n!)=e\}.$$

Once the numbers $c_e$ are known, every partition profile contributes through repeated powers of the same local coefficient $g(e)$.

Step 2: Encode equal-base collisions by set partitions

Suppose several labeled slots share the same base. Their collision pattern is a set partition $\pi$ of $\{1,\dots,10\}$. Every block $B\in\pi$ represents one actual base, and the exponents attached to that base add up to

$$S_B=\sum_{i\in B}\alpha_i.$$

So after choosing $\pi$, the original labels only matter through the multiset of block sums $\{S_B\}$. Two different partitions that produce the same sorted block-sum vector have the same local generating function, so they can be merged by adding their Möbius weights. The implementations do exactly that.

Step 3: Count one prime with a generating function

Fix a prime $p$ with exponent $e=e_p$. If block $B$ receives $u_B\ge 0$ copies of $p$ in its common base, then the exponent balance is

$$\sum_{B\in\pi} S_Bu_B=e.$$

The number of nonnegative solutions is the coefficient of

$$G_\pi(x)=\prod_{B\in\pi}\frac{1}{1-x^{S_B}}=\sum_{m\ge 0} g_\pi(m)x^m.$$

Therefore $g_\pi(e)$ is exactly the number of ways a single prime of valuation $e$ can be distributed across the bases while respecting the collision pattern $\pi$.

Step 4: Recover distinct labeled bases by Möbius inversion

Counting every partition $\pi$ would allow equal bases. To force the ten labeled bases to be distinct, we use Möbius inversion on the partition lattice. Its weight is

$$\mu(\pi)=(-1)^{10-|\pi|}\prod_{B\in\pi}(|B|-1)!.$$

Because primes are independent, the total labeled contribution of $\pi$ is

$$W_\pi(n)=\prod_{p\le n} g_\pi(e_p)=\prod_{e\ge 1} g_\pi(e)^{c_e}.$$

Hence the number of factorizations with fully labeled and pairwise distinct bases is

$$L(n)=\sum_{\pi}\mu(\pi)\,W_\pi(n).$$

Step 5: Forget relabelings inside equal-exponent groups

The labeled model distinguishes the two exponent-$2$ slots, the three exponent-$3$ slots, and the four exponent-$4$ slots. The original problem does not. Every genuine restricted factorisation is therefore counted

$$2!\,3!\,4!=288$$

times inside $L(n)$. The desired answer is

$$F(n)=\frac{L(n)}{288}\pmod{10^9+7}.$$

Since the modulus is prime, the division is implemented by multiplying with the modular inverse of $288$.

Step 6: Evaluate $g_\pi(e)$ quickly for all relevant exponents

For small exponents, $g_\pi(e)$ is computed directly by unbounded knapsack on the block sums $S_B$. This is exactly coefficient extraction for $G_\pi(x)$.

For large exponents, the denominator of $G_\pi(x)$ has degree at most

$$1+2+2+3+3+3+4+4+4+4=30,$$

so if

$$Q_\pi(x)=\prod_{B\in\pi}(1-x^{S_B})=1+q_1x+\cdots+q_{30}x^{30},$$

then the coefficients satisfy the linear recurrence

$$g_\pi(m)=-\sum_{j=1}^{30}q_j\,g_\pi(m-j)\qquad(m\ge 30).$$

This fixed order is the key reason the implementation can jump to distant coefficients quickly instead of extending the dynamic program all the way to $v_2(n!)$.

Worked Example: the partition with no collisions

If every labeled base is distinct, all ten blocks are singletons and the local series becomes

$$G(x)=\frac{1}{(1-x)(1-x^2)^2(1-x^3)^3(1-x^4)^4}.$$

Now read off a few coefficients. We have $g(1)=1$, because an exponent $1$ can only come from the first slot. We also have $g(2)=3$, because exponent $2$ can be formed either as $2\cdot 1$ in the first slot or as one copy of either of the two square slots.

This tiny example illustrates two implementation ideas at once. First, the local coefficients are genuine coin-change counts for the block sums. Second, if a partition profile had all block sums divisible by some $d>1$, then $g_\pi(1)=0$. For the target input $n=10^6$, primes with valuation $1$ do occur in $n!$, so such profiles can be discarded immediately.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical pipeline. They first enumerate set partitions of the ten exponent slots, convert each partition to its sorted block-sum vector, merge profiles with identical vectors, and store the corresponding Möbius weight. They also keep the gcd of each profile so that obviously impossible profiles can be pruned when exponent $1$ appears in the histogram.

Next, the implementation sieves primes up to $n$, applies Legendre's formula to every prime, and compresses the result into the frequency table $c_e$. For each surviving profile it computes all small coefficients $g_\pi(e)$ by dynamic programming, then uses the degree-$30$ recurrence to evaluate larger exponents. The contribution $\mu(\pi)\prod_e g_\pi(e)^{c_e}$ is accumulated modulo $10^9+7$, and the final sum is multiplied by the modular inverse of $288$. The C++ version parallelizes independent profile evaluations, the Java version keeps the same arithmetic sequentially, and the Python version delegates to the same C++ core.

Complexity Analysis

The partition-lattice work depends only on the fixed exponent multiset $(1,2,2,3,3,3,4,4,4,4)$, so it is constant-size preprocessing for this problem. After that, building the prime table up to $n$ costs $O(n\log\log n)$ time and $O(n)$ memory. Let $P$ be the number of distinct block-sum profiles after aggregation; $P$ is again a problem-specific constant. The dynamic-programming phase for each profile only goes up to about $\sqrt n$, while larger coefficients are recovered from an order-$30$ recurrence in logarithmic time. Consequently the overall asymptotic cost is dominated by the sieve, with $O(n\log\log n)$ time and $O(n)$ memory, and the partition-profile work contributes a manageable constant factor.

Footnotes and References

Problem page: https://projecteuler.net/problem=636
Legendre's formula: Wikipedia - Legendre's formula
Set partitions and partition lattices: Wikipedia - Partition of a set
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Generating functions: Wikipedia - Generating function
Linear recurrences with constant coefficients: Wikipedia - Linear recurrence with constant coefficients

Problem 636 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <map>
#include <numeric>
#include <thread>
#include <vector>

using namespace std;

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

// Count labeled tuples using inclusion-exclusion on equal bases; each partition
// yields block sums S and a generating function 1/∏(1 - x^S). Then divide by
// 2!3!4! to forget permutations of identical exponents.

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 int mod_pow(int base, long long exp) {
    long long res = 1;
    long long b = base % MOD;
    while (exp > 0) {
        if (exp & 1LL) res = (res * b) % MOD;
        b = (b * b) % MOD;
        exp >>= 1LL;
    }
    return static_cast<int>(res);
}

struct Entry {
    vector<int> sums;
    int weight_mod;
    int gcd;
};

static vector<Entry> build_entries() {
    const vector<int> exps = {1, 2, 2, 3, 3, 3, 4, 4, 4, 4};
    const int n = static_cast<int>(exps.size());

    long long fact[11];
    fact[0] = 1;
    for (int i = 1; i <= 10; ++i) fact[i] = fact[i - 1] * i;

    map<vector<int>, long long> weight_map;
    vector<int> sums;
    vector<int> sizes;

    auto dfs = [&](auto&& self, int idx) -> void {
        if (idx == n) {
            vector<int> key = sums;
            sort(key.begin(), key.end());
            int blocks = static_cast<int>(sizes.size());
            long long mu = ((n - blocks) % 2 == 0) ? 1 : -1;
            for (int sz : sizes) mu *= fact[sz - 1];
            weight_map[key] += mu;
            return;
        }

        for (size_t i = 0; i < sums.size(); ++i) {
            sums[i] += exps[idx];
            sizes[i] += 1;
            self(self, idx + 1);
            sizes[i] -= 1;
            sums[i] -= exps[idx];
        }

        sums.push_back(exps[idx]);
        sizes.push_back(1);
        self(self, idx + 1);
        sums.pop_back();
        sizes.pop_back();
    };

    dfs(dfs, 0);

    vector<Entry> entries;
    entries.reserve(weight_map.size());
    for (const auto& kv : weight_map) {
        long long w = kv.second;
        if (w == 0) continue;
        int g = 0;
        for (int s : kv.first) g = std::gcd(g, s);
        int wmod = static_cast<int>(w % MOD);
        if (wmod < 0) wmod += MOD;
        entries.push_back({kv.first, wmod, g});
    }
    return entries;
}

static const vector<Entry>& get_entries() {
    static vector<Entry> entries = build_entries();
    return entries;
}

static vector<int> sieve_primes(int n) {
    vector<int> primes;
    if (n < 2) return primes;
    vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i <= n; ++i) {
        if (!is_prime[i]) continue;
        primes.push_back(i);
        if (static_cast<long long>(i) * i <= n) {
            for (int j = i * i; j <= n; j += i) is_prime[j] = false;
        }
    }
    return primes;
}

struct ExpCounts {
    vector<pair<int, int>> small;
    vector<pair<int, int>> big;
    int maxSmall;
    bool hasExp1;
};

static ExpCounts build_exp_counts(int n) {
    vector<int> primes = sieve_primes(n);
    vector<int> cnt(n + 1, 0);
    for (int p : primes) {
        int e = 0;
        int m = n;
        while (m > 0) {
            m /= p;
            e += m;
        }
        cnt[e] += 1;
    }

    int maxSmall = max(30, static_cast<int>(sqrt(static_cast<double>(n))));

    vector<pair<int, int>> small;
    vector<pair<int, int>> big;
    for (int e = 1; e <= n; ++e) {
        if (cnt[e] == 0) continue;
        if (e <= maxSmall) {
            small.emplace_back(e, cnt[e]);
        } else {
            big.emplace_back(e, cnt[e]);
        }
    }

    return {small, big, maxSmall, (n >= 2 ? cnt[1] > 0 : false)};
}

using Vec = array<int, K>;

static Vec combine_vecs(const Vec& a, const Vec& b, const Vec& rec) {
    long long tmp[2 * K] = {};
    for (int i = 0; i < K; ++i) {
        if (a[i] == 0) continue;
        for (int j = 0; j < K; ++j) {
            if (b[j] == 0) continue;
            tmp[i + j] = (tmp[i + j] + 1LL * a[i] * b[j]) % MOD;
        }
    }

    for (int i = 2 * K - 2; i >= K; --i) {
        long long val = tmp[i];
        if (val == 0) continue;
        for (int j = 1; j <= K; ++j) {
            tmp[i - j] = (tmp[i - j] + val * rec[j - 1]) % MOD;
        }
    }

    Vec res{};
    for (int i = 0; i < K; ++i) res[i] = static_cast<int>(tmp[i]);
    return res;
}

static int linear_rec(const Vec& rec, const int* init, long long n) {
    if (n < K) return init[n];
    Vec pol{};
    Vec base{};
    pol[0] = 1;
    base[1] = 1; // x

    long long m = n;
    while (m > 0) {
        if (m & 1LL) pol = combine_vecs(pol, base, rec);
        base = combine_vecs(base, base, rec);
        m >>= 1LL;
    }

    long long ans = 0;
    for (int i = 0; i < K; ++i) {
        ans = (ans + 1LL * pol[i] * init[i]) % MOD;
    }
    return static_cast<int>(ans);
}

static int compute_contribution(const Entry& entry,
                                const vector<pair<int, int>>& small,
                                const vector<pair<int, int>>& big,
                                int maxSmall,
                                vector<int>& dp) {
    // Unbounded knapsack for g(e) up to sqrt(n); bigger exponents use the recurrence.
    fill(dp.begin(), dp.end(), 0);
    dp[0] = 1;
    for (int s : entry.sums) {
        for (int e = s; e <= maxSmall; ++e) {
            dp[e] = addmod(dp[e], dp[e - s]);
        }
    }

    long long res = 1;
    for (const auto& pr : small) {
        int e = pr.first;
        int cnt = pr.second;
        int val = dp[e];
        if (val == 0) return 0;
        res = (res * mod_pow(val, cnt)) % MOD;
    }

    if (big.empty()) {
        return static_cast<int>(res * entry.weight_mod % MOD);
    }

    array<int, K + 1> C{};
    C[0] = 1;
    int deg = 0;
    for (int s : entry.sums) {
        for (int i = deg; i >= 0; --i) {
            C[i + s] = submod(C[i + s], C[i]);
        }
        deg += s;
    }

    Vec rec{};
    for (int i = 1; i <= K; ++i) {
        rec[i - 1] = submod(0, C[i]);
    }

    const int* init = dp.data();
    for (const auto& pr : big) {
        int e = pr.first;
        int cnt = pr.second;
        int val = linear_rec(rec, init, e);
        if (val == 0) return 0;
        res = (res * mod_pow(val, cnt)) % MOD;
    }

    return static_cast<int>(res * entry.weight_mod % MOD);
}

static int compute_F(int n, unsigned threads) {
    ExpCounts counts = build_exp_counts(n);
    const auto& entries = get_entries();

    vector<int> active;
    active.reserve(entries.size());
    if (counts.hasExp1) {
        for (size_t i = 0; i < entries.size(); ++i) {
            if (entries[i].gcd == 1) active.push_back(static_cast<int>(i));
        }
    } else {
        for (size_t i = 0; i < entries.size(); ++i) {
            active.push_back(static_cast<int>(i));
        }
    }

    if (active.empty()) return 0;

    unsigned tcount = threads == 0 ? 1u : threads;
    if (tcount > active.size()) tcount = static_cast<unsigned>(active.size());

    const int inv288 = mod_pow(288, MOD - 2);

    if (tcount <= 1 || active.size() < 64) {
        vector<int> dp(counts.maxSmall + 1);
        long long total = 0;
        for (int idx : active) {
            total += compute_contribution(entries[idx], counts.small, counts.big, counts.maxSmall, dp);
            if (total >= MOD) total -= MOD;
        }
        return static_cast<int>(total * inv288 % MOD);
    }

    atomic<size_t> next(0);
    vector<long long> partial(tcount, 0);
    vector<thread> pool;
    pool.reserve(tcount);

    for (unsigned t = 0; t < tcount; ++t) {
        pool.emplace_back([&, t]() {
            vector<int> dp(counts.maxSmall + 1);
            long long local = 0;
            while (true) {
                size_t pos = next.fetch_add(1);
                if (pos >= active.size()) break;
                int idx = active[pos];
                local += compute_contribution(entries[idx], counts.small, counts.big, counts.maxSmall, dp);
                if (local >= MOD) local -= MOD;
            }
            partial[t] = local % MOD;
        });
    }

    for (auto& th : pool) th.join();

    long long total = 0;
    for (long long v : partial) {
        total += v;
        if (total >= MOD) total -= MOD;
    }

    return static_cast<int>(total * inv288 % MOD);
}

static bool run_validation(unsigned threads) {
    struct Case {
        int n;
        int expected;
    };

    const Case cases[] = {
        {25, 4933},
        {100, 693952493},
        {1000, 6364496},
    };

    bool ok = true;
    for (const auto& c : cases) {
        int got = compute_F(c.n, threads);
        if (got != c.expected) {
            cerr << "Validation failed for n=" << c.n
                 << ": got " << got << ", expected " << c.expected << "\n";
            ok = false;
        }
    }

    if (ok) {
        cerr << "Validation checkpoints passed.\n";
    }
    return ok;
}

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

    int n = 1'000'000;
    unsigned threads = thread::hardware_concurrency();
    if (threads == 0) threads = 1;
    threads = min(threads, 8u);
    bool validate = true;

    // Optional CLI: ./a.out [n] [threads] [validate(0/1)]
    if (argc >= 2) n = stoi(argv[1]);
    if (argc >= 3) threads = static_cast<unsigned>(stoul(argv[2]));
    if (argc >= 4) validate = (stoi(argv[3]) != 0);

    if (threads == 0) threads = 1;

    if (validate && !run_validation(min(threads, 4u))) {
        return 1;
    }

    cout << compute_F(n, threads) << "\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

import java.util.*;

public class Euler636 {
    static final int MOD = 1000000007;
    static final int K = 30;

    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 int modPow(long base, long exp) {
        long res = 1;
        long b = base % MOD;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = (res * b) % MOD;
            b = (b * b) % MOD;
            exp >>= 1;
        }
        return (int) res;
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static class Entry {
        ArrayList<Integer> sums;
        int weightMod;
        int gcdVal;

        Entry(ArrayList<Integer> sums, int weightMod, int gcdVal) {
            this.sums = sums;
            this.weightMod = weightMod;
            this.gcdVal = gcdVal;
        }
    }

    static ArrayList<Entry> entriesCache = null;

    static ArrayList<Entry> buildEntries() {
        if (entriesCache != null)
            return entriesCache;
        int[] exps = { 1, 2, 2, 3, 3, 3, 4, 4, 4, 4 };
        int n = exps.length;
        long[] fact = new long[11];
        fact[0] = 1;
        for (int i = 1; i <= 10; ++i)
            fact[i] = fact[i - 1] * i;

        HashMap<ArrayList<Integer>, Long> weightMap = new HashMap<>();
        ArrayList<Integer> sums = new ArrayList<>();
        ArrayList<Integer> sizes = new ArrayList<>();

        class Dfs {
            void run(int idx) {
                if (idx == n) {
                    ArrayList<Integer> key = new ArrayList<>(sums);
                    Collections.sort(key);
                    int blocks = sizes.size();
                    long mu = ((n - blocks) % 2 == 0) ? 1 : -1;
                    for (int sz : sizes)
                        mu *= fact[sz - 1];
                    weightMap.put(key, weightMap.getOrDefault(key, 0L) + mu);
                    return;
                }
                for (int i = 0; i < sums.size(); ++i) {
                    sums.set(i, sums.get(i) + exps[idx]);
                    sizes.set(i, sizes.get(i) + 1);
                    run(idx + 1);
                    sizes.set(i, sizes.get(i) - 1);
                    sums.set(i, sums.get(i) - exps[idx]);
                }
                sums.add(exps[idx]);
                sizes.add(1);
                run(idx + 1);
                sums.remove(sums.size() - 1);
                sizes.remove(sizes.size() - 1);
            }
        }
        new Dfs().run(0);

        ArrayList<Entry> entries = new ArrayList<>();
        for (Map.Entry<ArrayList<Integer>, Long> kv : weightMap.entrySet()) {
            long w = kv.getValue();
            if (w == 0)
                continue;
            int g = 0;
            for (int s : kv.getKey())
                g = gcd(g, s);
            int wmod = (int) (w % MOD);
            if (wmod < 0)
                wmod += MOD;
            entries.add(new Entry(kv.getKey(), wmod, g));
        }
        entriesCache = entries;
        return entries;
    }

    static ArrayList<Integer> sievePrimes(int n) {
        ArrayList<Integer> primes = new ArrayList<>();
        if (n < 2)
            return primes;
        boolean[] isPrime = new boolean[n + 1];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int i = 2; i <= n; ++i) {
            if (!isPrime[i])
                continue;
            primes.add(i);
            if ((long) i * i <= n) {
                for (int j = i * i; j <= n; j += i)
                    isPrime[j] = false;
            }
        }
        return primes;
    }

    static class ExpCounts {
        ArrayList<int[]> small = new ArrayList<>();
        ArrayList<int[]> big = new ArrayList<>();
        int maxSmall;
        boolean hasExp1;
    }

    static ExpCounts buildExpCounts(int n) {
        ArrayList<Integer> primes = sievePrimes(n);
        int[] cnt = new int[n + 1];
        for (int p : primes) {
            int e = 0;
            int m = n;
            while (m > 0) {
                m /= p;
                e += m;
            }
            cnt[e]++;
        }

        int maxSmall = Math.max(30, (int) Math.sqrt(n));
        ExpCounts res = new ExpCounts();
        res.maxSmall = maxSmall;
        res.hasExp1 = n >= 2 && cnt[1] > 0;

        for (int e = 1; e <= n; ++e) {
            if (cnt[e] == 0)
                continue;
            if (e <= maxSmall)
                res.small.add(new int[] { e, cnt[e] });
            else
                res.big.add(new int[] { e, cnt[e] });
        }
        return res;
    }

    static int[] combineVecs(int[] a, int[] b, int[] rec) {
        long[] tmp = new long[2 * K];
        for (int i = 0; i < K; ++i) {
            if (a[i] == 0)
                continue;
            for (int j = 0; j < K; ++j) {
                if (b[j] == 0)
                    continue;
                tmp[i + j] = (tmp[i + j] + (long) a[i] * b[j]) % MOD;
            }
        }
        for (int i = 2 * K - 2; i >= K; --i) {
            long val = tmp[i];
            if (val == 0)
                continue;
            for (int j = 1; j <= K; ++j) {
                tmp[i - j] = (tmp[i - j] + val * rec[j - 1]) % MOD;
            }
        }
        int[] res = new int[K];
        for (int i = 0; i < K; ++i)
            res[i] = (int) tmp[i];
        return res;
    }

    static int linearRec(int[] rec, int[] init, int n) {
        if (n < K)
            return init[n];
        int[] pol = new int[K];
        int[] base = new int[K];
        pol[0] = 1;
        base[1] = 1;

        int m = n;
        while (m > 0) {
            if ((m & 1) == 1)
                pol = combineVecs(pol, base, rec);
            base = combineVecs(base, base, rec);
            m >>= 1;
        }

        long ans = 0;
        for (int i = 0; i < K; ++i) {
            ans = (ans + (long) pol[i] * init[i]) % MOD;
        }
        return (int) ans;
    }

    static int computeContribution(Entry entry, ExpCounts counts, int[] dp) {
        Arrays.fill(dp, 0);
        dp[0] = 1;
        for (int s : entry.sums) {
            for (int e = s; e <= counts.maxSmall; ++e) {
                dp[e] = addmod(dp[e], dp[e - s]);
            }
        }

        long res = 1;
        for (int[] pr : counts.small) {
            int val = dp[pr[0]];
            if (val == 0)
                return 0;
            res = (res * modPow(val, pr[1])) % MOD;
        }

        if (counts.big.isEmpty()) {
            return (int) (res * entry.weightMod % MOD);
        }

        int[] C = new int[K + 1];
        C[0] = 1;
        int deg = 0;
        for (int s : entry.sums) {
            for (int i = deg; i >= 0; --i) {
                C[i + s] = submod(C[i + s], C[i]);
            }
            deg += s;
        }

        int[] rec = new int[K];
        for (int i = 1; i <= K; ++i) {
            rec[i - 1] = submod(0, C[i]);
        }

        for (int[] pr : counts.big) {
            int val = linearRec(rec, dp, pr[0]);
            if (val == 0)
                return 0;
            res = (res * modPow(val, pr[1])) % MOD;
        }

        return (int) (res * entry.weightMod % MOD);
    }

    public static String solve() {
        int n = 1000000;
        ExpCounts counts = buildExpCounts(n);
        ArrayList<Entry> entries = buildEntries();

        ArrayList<Integer> active = new ArrayList<>();
        if (counts.hasExp1) {
            for (int i = 0; i < entries.size(); ++i) {
                if (entries.get(i).gcdVal == 1)
                    active.add(i);
            }
        } else {
            for (int i = 0; i < entries.size(); ++i)
                active.add(i);
        }

        if (active.isEmpty())
            return "0";

        int inv288 = modPow(288, MOD - 2);
        int[] dp = new int[counts.maxSmall + 1];
        long total = 0;

        for (int idx : active) {
            total += computeContribution(entries.get(idx), counts, dp);
            if (total >= MOD)
                total -= MOD;
        }

        return Integer.toString((int) (total * inv288 % MOD));
    }

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