Problem 548: Gozinta Chains

Project Euler Problem 548 Solution

EulerSolve provides an optimized solution for Project Euler Problem 548, Gozinta Chains, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A gozinta chain for $n$ is a sequence $$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$ in which every term divides the next one. Let $g(n)$ be the number of such chains. The task is to sum all integers $n\le 10^{16}$ satisfying $$g(n)=n.$$ The decisive observation is that $g(n)$ depends only on the sorted prime-exponent pattern of $n$, not on the actual primes. That turns a hopeless search over all integers up to $10^{16}$ into a manageable search over exponent patterns. Mathematical Approach Write the prime factorization in non-increasing exponent order: $$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$ Let $\mathbf a=(a_1,\dots,a_r)$ and $\omega=a_1+\cdots+a_r$. Step 1: Replace the Integer by Its Exponent Pattern Every divisor of $n$ has the form $$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$ So a gozinta chain is equivalent to a coordinatewise non-decreasing sequence of exponent vectors starting at $(0,\dots,0)$ and ending at $(a_1,\dots,a_r)$. Because that description uses only the exponents, the chain count is determined entirely by $\mathbf a$. We may therefore write $$g(n)=G(\mathbf a),$$ where $G$ is a function on exponent patterns. This is the key reduction used by the implementation....

Detailed mathematical approach

Problem Summary

A gozinta chain for $n$ is a sequence

$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n$$

in which every term divides the next one. Let $g(n)$ be the number of such chains. The task is to sum all integers $n\le 10^{16}$ satisfying

$$g(n)=n.$$

The decisive observation is that $g(n)$ depends only on the sorted prime-exponent pattern of $n$, not on the actual primes. That turns a hopeless search over all integers up to $10^{16}$ into a manageable search over exponent patterns.

Mathematical Approach

Write the prime factorization in non-increasing exponent order:

$$n=\prod_{i=1}^{r} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_r \gt 0.$$

Let $\mathbf a=(a_1,\dots,a_r)$ and $\omega=a_1+\cdots+a_r$.

Step 1: Replace the Integer by Its Exponent Pattern

Every divisor of $n$ has the form

$$\prod_{i=1}^{r} p_i^{b_i},\qquad 0\le b_i\le a_i.$$

So a gozinta chain is equivalent to a coordinatewise non-decreasing sequence of exponent vectors starting at $(0,\dots,0)$ and ending at $(a_1,\dots,a_r)$.

Because that description uses only the exponents, the chain count is determined entirely by $\mathbf a$. We may therefore write

$$g(n)=G(\mathbf a),$$

where $G$ is a function on exponent patterns. This is the key reduction used by the implementation.

Step 2: Count Chains with a Fixed Number of Steps

Suppose the chain has exactly $k$ nontrivial divisibility steps:

$$1=d_0 \lt d_1 \lt \cdots \lt d_k=n.$$

Define the step ratios $q_t=d_t/d_{t-1}$. Then each $q_t \gt 1$ and

$$q_1q_2\cdots q_k=n.$$

For one prime $p_i^{a_i}$, distribute the exponent $a_i$ across those $k$ ratios:

$$e_{i,1}+e_{i,2}+\cdots+e_{i,k}=a_i,\qquad e_{i,t}\ge 0.$$

The number of weak compositions of $a_i$ into $k$ parts is

$$\binom{a_i+k-1}{k-1}.$$

Independence across primes gives

$$T_k(\mathbf a)=\prod_{i=1}^{r}\binom{a_i+k-1}{k-1}.$$

This counts ordered $k$-step decompositions, but it still allows some ratios to equal $1$, so it overcounts genuine chains.

Step 3: Use Inclusion-Exclusion to Remove Empty Steps

A genuine gozinta chain has no empty step. If exactly $j$ positions among $k$ are active, the exponent assignments are counted by $T_j(\mathbf a)$, and there are

$$\binom{k}{j}$$

ways to choose those active positions. Inclusion-exclusion over the empty positions yields the number of chains with exactly $k$ nontrivial steps:

$$F_k(\mathbf a)=\sum_{j=1}^{k}(-1)^{k-j}\binom{k}{j}T_j(\mathbf a).$$

This alternating sum is the core combinatorial formula in the C++, Python, and Java implementations.

Step 4: Sum over All Admissible Lengths

Each nontrivial step increases the total exponent sum by at least $1$, so the number of steps cannot exceed

$$\omega=a_1+\cdots+a_r.$$

Therefore the total number of gozinta chains for pattern $\mathbf a$ is

$$G(\mathbf a)=\sum_{k=1}^{\omega}F_k(\mathbf a).$$

At this point the original divisor-chain problem has been reduced to exact finite sums of binomial coefficients.

Step 5: Convert $g(n)=n$ into a Pattern Fixed-Point Test

Let

$$m=G(\mathbf a).$$

If the exponent pattern of $m$ is again $\mathbf a$, then $m$ itself has pattern $\mathbf a$, so

$$g(m)=G(\mathbf a)=m.$$

Conversely, every fixed point $n$ produces its own exponent pattern and passes this test. So the search does not need to enumerate all integers; it only needs to enumerate exponent patterns, compute $G(\mathbf a)$, factor the resulting integer, and compare patterns.

Step 6: Prune the Search with the Smallest Representative

Among all integers having exponent pattern $\mathbf a$, the smallest one is obtained by matching the largest exponent with the smallest prime, the next largest exponent with the next prime, and so on:

$$n_{\min}(\mathbf a)=2^{a_1}3^{a_2}5^{a_3}\cdots.$$

$$n_{\min}(\mathbf a)\gt 10^{16},$$

then no integer with that pattern can contribute and the entire branch is discarded. There is a second useful rejection test:

$$G(\mathbf a)\lt n_{\min}(\mathbf a)\implies \text{pattern}(G(\mathbf a))\ne \mathbf a.$$

So factorization is only attempted when the candidate is large enough even to have the required pattern.

Worked Example

Take $\mathbf a=(4,1)$, the exponent pattern of $48=2^4\cdot 3$. Here $\omega=5$. First compute

$$\begin{aligned} T_1&=\binom{4}{0}\binom{1}{0}=1,\\ T_2&=\binom{5}{1}\binom{2}{1}=10,\\ T_3&=\binom{6}{2}\binom{3}{2}=45,\\ T_4&=\binom{7}{3}\binom{4}{3}=140,\\ T_5&=\binom{8}{4}\binom{5}{4}=350. \end{aligned}$$

Then inclusion-exclusion gives

$$\begin{aligned} F_1&=1,\\ F_2&=-2T_1+T_2=8,\\ F_3&=3T_1-3T_2+T_3=18,\\ F_4&=-4T_1+6T_2-4T_3+T_4=16,\\ F_5&=5T_1-10T_2+10T_3-5T_4+T_5=5. \end{aligned}$$

$$G(4,1)=1+8+18+16+5=48.$$

Since $48$ has exponent pattern $(4,1)$, it is a true fixed point and must be included in the final sum. By contrast, $(2,1)$ gives $G(2,1)=8$, whose pattern is $(3)$, so $12$ is not a fixed point.

How the Code Works

The C++, Python, and Java implementations recursively enumerate exponent patterns in non-increasing order. The search starts with a safe upper bound on the first exponent and extends the pattern one prime at a time. Every recursive step updates the smallest possible integer with that pattern, so branches exceeding $10^{16}$ are cut immediately.

For each surviving pattern, the implementation builds the needed binomial coefficients exactly, evaluates the product terms $T_k$, and applies the inclusion-exclusion formula to obtain $F_k$ and then $G(\mathbf a)$. The running total is stopped as soon as it exceeds the global limit, because larger values can never satisfy the problem condition.

Only candidates $m$ with $m\ge n_{\min}(\mathbf a)$ are factored. Their prime exponents are extracted, sorted in non-increasing order, and compared with the original pattern. Previously analyzed candidates are cached, valid fixed points are deduplicated, and the final sum is accumulated with wide integer arithmetic.

Complexity Analysis

Let $P$ be the number of exponent patterns whose smallest representative does not exceed $10^{16}$. For a pattern $\mathbf a$ with weight $\omega=\sum a_i$ and length $r$, evaluating $G(\mathbf a)$ costs $O(\omega^2+r\omega)$ arithmetic operations: $O(\omega^2)$ for the binomial table and the inclusion-exclusion layers, and $O(r\omega)$ for the products $T_k$.

The recursive enumeration is practical because the lower-bound pruning keeps $P$ small. The factorization stage uses deterministic Miller-Rabin primality testing for 64-bit integers together with Pollard-Rho splitting. That stage has no simple closed worst-case formula, but below $10^{16}$ it is fast in practice and runs only on the comparatively small set of surviving candidates. Memory usage is modest: $O(\omega^2)$ while processing one pattern, plus the cache of analyzed candidates and the set of discovered fixed points.

Footnotes and References

Problem page: https://projecteuler.net/problem=548
Stars and bars / weak compositions: Wikipedia — Stars and bars
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
Miller-Rabin primality test: Wikipedia — Miller-Rabin primality test
Pollard-Rho factorization: Wikipedia — Pollard's rho algorithm

Problem 548 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;
using i128 = __int128_t;

static constexpr u64 kLimit = 10'000'000'000'000'000ULL;  // 1e16

static u64 gcd_u64(u64 a, u64 b) {
    while (b) {
        u64 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

static u64 mul_mod_u64(u64 a, u64 b, u64 mod) { return static_cast<u64>((static_cast<u128>(a) * b) % mod); }

static u64 pow_mod_u64(u64 a, u64 d, u64 mod) {
    u64 r = 1;
    while (d) {
        if (d & 1) r = mul_mod_u64(r, a, mod);
        a = mul_mod_u64(a, a, mod);
        d >>= 1;
    }
    return r;
}

static bool is_prime_u64(u64 n) {
    if (n < 2) return false;
    static u64 small_primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (u64 p : small_primes) {
        if (n == p) return true;
        if (n % p == 0) return false;
    }

    u64 d = n - 1, s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        ++s;
    }

    auto witness = [&](u64 a) -> bool {
        if (a % n == 0) return false;
        u64 x = pow_mod_u64(a, d, n);
        if (x == 1 || x == n - 1) return false;
        for (u64 i = 1; i < s; ++i) {
            x = mul_mod_u64(x, x, n);
            if (x == n - 1) return false;
        }
        return true;
    };

    static u64 bases[] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL};
    for (u64 a : bases) {
        if (witness(a)) return false;
    }
    return true;
}

static u64 splitmix64(u64& x) {
    u64 z = (x += 0x9e3779b97f4a7c15ULL);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
    return z ^ (z >> 31);
}

static thread_local u64 rng_state = 0x123456789abcdef0ULL;

static u64 rand_u64(u64 lo, u64 hi) {
    u64 r = splitmix64(rng_state);
    return lo + (hi > lo ? (r % (hi - lo + 1)) : 0);
}

static u64 pollard_rho(u64 n) {
    if ((n & 1ULL) == 0) return 2;
    if (n % 3ULL == 0) return 3;

    while (true) {
        u64 c = rand_u64(1, n - 1);
        u64 x = rand_u64(0, n - 1);
        u64 y = x;
        u64 d = 1;

        auto f = [&](u64 v) { return (mul_mod_u64(v, v, n) + c) % n; };

        while (d == 1) {
            x = f(x);
            y = f(f(y));
            u64 diff = (x > y) ? (x - y) : (y - x);
            d = gcd_u64(diff, n);
        }
        if (d != n) return d;
    }
}

static void factor_rec(u64 n, std::vector<u64>& out) {
    if (n == 1) return;
    if (is_prime_u64(n)) {
        out.push_back(n);
        return;
    }
    u64 d = pollard_rho(n);
    factor_rec(d, out);
    factor_rec(n / d, out);
}

static std::vector<int> exponent_pattern(u64 n) {
    if (n == 1) return {};
    std::vector<u64> fac;
    factor_rec(n, fac);
    std::sort(fac.begin(), fac.end());
    std::vector<int> exps;
    for (std::size_t i = 0; i < fac.size();) {
        std::size_t j = i;
        while (j < fac.size() && fac[j] == fac[i]) ++j;
        exps.push_back(static_cast<int>(j - i));
        i = j;
    }
    std::sort(exps.begin(), exps.end(), std::greater<int>());
    return exps;
}

static u128 binom_u128(int n, int k) {
    if (k < 0 || k > n) return 0;
    k = std::min(k, n - k);
    u128 r = 1;
    for (int i = 1; i <= k; ++i) {
        r = (r * static_cast<u128>(n - k + i)) / static_cast<u128>(i);
    }
    return r;
}

static u64 g_from_pattern(const std::vector<int>& exps, u64 limit) {
    if (exps.empty()) return 1;
    int omega = 0;
    for (int e : exps) omega += e;

    std::vector<std::vector<u128>> C(omega + 1, std::vector<u128>(omega + 1, 0));
    C[0][0] = 1;
    for (int n = 1; n <= omega; ++n) {
        C[n][0] = C[n][n] = 1;
        for (int k = 1; k < n; ++k) C[n][k] = C[n - 1][k - 1] + C[n - 1][k];
    }

    std::vector<u128> total(static_cast<std::size_t>(omega + 1), 0);
    for (int k = 1; k <= omega; ++k) {
        u128 prod = 1;
        for (int a : exps) {
            prod *= binom_u128(a + k - 1, k - 1);
        }
        total[static_cast<std::size_t>(k)] = prod;
    }

    u128 g = 0;
    for (int k = 1; k <= omega; ++k) {
        i128 fk = 0;
        for (int i = 1; i <= k; ++i) {
            const u128 term = C[k][i] * total[static_cast<std::size_t>(i)];
            if (((k - i) & 1) != 0) fk -= static_cast<i128>(term);
            else fk += static_cast<i128>(term);
        }
        g += static_cast<u128>(fk);
        if (g > static_cast<u128>(limit)) return limit + 1;
    }
    return static_cast<u64>(g);
}

static std::string to_string_u128(u128 x) {
    if (x == 0) return "0";
    std::string s;
    while (x > 0) {
        int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

struct Enumerator {
    std::vector<u64> primes;
    std::vector<int> cur;
    std::unordered_map<u64, std::vector<int>> pattern_cache;
    std::vector<u64> solutions;

    Enumerator() {
        // Enough primes for any exponent pattern with minimal value <= 1e16.
        primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71};
        pattern_cache.reserve(1 << 16);
        solutions.reserve(64);
    }

    const std::vector<int>& cached_pattern(u64 n) {
        auto it = pattern_cache.find(n);
        if (it != pattern_cache.end()) return it->second;
        std::vector<int> p = exponent_pattern(n);
        auto [ins_it, _] = pattern_cache.emplace(n, std::move(p));
        return ins_it->second;
    }

    void dfs(int last_exp, u128 min_n) {
        const u64 g = g_from_pattern(cur, kLimit);
        if (g > kLimit) return;

        if (static_cast<u128>(g) >= min_n) {
            const std::vector<int>& pat = cached_pattern(g);
            if (pat == cur) solutions.push_back(g);
        }

        const std::size_t idx = cur.size();
        if (idx >= primes.size()) return;
        const u64 p = primes[idx];

        u128 pow = 1;
        for (int e = 1; e <= last_exp; ++e) {
            pow *= static_cast<u128>(p);
            const u128 next_min = min_n * pow;
            if (next_min > kLimit) break;
            cur.push_back(e);
            dfs(e, next_min);
            cur.pop_back();
        }
    }

    std::vector<u64> run() {
        cur.clear();
        dfs(53, 1);
        std::sort(solutions.begin(), solutions.end());
        solutions.erase(std::unique(solutions.begin(), solutions.end()), solutions.end());
        return solutions;
    }
};

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    // Statement checkpoints.
    assert(g_from_pattern(exponent_pattern(12), kLimit) == 8ULL);
    assert(g_from_pattern(exponent_pattern(48), kLimit) == 48ULL);
    assert(g_from_pattern(exponent_pattern(120), kLimit) == 132ULL);

    Enumerator en;
    const std::vector<u64> sols = en.run();

    u128 sum = 0;
    for (u64 x : sols) sum += static_cast<u128>(x);
    std::cout << to_string_u128(sum) << '\n';
    return 0;
}

Python

import math

kLimit = 10**16

def gcd_u64(a, b):
    return math.gcd(a, b)

def is_prime_u64(n):
    if n < 2: return False
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    for p in small_primes:
        if n == p: return True
        if n % p == 0: return False
        
    d = n - 1
    s = 0
    while (d & 1) == 0:
        d >>= 1
        s += 1
        
    def witness(a):
        if a % n == 0: return False
        x = pow(a, d, n)
        if x == 1 or x == n - 1: return False
        for _ in range(1, s):
            x = (x * x) % n
            if x == n - 1: return False
        return True
        
    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
    for a in bases:
        if witness(a): return False
    return True

rng_state = 0x123456789abcdef0

def splitmix64():
    global rng_state
    rng_state = (rng_state + 0x9e3779b97f4a7c15) & 0xffffffffffffffff
    z = rng_state
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9 & 0xffffffffffffffff
    z = (z ^ (z >> 27)) * 0x94d049bb133111eb & 0xffffffffffffffff
    return z ^ (z >> 31)

def rand_u64(lo, hi):
    r = splitmix64()
    if hi > lo:
        return lo + (r % (hi - lo + 1))
    return lo

def pollard_rho(n):
    if (n & 1) == 0: return 2
    if n % 3 == 0: return 3
    
    while True:
        c = rand_u64(1, n - 1)
        x = rand_u64(0, n - 1)
        y = x
        d = 1
        
        def f(v):
            return (v * v + c) % n
            
        while d == 1:
            x = f(x)
            y = f(f(y))
            diff = x - y if x > y else y - x
            d = gcd_u64(diff, n)
            
        if d != n: return d

def factor_rec(n, out):
    if n == 1: return
    if is_prime_u64(n):
        out.append(n)
        return
    d = pollard_rho(n)
    factor_rec(d, out)
    factor_rec(n // d, out)

def exponent_pattern(n):
    if n == 1: return ()
    fac = []
    factor_rec(n, fac)
    fac.sort()
    
    exps = []
    i = 0
    while i < len(fac):
        j = i
        while j < len(fac) and fac[j] == fac[i]:
            j += 1
        exps.append(j - i)
        i = j
    exps.sort(reverse=True)
    return tuple(exps)

def binom(n, k):
    if k < 0 or k > n: return 0
    k = min(k, n - k)
    r = 1
    for i in range(1, k + 1):
        r = r * (n - k + i) // i
    return r

def g_from_pattern(exps, limit):
    if not exps: return 1
    omega = sum(exps)
    
    C = [[0] * (omega + 1) for _ in range(omega + 1)]
    C[0][0] = 1
    for n in range(1, omega + 1):
        C[n][0] = C[n][n] = 1
        for k in range(1, n):
            C[n][k] = C[n - 1][k - 1] + C[n - 1][k]
            
    total = [0] * (omega + 1)
    for k in range(1, omega + 1):
        prod = 1
        for a in exps:
            prod *= binom(a + k - 1, k - 1)
        total[k] = prod
        
    g = 0
    for k in range(1, omega + 1):
        fk = 0
        for i in range(1, k + 1):
            term = C[k][i] * total[i]
            if ((k - i) & 1) != 0:
                fk -= term
            else:
                fk += term
        g += fk
        if g > limit: return limit + 1
    return g

class Enumerator:
    def __init__(self):
        self.primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]
        self.pattern_cache = {}
        self.solutions = set()
        self.cur = []
        
    def cached_pattern(self, n):
        if n in self.pattern_cache:
            return self.pattern_cache[n]
        p = exponent_pattern(n)
        self.pattern_cache[n] = p
        return p
        
    def dfs(self, last_exp, min_n):
        t_cur = tuple(self.cur)
        g = g_from_pattern(t_cur, kLimit)
        if g > kLimit: return
        
        if g >= min_n:
            pat = self.cached_pattern(g)
            if pat == t_cur:
                self.solutions.add(g)
                
        idx = len(self.cur)
        if idx >= len(self.primes): return
        p = self.primes[idx]
        
        pow_val = 1
        for e in range(1, last_exp + 1):
            pow_val *= p
            next_min = min_n * pow_val
            if next_min > kLimit: break
            self.cur.append(e)
            self.dfs(e, next_min)
            self.cur.pop()
            
    def run(self):
        self.cur = []
        self.dfs(53, 1)
        return sorted(list(self.solutions))

def solve():
    en = Enumerator()
    sols = en.run()
    return str(sum(sols))

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;

public class Euler548 {
    static final long kLimit = 10000000000000000L; // 1e16

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

    static long mulModU64(long a, long b, long mod) {
        long res = 0;
        a %= mod;
        while (b > 0) {
            if ((b & 1) != 0) {
                res = (res + a) % mod;
            }
            a = (a << 1) % mod;
            b >>= 1;
        }
        return res;
    }

    static long powModU64(long a, long d, long mod) {
        long r = 1;
        while (d > 0) {
            if ((d & 1) != 0)
                r = mulModU64(r, a, mod);
            a = mulModU64(a, a, mod);
            d >>= 1;
        }
        return r;
    }

    static boolean isPrimeU64(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }

        long d = n - 1;
        long s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }

        long[] bases = { 2L, 325L, 9375L, 28178L, 450775L, 9780504L, 1795265022L };
        for (long a : bases) {
            if (a % n == 0)
                continue;
            long x = powModU64(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean possiblePrime = false;
            for (long i = 1; i < s; i++) {
                x = mulModU64(x, x, n);
                if (x == n - 1) {
                    possiblePrime = true;
                    break;
                }
            }
            if (!possiblePrime)
                return false;
        }
        return true;
    }

    static long rngState = 0x123456789abcdef0L;

    static long splitmix64() {
        rngState += 0x9e3779b97f4a7c15L;
        long z = rngState;
        z = (z ^ (z >>> 30)) * 0xbf58476d1ce4e5b9L;
        z = (z ^ (z >>> 27)) * 0x94d049bb133111ebL;
        return z ^ (z >>> 31);
    }

    static long randU64(long lo, long hi) {
        long r = splitmix64();
        if (hi > lo) {
            return lo + Long.remainderUnsigned(r, hi - lo + 1);
        }
        return lo;
    }

    static long pollardRho(long n) {
        if ((n & 1) == 0)
            return 2;
        if (n % 3 == 0)
            return 3;

        while (true) {
            long c = randU64(1, n - 1);
            long x = randU64(0, n - 1);
            long y = x;
            long d = 1;

            while (d == 1) {
                x = (mulModU64(x, x, n) + c) % n;
                y = (mulModU64(y, y, n) + c) % n;
                y = (mulModU64(y, y, n) + c) % n;
                long diff = (x > y) ? (x - y) : (y - x);
                d = gcdU64(diff, n);
            }
            if (d != n)
                return d;
        }
    }

    static void factorRec(long n, List<Long> out) {
        if (n == 1)
            return;
        if (isPrimeU64(n)) {
            out.add(n);
            return;
        }
        long d = pollardRho(n);
        factorRec(d, out);
        factorRec(n / d, out);
    }

    static List<Integer> exponentPattern(long n) {
        if (n == 1)
            return Collections.emptyList();
        List<Long> fac = new ArrayList<>();
        factorRec(n, fac);
        Collections.sort(fac);
        List<Integer> exps = new ArrayList<>();
        for (int i = 0; i < fac.size();) {
            int j = i;
            while (j < fac.size() && fac.get(j).equals(fac.get(i)))
                j++;
            exps.add(j - i);
            i = j;
        }
        exps.sort(Collections.reverseOrder());
        return exps;
    }

    static java.math.BigInteger binomU128(int n, int k) {
        if (k < 0 || k > n)
            return java.math.BigInteger.ZERO;
        k = Math.min(k, n - k);
        java.math.BigInteger r = java.math.BigInteger.ONE;
        for (int i = 1; i <= k; i++) {
            r = r.multiply(java.math.BigInteger.valueOf(n - k + i))
                    .divide(java.math.BigInteger.valueOf(i));
        }
        return r;
    }

    static long gFromPattern(List<Integer> exps, long limit) {
        if (exps.isEmpty())
            return 1;
        int omega = 0;
        for (int e : exps)
            omega += e;

        java.math.BigInteger[][] C = new java.math.BigInteger[omega + 1][omega + 1];
        C[0][0] = java.math.BigInteger.ONE;
        for (int n = 1; n <= omega; n++) {
            C[n][0] = java.math.BigInteger.ONE;
            C[n][n] = java.math.BigInteger.ONE;
            for (int k = 1; k < n; k++) {
                C[n][k] = C[n - 1][k - 1].add(C[n - 1][k]);
            }
        }

        java.math.BigInteger[] total = new java.math.BigInteger[omega + 1];
        java.math.BigInteger limitBi = java.math.BigInteger.valueOf(limit);

        for (int k = 1; k <= omega; k++) {
            java.math.BigInteger prod = java.math.BigInteger.ONE;
            for (int a : exps) {
                prod = prod.multiply(binomU128(a + k - 1, k - 1));
            }
            total[k] = prod;
        }

        java.math.BigInteger g = java.math.BigInteger.ZERO;
        for (int k = 1; k <= omega; k++) {
            java.math.BigInteger fk = java.math.BigInteger.ZERO;
            for (int i = 1; i <= k; i++) {
                java.math.BigInteger term = C[k][i].multiply(total[i]);
                if (((k - i) & 1) != 0) {
                    fk = fk.subtract(term);
                } else {
                    fk = fk.add(term);
                }
            }
            g = g.add(fk);
            if (g.compareTo(limitBi) > 0)
                return limit + 1;
        }
        return g.longValueExact();
    }

    static class Enumerator {
        long[] primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71 };
        Map<Long, List<Integer>> patternCache = new HashMap<>();
        Set<Long> solutions = new HashSet<>();
        List<Integer> cur = new ArrayList<>();

        List<Integer> cachedPattern(long n) {
            if (patternCache.containsKey(n))
                return patternCache.get(n);
            List<Integer> p = exponentPattern(n);
            patternCache.put(n, p);
            return p;
        }

        void dfs(int lastExp, java.math.BigInteger minN) {
            long g = gFromPattern(cur, kLimit);
            if (g > kLimit)
                return;

            if (java.math.BigInteger.valueOf(g).compareTo(minN) >= 0) {
                List<Integer> pat = cachedPattern(g);
                if (pat.equals(cur))
                    solutions.add(g);
            }

            int idx = cur.size();
            if (idx >= primes.length)
                return;
            long p = primes[idx];

            java.math.BigInteger pow = java.math.BigInteger.ONE;
            java.math.BigInteger limitBi = java.math.BigInteger.valueOf(kLimit);
            for (int e = 1; e <= lastExp; e++) {
                pow = pow.multiply(java.math.BigInteger.valueOf(p));
                java.math.BigInteger nextMin = minN.multiply(pow);
                if (nextMin.compareTo(limitBi) > 0)
                    break;
                cur.add(e);
                dfs(e, nextMin);
                cur.remove(cur.size() - 1);
            }
        }

        List<Long> run() {
            cur.clear();
            dfs(53, java.math.BigInteger.ONE);
            List<Long> sorted = new ArrayList<>(solutions);
            Collections.sort(sorted);
            return sorted;
        }
    }

    public static String solve() {
        Enumerator en = new Enumerator();
        List<Long> sols = en.run();
        java.math.BigInteger sum = java.math.BigInteger.ZERO;
        for (long x : sols) {
            sum = sum.add(java.math.BigInteger.valueOf(x));
        }
        return sum.toString();
    }

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