Problem 598: Split Divisibilities

Project Euler Problem 598 Solution

EulerSolve provides an optimized solution for Project Euler Problem 598, Split Divisibilities, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We study factor pairs of $100!$. For every pair $(a,b)$ with $ab=100!$ and $a\le b$, we ask whether $a$ and $b$ have the same number of positive divisors. If $\tau$ denotes the divisor-counting function, the required quantity is $$C(100!)=\#\left\{(a,b):ab=100!,\ a\le b,\ \tau(a)=\tau(b)\right\}.$$ A naive scan over all divisors of $100!$ is far too large. The successful approach works with prime exponents and balances the prime factors of $\tau(a)$ and $\tau(b)$ instead of comparing $a$ and $b$ directly. Mathematical Approach Write the prime factorization of $100!$ as $$100!=\prod_{p\le 100} p^{e_p},\qquad e_p=\sum_{k\ge 1}\left\lfloor\frac{100}{p^k}\right\rfloor.$$ Each prime $p$ contributes an exponent $e_p$, and every factor pair of $100!$ comes from deciding how much of that exponent belongs to $a$ and how much belongs to $b$. Step 1: Encode Every Factor Pair by Exponent Choices Choose integers $x_p$ with $0\le x_p\le e_p$ and define $$a=\prod_{p\le 100} p^{x_p},\qquad b=\prod_{p\le 100} p^{e_p-x_p}.$$ Every ordered factorization $ab=100!$ arises exactly once in this way. The divisor-count formulas become $$\tau(a)=\prod_{p\le 100}(x_p+1),\qquad \tau(b)=\prod_{p\le 100}(e_p-x_p+1).$$ So the problem reduces to choosing one integer from each interval $0,\dots,e_p$ so that these two products are equal....

Detailed mathematical approach

Problem Summary

We study factor pairs of $100!$. For every pair $(a,b)$ with $ab=100!$ and $a\le b$, we ask whether $a$ and $b$ have the same number of positive divisors. If $\tau$ denotes the divisor-counting function, the required quantity is

$$C(100!)=\#\left\{(a,b):ab=100!,\ a\le b,\ \tau(a)=\tau(b)\right\}.$$

A naive scan over all divisors of $100!$ is far too large. The successful approach works with prime exponents and balances the prime factors of $\tau(a)$ and $\tau(b)$ instead of comparing $a$ and $b$ directly.

Mathematical Approach

Write the prime factorization of $100!$ as

$$100!=\prod_{p\le 100} p^{e_p},\qquad e_p=\sum_{k\ge 1}\left\lfloor\frac{100}{p^k}\right\rfloor.$$

Each prime $p$ contributes an exponent $e_p$, and every factor pair of $100!$ comes from deciding how much of that exponent belongs to $a$ and how much belongs to $b$.

Step 1: Encode Every Factor Pair by Exponent Choices

Choose integers $x_p$ with $0\le x_p\le e_p$ and define

$$a=\prod_{p\le 100} p^{x_p},\qquad b=\prod_{p\le 100} p^{e_p-x_p}.$$

Every ordered factorization $ab=100!$ arises exactly once in this way. The divisor-count formulas become

$$\tau(a)=\prod_{p\le 100}(x_p+1),\qquad \tau(b)=\prod_{p\le 100}(e_p-x_p+1).$$

So the problem reduces to choosing one integer from each interval $0,\dots,e_p$ so that these two products are equal.

Step 2: Turn $\tau(a)=\tau(b)$ into Prime-Valuation Balances

Two positive integers are equal if and only if every prime appears with the same exponent in both. Therefore $\tau(a)=\tau(b)$ is equivalent to

$$\sum_{p\le 100}\nu_q(x_p+1)=\sum_{p\le 100}\nu_q(e_p-x_p+1).$$

This identity must hold for every prime $q$.

Define the balance at prime $q$ by

$$D_q=\sum_{p\le 100}\left(\nu_q(x_p+1)-\nu_q(e_p-x_p+1)\right).$$

Then a choice of exponents is valid exactly when

$$D_q=0.$$

This must again hold for every prime $q$.

This reformulation is the core of the method: we count exponent assignments whose valuation-balance vector is identically zero.

Step 3: Why $2,3,5,7$ Are the Exceptional Primes

Legendre's formula gives

$$e_2=97,\qquad e_3=48,\qquad e_5=24,\qquad e_7=16,$$

while every prime $p>7$ satisfies

$$e_p\le 9.$$

Hence, for $p>7$, the numbers $x_p+1$ and $e_p-x_p+1$ always lie between $1$ and $10$. Their prime factors can only belong to

$$\{2,3,5,7\}.$$

That means every prime $p>7$ can affect only the four balances $D_2,D_3,D_5,D_7$. It can never create or cancel a balance at any larger prime $q>7$.

Step 4: Split the Search into Two Independent Parts

The previous observation implies that every balance $D_q$ with $q>7$ must already vanish inside the contribution coming from $p\in\{2,3,5,7\}$. The search therefore splits naturally.

First, enumerate all choices for the four exceptional primes $2,3,5,7$. For each choice, compute the full balance vector $(D_q)$. Since $x_2+1$ can be as large as $98$, this stage may involve primes such as $11,13,\dots,97$. Keep only those assignments for which every coordinate with $q>7$ is already zero.

After that filter, only the four coordinates

$$\left(D_2,D_3,D_5,D_7\right)$$

still matter. Second, process the primes $p>7$. Each such prime contributes one small four-dimensional delta vector, because only $2,3,5,7$ can appear in the corresponding divisor-count factors.

Step 5: Match Complementary Four-Dimensional States

Let $\mathbf d=(d_2,d_3,d_5,d_7)$ be a state produced by the exceptional primes after the large-prime balances have been forced to zero. Let $\mathbf r$ be a state produced by the remaining primes. A complete exponent assignment is valid exactly when

$$\mathbf d+\mathbf r=\mathbf 0.$$

So the ordered answer is obtained by counting, for every state from the exceptional-prime side, how many states from the remaining-prime side are its exact negative.

Step 6: Convert Ordered Solutions to the Requested Unordered Pairs

The construction above counts ordered pairs $(a,b)$. The problem asks for unordered pairs with $a\le b$. Swapping $a$ and $b$ replaces every $x_p$ by $e_p-x_p$, so non-fixed solutions come in symmetric pairs. Therefore

$$C(100!)=\frac{M+F}{2},$$

where $M$ is the ordered count and $F=1$ only if $100!$ is a perfect square. Here $F=0$, because not all exponents $e_p$ are even.

Worked Example: $10!$

The same machinery is easy to inspect on

$$10!=2^8\cdot 3^4\cdot 5^2\cdot 7.$$

In this smaller case there are no primes $p>7$, so only the exceptional stage remains. The primes $5$ and $7$ together can produce the six balances

$$(-1,-1,0,0),\ (-1,0,0,0),\ (-1,1,0,0),\ (1,-1,0,0),\ (1,0,0,0),\ (1,1,0,0).$$

The primes $2$ and $3$ produce the opposite balances with multiplicities $2,1,1,2$ on the four vectors

$$(-1,0,0,0),\ (-1,1,0,0),\ (1,-1,0,0),\ (1,0,0,0),$$

and no matches for $(-1,-1,0,0)$ or $(1,1,0,0)$. Hence the ordered count is

$$M=2+1+1+2=6.$$

Since $10!$ is not a square, $F=0$, so

$$C(10!)=\frac{6}{2}=3.$$

This is the small checkpoint reproduced by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first generate all primes up to $100$ and compute every exponent $e_p$ with Legendre's formula. Next they precompute the prime factorizations of all integers from $1$ to $\max(e_p)+1$. That turns every quantity of the form $\nu_q(x_p+1)$ or $\nu_q(e_p-x_p+1)$ into a fast table lookup.

For the exceptional primes $2,3,5,7$, the implementation enumerates every possible complementary pair $(x_p+1,e_p-x_p+1)$, builds the corresponding balance vectors, and merges them in two batches before forming the full four-prime combinations. Only combinations whose coordinates for primes $q>7$ are all zero are kept, and their remaining four-coordinate balances are stored with multiplicities.

For each prime $p>7$, the implementation computes all possible deltas on $(D_2,D_3,D_5,D_7)$ and updates a sparse map from four-dimensional balance states to counts. The ordered total is the sum of products of matching complementary states from the two stages. A final symmetry correction converts that ordered total into the requested count with $a\le b$.

Complexity Analysis

Let

$$A=(e_2+1)(e_3+1)(e_5+1)(e_7+1).$$

Enumerating the exceptional primes costs $O(A\cdot B)$, where $B$ is the number of prime coordinates that may appear in the divisor-count factors. For $100!$, this is practical because only four primes have large exponent ranges.

For the remaining primes, if $S_k$ denotes the number of distinct four-dimensional states after processing the first $k$ primes greater than $7$, then the dynamic program costs

$$O\left(\sum_k S_{k-1}(e_{p_k}+1)\right)$$

time and $O(\max_k S_k)$ memory. Since every prime $p>7$ has $e_p+1\le 10$, the branching factor is tiny and the sparse state space stays manageable.

Footnotes and References

Problem page: https://projecteuler.net/problem=598
Divisor function: Wikipedia - Divisor function
Legendre's formula: Wikipedia - Legendre's formula
Prime factorization: Wikipedia - Prime factorization
$p$-adic valuation: Wikipedia - $p$-adic valuation

Problem 598 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>


using u64 = std::uint64_t;
using u128 = unsigned __int128;

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

static int exp_in_fact(int n, int p) {
    int e = 0;
    while (n) {
        n /= p;
        e += n;
    }
    return e;
}

static u64 pack4(int d2, int d3, int d5, int d7) {
    static constexpr int OFF = 128;
    const u64 a = (u64)(std::uint16_t)(d2 + OFF);
    const u64 b = (u64)(std::uint16_t)(d3 + OFF);
    const u64 c = (u64)(std::uint16_t)(d5 + OFF);
    const u64 d = (u64)(std::uint16_t)(d7 + OFF);
    return a | (b << 16) | (c << 32) | (d << 48);
}

static void unpack4(u64 key, int& d2, int& d3, int& d5, int& d7) {
    static constexpr int OFF = 128;
    d2 = (int)(std::uint16_t)(key & 0xFFFFULL) - OFF;
    d3 = (int)(std::uint16_t)((key >> 16) & 0xFFFFULL) - OFF;
    d5 = (int)(std::uint16_t)((key >> 32) & 0xFFFFULL) - OFF;
    d7 = (int)(std::uint16_t)((key >> 48) & 0xFFFFULL) - OFF;
}

static u128 C_factorial(int n) {
    const auto ps = primes_up_to(n);

    std::vector<int> exps;
    exps.reserve(ps.size());
    for (int p : ps) exps.push_back(exp_in_fact(n, p));

    bool all_even = true;
    for (int e : exps)
        if (e & 1) {
            all_even = false;
            break;
        }

    const std::array<int, 4> bigP{2, 3, 5, 7};

    for (size_t i = 0; i < ps.size(); ++i) {
        if (ps[i] > 7) assert(exps[i] <= 9);
    }

    int max_val = 1;
    for (size_t i = 0; i < ps.size(); ++i) {
        max_val = std::max(max_val, exps[i] + 1);
    }

    const auto basis = primes_up_to(max_val);
    const int P = (int)basis.size();
    std::unordered_map<int, int> p2idx;
    p2idx.reserve(basis.size());
    for (int i = 0; i < P; ++i) p2idx[basis[i]] = i;

    std::vector<std::vector<std::int8_t>> fac((size_t)max_val + 1, std::vector<std::int8_t>((size_t)P, 0));
    for (int x = 2; x <= max_val; ++x) {
        int m = x;
        for (int i = 0; i < P; ++i) {
            const int p = basis[i];
            if (p * p > m) break;
            while (m % p == 0) {
                ++fac[(size_t)x][(size_t)i];
                m /= p;
            }
        }
        if (m > 1) {
            const int i = p2idx[m];
            ++fac[(size_t)x][(size_t)i];
        }
    }

    const int idx2 = p2idx[2];
    const int idx3 = p2idx[3];
    const int idx5 = p2idx[5];
    const int idx7 = p2idx[7];

    std::vector<int> big_indices;
    for (int i = 0; i < P; ++i) {
        if (basis[i] > 7) big_indices.push_back(i);
    }

    std::vector<std::vector<std::vector<std::int8_t>>> diffs;
    diffs.reserve(4);
    for (int bp : bigP) {
        if (bp > n) {
            diffs.push_back({});
            continue;
        }
        int e = exp_in_fact(n, bp);
        const int E = e + 2;
        std::vector<std::vector<std::int8_t>> v;
        v.reserve((size_t)(E - 1));
        for (int y = 1; y <= e + 1; ++y) {
            const int z = E - y;
            std::vector<std::int8_t> dv((size_t)P, 0);
            for (int i = 0; i < P; ++i) dv[(size_t)i] = (std::int8_t)(fac[(size_t)y][(size_t)i] - fac[(size_t)z][(size_t)i]);
            v.push_back(std::move(dv));
        }
        diffs.push_back(std::move(v));
    }

    std::vector<std::vector<std::int8_t>> pair01;
    std::vector<std::vector<std::int8_t>> pair23;

    pair01.reserve(diffs[0].size() * diffs[1].size());
    for (const auto& a : diffs[0]) {
        for (const auto& b : diffs[1]) {
            std::vector<std::int8_t> s((size_t)P, 0);
            for (int i = 0; i < P; ++i) s[(size_t)i] = (std::int8_t)(a[(size_t)i] + b[(size_t)i]);
            pair01.push_back(std::move(s));
        }
    }

    pair23.reserve(diffs[2].size() * diffs[3].size());
    for (const auto& a : diffs[2]) {
        for (const auto& b : diffs[3]) {
            std::vector<std::int8_t> s((size_t)P, 0);
            for (int i = 0; i < P; ++i) s[(size_t)i] = (std::int8_t)(a[(size_t)i] + b[(size_t)i]);
            pair23.push_back(std::move(s));
        }
    }

    std::unordered_map<u64, u64> map1;
    map1.reserve(10000);

    for (const auto& a : pair01) {
        for (const auto& b : pair23) {
            bool ok = true;
            for (int i : big_indices) {
                if ((int)a[(size_t)i] + (int)b[(size_t)i] != 0) {
                    ok = false;
                    break;
                }
            }
            if (!ok) continue;

            const int d2 = (int)a[(size_t)idx2] + (int)b[(size_t)idx2];
            const int d3 = (int)a[(size_t)idx3] + (int)b[(size_t)idx3];
            const int d5 = (int)a[(size_t)idx5] + (int)b[(size_t)idx5];
            const int d7 = (int)a[(size_t)idx7] + (int)b[(size_t)idx7];
            ++map1[pack4(d2, d3, d5, d7)];
        }
    }

    std::unordered_map<u64, u64> map2;
    map2.reserve(30000);
    map2[pack4(0, 0, 0, 0)] = 1;

    for (size_t i = 0; i < ps.size(); ++i) {
        const int p = ps[i];
        if (p <= 7) continue;
        const int e = exps[i];
        const int E = e + 2;

        std::vector<std::array<int, 4>> deltas;
        deltas.reserve((size_t)(E - 1));
        for (int y = 1; y <= e + 1; ++y) {
            const int z = E - y;
            deltas.push_back({(int)fac[(size_t)y][(size_t)idx2] - (int)fac[(size_t)z][(size_t)idx2],
                              (int)fac[(size_t)y][(size_t)idx3] - (int)fac[(size_t)z][(size_t)idx3],
                              (int)fac[(size_t)y][(size_t)idx5] - (int)fac[(size_t)z][(size_t)idx5],
                              (int)fac[(size_t)y][(size_t)idx7] - (int)fac[(size_t)z][(size_t)idx7]});
        }

        std::unordered_map<u64, u64> nxt;
        nxt.reserve(map2.size() * deltas.size());
        for (const auto& kv : map2) {
            int d2, d3, d5, d7;
            unpack4(kv.first, d2, d3, d5, d7);
            const u64 cnt = kv.second;
            for (const auto& dd : deltas) {
                const u64 nk = pack4(d2 + dd[0], d3 + dd[1], d5 + dd[2], d7 + dd[3]);
                nxt[nk] += cnt;
            }
        }
        map2.swap(nxt);
    }

    u128 M = 0;
    for (const auto& kv : map1) {
        int d2, d3, d5, d7;
        unpack4(kv.first, d2, d3, d5, d7);
        const u64 need = pack4(-d2, -d3, -d5, -d7);
        auto it = map2.find(need);
        if (it == map2.end()) continue;
        M += (u128)kv.second * (u128)it->second;
    }

    const u128 fixed = all_even ? 1 : 0;
    return (M + fixed) / 2;
}

static void print_u128(u128 x) {
    if (x == 0) {
        std::cout << '0';
        return;
    }
    char buf[64];
    int n = 0;
    while (x > 0) {
        buf[n++] = (char)('0' + (u64)(x % 10));
        x /= 10;
    }
    while (n--) std::cout << buf[n];
}

int main() {
    if (C_factorial(10) != 3) {
        std::cerr << "Validation failed: C(10!)\n";
        return 1;
    }

    const u128 ans = C_factorial(100);
    print_u128(ans);
    std::cout << "\n";
    return 0;
}

Python

import sys
import math

def primes_up_to(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(math.sqrt(n)) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False
    return [i for i in range(2, n + 1) if is_prime[i]]

def exp_in_fact(n, p):
    e = 0
    while n > 0:
        n //= p
        e += n
    return e

def solve_factorial(n):
    ps = primes_up_to(n)
    exps = [exp_in_fact(n, p) for p in ps]
    
    all_even = all(e % 2 == 0 for e in exps)
    
    bigP = [2, 3, 5, 7]
    
    max_val = max(exps) + 1 if exps else 1
    basis = primes_up_to(max_val)
    P = len(basis)
    p2idx = {p: i for i, p in enumerate(basis)}
    
    fac = [[0] * P for _ in range(max_val + 1)]
    for x in range(2, max_val + 1):
        m = x
        for i, p in enumerate(basis):
            if p * p > m: break
            while m % p == 0:
                fac[x][i] += 1
                m //= p
        if m > 1:
            fac[x][p2idx[m]] += 1
            
    idx2 = p2idx.get(2, -1)
    idx3 = p2idx.get(3, -1)
    idx5 = p2idx.get(5, -1)
    idx7 = p2idx.get(7, -1)
    
    big_indices = [i for i, p in enumerate(basis) if p > 7]
    
    diffs = []
    for bp in bigP:
        if bp > n:
            diffs.append([])
            continue
            
        e = exp_in_fact(n, bp)
        E = e + 2
        v = []
        for y in range(1, e + 2):
            z = E - y
            dv = [fac[y][i] - fac[z][i] for i in range(P)]
            v.append(dv)
        diffs.append(v)
        
    while len(diffs) < 4:
        diffs.append([])
        
    pair01 = []
    for a in diffs[0]:
        for b in diffs[1]:
            s = [a[i] + b[i] for i in range(P)]
            pair01.append(s)
            
    pair23 = []
    for a in diffs[2]:
        for b in diffs[3]:
            s = [a[i] + b[i] for i in range(P)]
            pair23.append(s)
            
    map1 = {}
    for a in pair01:
        for b in pair23:
            ok = True
            for i in big_indices:
                if a[i] + b[i] != 0:
                    ok = False
                    break
            if not ok: continue
            
            d2 = a[idx2] + b[idx2]
            d3 = a[idx3] + b[idx3]
            d5 = a[idx5] + b[idx5]
            d7 = a[idx7] + b[idx7]
            key = (d2, d3, d5, d7)
            map1[key] = map1.get(key, 0) + 1
            
    map2 = {}
    map2[(0, 0, 0, 0)] = 1
    
    for i in range(len(ps)):
        p = ps[i]
        if p <= 7: continue
        e = exps[i]
        E = e + 2
        
        deltas = []
        for y in range(1, e + 2):
            z = E - y
            deltas.append((
                fac[y][idx2] - fac[z][idx2],
                fac[y][idx3] - fac[z][idx3],
                fac[y][idx5] - fac[z][idx5],
                fac[y][idx7] - fac[z][idx7]
            ))
            
        nxt = {}
        for kv, cnt in map2.items():
            d2, d3, d5, d7 = kv
            for dd in deltas:
                nk = (d2 + dd[0], d3 + dd[1], d5 + dd[2], d7 + dd[3])
                nxt[nk] = nxt.get(nk, 0) + cnt
        map2 = nxt
        
    M = 0
    for kv, cnt in map1.items():
        d2, d3, d5, d7 = kv
        need = (-d2, -d3, -d5, -d7)
        if need in map2:
            M += cnt * map2[need]
            
    fixed = 1 if all_even else 0
    ans = (M + fixed) // 2
    return ans
    
def solve():
    return str(solve_factorial(100))

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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler598 {

    static List<Integer> primesUpTo(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int i = 2; i * i <= n; i++) {
            if (!isPrime[i])
                continue;
            for (int j = i * i; j <= n; j += i)
                isPrime[j] = false;
        }
        List<Integer> ps = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                ps.add(i);
        }
        return ps;
    }

    static int expInFact(int n, int p) {
        int e = 0;
        while (n > 0) {
            n /= p;
            e += n;
        }
        return e;
    }

    static long pack4(int d2, int d3, int d5, int d7) {
        int OFF = 128;
        long a = d2 + OFF;
        long b = d3 + OFF;
        long c = d5 + OFF;
        long d = d7 + OFF;
        return (a & 0xFFFF) | ((b & 0xFFFF) << 16) | ((c & 0xFFFF) << 32) | ((d & 0xFFFF) << 48);
    }

    static int[] unpack4(long key) {
        int OFF = 128;
        int d2 = (int) (key & 0xFFFF) - OFF;
        int d3 = (int) ((key >> 16) & 0xFFFF) - OFF;
        int d5 = (int) ((key >> 32) & 0xFFFF) - OFF;
        int d7 = (int) ((key >> 48) & 0xFFFF) - OFF;
        return new int[] { d2, d3, d5, d7 };
    }

    static String solveFactorial(int n) {
        List<Integer> ps = primesUpTo(n);
        int[] exps = new int[ps.size()];
        for (int i = 0; i < ps.size(); i++) {
            exps[i] = expInFact(n, ps.get(i));
        }

        boolean allEven = true;
        for (int e : exps) {
            if (e % 2 != 0) {
                allEven = false;
                break;
            }
        }

        int[] bigP = { 2, 3, 5, 7 };
        int maxVal = 1;
        for (int e : exps) {
            maxVal = Math.max(maxVal, e + 1);
        }

        List<Integer> basis = primesUpTo(maxVal);
        int P = basis.size();
        Map<Integer, Integer> p2idx = new HashMap<>();
        for (int i = 0; i < P; i++)
            p2idx.put(basis.get(i), i);

        byte[][] fac = new byte[maxVal + 1][P];
        for (int x = 2; x <= maxVal; x++) {
            int m = x;
            for (int i = 0; i < P; i++) {
                int p = basis.get(i);
                if (p * p > m)
                    break;
                while (m % p == 0) {
                    fac[x][i]++;
                    m /= p;
                }
            }
            if (m > 1) {
                fac[x][p2idx.get(m)]++;
            }
        }

        int idx2 = p2idx.getOrDefault(2, -1);
        int idx3 = p2idx.getOrDefault(3, -1);
        int idx5 = p2idx.getOrDefault(5, -1);
        int idx7 = p2idx.getOrDefault(7, -1);

        List<Integer> bigIndices = new ArrayList<>();
        for (int i = 0; i < P; i++) {
            if (basis.get(i) > 7)
                bigIndices.add(i);
        }

        List<List<byte[]>> diffs = new ArrayList<>();
        for (int bp : bigP) {
            if (bp > n) {
                diffs.add(new ArrayList<>());
                continue;
            }
            int e = expInFact(n, bp);
            int E = e + 2;
            List<byte[]> v = new ArrayList<>();
            for (int y = 1; y <= e + 1; y++) {
                int z = E - y;
                byte[] dv = new byte[P];
                for (int i = 0; i < P; i++) {
                    dv[i] = (byte) (fac[y][i] - fac[z][i]);
                }
                v.add(dv);
            }
            diffs.add(v);
        }

        List<byte[]> pair01 = new ArrayList<>();
        for (byte[] a : diffs.get(0)) {
            for (byte[] b : diffs.get(1)) {
                byte[] s = new byte[P];
                for (int i = 0; i < P; i++)
                    s[i] = (byte) (a[i] + b[i]);
                pair01.add(s);
            }
        }

        List<byte[]> pair23 = new ArrayList<>();
        for (byte[] a : diffs.get(2)) {
            for (byte[] b : diffs.get(3)) {
                byte[] s = new byte[P];
                for (int i = 0; i < P; i++)
                    s[i] = (byte) (a[i] + b[i]);
                pair23.add(s);
            }
        }

        Map<Long, Long> map1 = new HashMap<>();
        for (byte[] a : pair01) {
            for (byte[] b : pair23) {
                boolean ok = true;
                for (int i : bigIndices) {
                    if (a[i] + b[i] != 0) {
                        ok = false;
                        break;
                    }
                }
                if (!ok)
                    continue;

                int d2 = a[idx2] + b[idx2];
                int d3 = a[idx3] + b[idx3];
                int d5 = a[idx5] + b[idx5];
                int d7 = a[idx7] + b[idx7];
                long key = pack4(d2, d3, d5, d7);
                map1.put(key, map1.getOrDefault(key, 0L) + 1L);
            }
        }

        Map<Long, Long> map2 = new HashMap<>();
        map2.put(pack4(0, 0, 0, 0), 1L);

        for (int i = 0; i < ps.size(); i++) {
            int p = ps.get(i);
            if (p <= 7)
                continue;
            int e = exps[i];
            int E = e + 2;

            List<int[]> deltas = new ArrayList<>();
            for (int y = 1; y <= e + 1; y++) {
                int z = E - y;
                deltas.add(new int[] {
                        fac[y][idx2] - fac[z][idx2],
                        fac[y][idx3] - fac[z][idx3],
                        fac[y][idx5] - fac[z][idx5],
                        fac[y][idx7] - fac[z][idx7]
                });
            }

            Map<Long, Long> nxt = new HashMap<>();
            for (Map.Entry<Long, Long> kv : map2.entrySet()) {
                long key = kv.getKey();
                long cnt = kv.getValue();
                int[] ds = unpack4(key);
                for (int[] dd : deltas) {
                    long nk = pack4(ds[0] + dd[0], ds[1] + dd[1], ds[2] + dd[2], ds[3] + dd[3]);
                    nxt.put(nk, nxt.getOrDefault(nk, 0L) + cnt);
                }
            }
            map2 = nxt;
        }

        java.math.BigInteger M = java.math.BigInteger.ZERO;
        for (Map.Entry<Long, Long> kv : map1.entrySet()) {
            int[] ds = unpack4(kv.getKey());
            long need = pack4(-ds[0], -ds[1], -ds[2], -ds[3]);
            Long cnt2 = map2.get(need);
            if (cnt2 != null) {
                java.math.BigInteger terms = java.math.BigInteger.valueOf(kv.getValue())
                        .multiply(java.math.BigInteger.valueOf(cnt2));
                M = M.add(terms);
            }
        }

        java.math.BigInteger fixed = allEven ? java.math.BigInteger.ONE : java.math.BigInteger.ZERO;
        return M.add(fixed).divide(java.math.BigInteger.valueOf(2)).toString();
    }

    public static String solve() {
        return solveFactorial(100);
    }

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