Problem 861: Products of Bi-Unitary Divisors

Project Euler Problem 861 Solution

EulerSolve provides an optimized solution for Project Euler Problem 861, Products of Bi-Unitary Divisors, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer $n$, let $b(n)$ denote the number of bi-unitary divisors of $n$. The product of all bi-unitary divisors has the form $n^{b(n)/2}$, so the integers counted by $$Q_k(N)=\#\{1\le n\le N : b(n)=2k\}$$ are exactly the integers whose bi-unitary-divisor product equals $n^k$. The task is to compute $$\sum_{k=2}^{10} Q_k(10^{12}).$$ A brute-force scan up to $10^{12}$ is impossible, so the solution works directly with prime exponents and counts only the exponent patterns that can produce the small targets $2k\in\{4,6,\dots,20\}$. Mathematical Approach Write the prime factorization of $n$ as $$n=\prod_{i=1}^r p_i^{e_i}.$$ The key observation is that the bi-unitary-divisor count depends only on the exponents $e_i$, and the code exploits that structure very aggressively. Step 1: Count bi-unitary divisors of one prime power For $p^e$, the ordinary divisors are $1,p,p^2,\dots,p^e$. In the bi-unitary setting, every divisor survives except the middle divisor $p^{e/2}$ when $e$ is even....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $b(n)$ denote the number of bi-unitary divisors of $n$. The product of all bi-unitary divisors has the form $n^{b(n)/2}$, so the integers counted by

$$Q_k(N)=\#\{1\le n\le N : b(n)=2k\}$$

are exactly the integers whose bi-unitary-divisor product equals $n^k$. The task is to compute

$$\sum_{k=2}^{10} Q_k(10^{12}).$$

A brute-force scan up to $10^{12}$ is impossible, so the solution works directly with prime exponents and counts only the exponent patterns that can produce the small targets $2k\in\{4,6,\dots,20\}$.

Mathematical Approach

Write the prime factorization of $n$ as

$$n=\prod_{i=1}^r p_i^{e_i}.$$

The key observation is that the bi-unitary-divisor count depends only on the exponents $e_i$, and the code exploits that structure very aggressively.

Step 1: Count bi-unitary divisors of one prime power

For $p^e$, the ordinary divisors are $1,p,p^2,\dots,p^e$. In the bi-unitary setting, every divisor survives except the middle divisor $p^{e/2}$ when $e$ is even. Therefore

$$b(p^e)= \begin{cases} e, & e \text{ is even},\\ e+1, & e \text{ is odd}, \end{cases}$$

which is the same as

$$b(p^e)=e+(e\bmod 2).$$

Because prime powers contribute independently, $b(n)$ is multiplicative across the prime-power factors:

$$b(n)=\prod_{i=1}^r \bigl(e_i+(e_i\bmod 2)\bigr).$$

Step 2: Split $n$ into a heavy part and a squarefree part

Separate primes with exponent at least $2$ from the primes that appear only once:

$$n=\left(\prod_{i=1}^s p_i^{e_i}\right)\left(\prod_{j=1}^t q_j\right),\qquad e_i\ge 2,$$

where the $q_j$ are distinct primes different from every $p_i$. Then

$$b(n)=\left(\prod_{i=1}^s \bigl(e_i+(e_i\bmod 2)\bigr)\right)2^t.$$

Call the first factor $c$. Once the heavy exponents are fixed, the only freedom left is the number $t$ of exponent-$1$ primes. To force $b(n)=2k$, we need

$$2k=c\,2^t.$$

So $2k/c$ must be a power of two. If that quotient is not a power of two, the chosen heavy part cannot contribute to $Q_k(N)$.

Step 3: Enumerate only the heavy parts that can matter

Each heavy prime contributes at least a factor $p^2$ to $n$, so the numeric product grows quickly. At the same time, each heavy exponent contributes at least $2$ to $b(n)$, so the bi-unitary-divisor count also grows quickly.

Because the largest target here is $2k=20$, only very small values of $c$ are relevant. The implementation therefore performs a depth-first search over prime powers $p^e$ with $e\ge 2$, multiplies their contributions into the heavy product and into $c$, and prunes the branch as soon as either

$$\text{heavy product}>N$$

or the current $c$ can no longer divide any target $2k$ with $2\le k\le 10$.

Step 4: Count the squarefree completion under the product bound

Fix a heavy part

$$a=\prod_{i=1}^s p_i^{e_i}$$

and suppose that $2k=c\,2^t$. Then every corresponding integer has the form

$$n=a\,q_1q_2\cdots q_t,$$

with distinct primes

$$q_1<q_2<\cdots<q_t,\qquad q_j\notin\{p_1,\dots,p_s\},$$

and the product restriction

$$a\,q_1q_2\cdots q_t\le N.$$

When $t=0$, the heavy part itself contributes one solution. When $t=1$, we only need to count primes $q\le N/a$ while excluding the heavy primes already used. When $t\ge 2$, the implementation recursively chooses increasing primes and updates the remaining bound after every choice.

The crucial optimization is the base case $t=1$: instead of scanning primes one by one, the algorithm queries the prime-counting function $\pi(x)$. To make those queries fast, it precomputes $\pi(x)$ on the standard set of arguments

$$x\in\{1,2,\dots,\lfloor\sqrt N\rfloor\}\cup\left\{\left\lfloor\frac{N}{i}\right\rfloor:1\le i\le \lfloor\sqrt N\rfloor\right\}.$$

Step 5: Compute cumulative counts and recover each $Q_k$

Instead of rebuilding the search from scratch for every $k$, the implementation first computes the cumulative quantity

$$C(k)=\sum_{j=2}^k Q_j(N).$$

Every heavy-part pattern can contribute to several targets $2j$, so the search tree is reused while accumulating $C(k)$. After that, the exact count is recovered by the telescoping identity

$$Q_k(N)=C(k)-C(k-1).$$

This is especially effective here because the targets are very small and heavily overlap.

Worked Example: Why $Q_2(100)=51$

Here the target condition is $b(n)=4$. We classify possibilities by the heavy factor $c$.

If $c=1$, then $t=2$, so $n=pq$ with distinct primes and $pq\le 100$. These are the squarefree semiprimes, and there are $30$ of them.

If $c=2$, then $t=1$. The only way to get $c=2$ from a heavy exponent is $p^2$, so we count numbers of the form $p^2q$ with $p\ne q$ prime and $p^2q\le 100$. Grouping by $p$ gives

$$8+4+2+1=15,$$

coming from $p=2,3,5,7$.

If $c=4$, then $t=0$, so the heavy part alone must satisfy $b(n)=4$. The possible shapes are

$$p^3,\qquad p^4,\qquad p^2q^2.$$

Under $100$, these contribute $2$, $2$, and $2$ cases respectively, so this branch contributes $6$.

Therefore

$$Q_2(100)=30+15+6=51,$$

which matches the checkpoint built into the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same number-theoretic strategy. The C++ and Java versions contain the full combinatorial search, while the Python version delegates to the compiled implementation, so the observable algorithm is the same in all three languages.

First, the implementation generates primes up to about $2\sqrt N$ and builds a fast $\pi(x)$ table on the complementary floor-division values needed by the recursion. Next, it enumerates every admissible heavy part made of prime powers $p^e$ with $e\ge 2$, carrying both the numeric contribution of that part and its multiplicative contribution to $b(n)$.

For each heavy part, the implementation tests every target $2j$ up to the required bound. If the remaining quotient is a power of two, it knows exactly how many exponent-$1$ primes must still be appended. Those squarefree completions are counted recursively with increasing primes, and the final one-prime case is handled by the precomputed prime-counting table instead of explicit iteration.

Finally, cumulative values are cached, individual $Q_k(N)$ values are recovered by subtraction, and the required sum over $k=2,\dots,10$ is produced.

Complexity Analysis

Let $R=\lfloor\sqrt N\rfloor$. The prime sieve and the prime-counting table both use $O(R)$ memory and near $O(R\log\log R)$ time. The search over heavy parts is far smaller than a full factorization sweep, because every new heavy prime costs at least a square factor and every new heavy exponent multiplies the bi-unitary-divisor factor by at least $2$.

For this particular problem the squarefree-completion recursion is very shallow: since $2k\le 20$, the power-of-two quotient is at most $16$, so at most four exponent-$1$ primes ever need to be appended. In practice, the runtime is dominated by prime preprocessing and a modest number of prime-counting queries, while the memory usage remains $O(\sqrt N)$.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=861
Bi-unitary divisor: MathWorld - Biunitary Divisor
Unitary divisor: MathWorld - Unitary Divisor
Prime-counting function: Wikipedia - Prime-counting function
Multiplicative function: Wikipedia - Multiplicative function

Problem 861 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <limits>
#include <vector>

using i64 = std::int64_t;
using u64 = std::uint64_t;

class Solver {
public:
    explicit Solver(u64 n) : N_(n) {
        rN_ = static_cast<int>(std::sqrt(static_cast<long double>(N_)));
        while (static_cast<u64>(rN_ + 1) * static_cast<u64>(rN_ + 1) <= N_) {
            ++rN_;
        }
        while (static_cast<u64>(rN_) * static_cast<u64>(rN_) > N_) {
            --rN_;
        }
        primes_.push_back(2);
        get_primes(2 * rN_);
        pc_ = lucy();
    }

    u64 Q(int k) const {
        if (k < 2) {
            return 0;
        }
        return cumulative(k) - cumulative(k - 1);
    }

    u64 cumulative(int k) const {
        if (k < 2) {
            return 0;
        }
        if (k < static_cast<int>(cumulative_cache_.size()) &&
            cumulative_cache_[static_cast<std::size_t>(k)] != std::numeric_limits<u64>::max()) {
            return cumulative_cache_[static_cast<std::size_t>(k)];
        }

        const int two_k = 2 * k;
        u64 ans = 0;
        std::vector<int> cp;

        auto power_of_two_index = [](u64 x) -> int {
            if (x == 1) return 0;
            if (x == 2) return 1;
            if (x == 4) return 2;
            if (x == 8) return 3;
            if (x == 16) return 4;
            return -1;
        };

        dfs(0, N_, 1, 1, cp, two_k, [&](u64 pw, u64 cb, const std::vector<int>& cur_cp) {
            for (int k0 = 2; k0 <= k; ++k0) {
                const u64 two_k0 = static_cast<u64>(2 * k0);
                if (two_k0 % cb != 0) {
                    continue;
                }
                const int k2 = power_of_two_index(two_k0 / cb);
                if (k2 < 0) {
                    continue;
                }
                if (k2 == 0) {
                    ans += 1;
                } else {
                    ans += rec(cur_cp, 0, N_ / pw, pw, k2);
                }
            }
        });

        if (k < static_cast<int>(cumulative_cache_.size())) {
            cumulative_cache_[static_cast<std::size_t>(k)] = ans;
        }
        return ans;
    }

private:
    u64 N_;
    int rN_;
    std::vector<int> primes_;
    std::vector<i64> pc_;
    mutable std::vector<u64> cumulative_cache_{32, std::numeric_limits<u64>::max()};

    static bool contains_prime(const std::vector<int>& cp, int p) {
        for (int q : cp) {
            if (q == p) {
                return true;
            }
        }
        return false;
    }

    i64 pc_at(int idx) const {
        if (idx >= 0) {
            return pc_[static_cast<std::size_t>(idx)];
        }
        return pc_[pc_.size() - static_cast<std::size_t>(-idx)];
    }

    void get_primes(int limit) {
        std::vector<int> sieve(limit + 1, 0);
        for (int i = 2; i <= limit; i += 2) {
            sieve[static_cast<std::size_t>(i)] = 2;
        }

        for (int i = 3; i <= limit; i += 2) {
            if (sieve[static_cast<std::size_t>(i)] == 0) {
                primes_.push_back(i);
                sieve[static_cast<std::size_t>(i)] = i;
                if (static_cast<i64>(i) * static_cast<i64>(i) <= limit) {
                    for (int j = i * i; j <= limit; j += 2 * i) {
                        if (sieve[static_cast<std::size_t>(j)] == 0) {
                            sieve[static_cast<std::size_t>(j)] = i;
                        }
                    }
                }
            }
        }
    }

    std::vector<i64> lucy() const {
        const int r = rN_;
        std::vector<i64> S;
        S.reserve(static_cast<std::size_t>(2 * r + 1));
        for (int i = 0; i <= r; ++i) {
            S.push_back(static_cast<i64>(i - 1));
        }
        for (int i = r; i >= 1; --i) {
            S.push_back(static_cast<i64>(N_ / static_cast<u64>(i) - 1));
        }

        auto at = [&](int idx) -> i64& {
            if (idx >= 0) {
                return S[static_cast<std::size_t>(idx)];
            }
            return S[S.size() - static_cast<std::size_t>(-idx)];
        };

        for (int p = 2; p <= r; ++p) {
            if (S[static_cast<std::size_t>(p)] <= S[static_cast<std::size_t>(p - 1)]) {
                continue;
            }

            const i64 sp = S[static_cast<std::size_t>(p - 1)];
            const u64 p2 = static_cast<u64>(p) * static_cast<u64>(p);

            for (int i = 1; i <= r; ++i) {
                if (N_ / static_cast<u64>(i) < p2) {
                    break;
                }
                const u64 ip = static_cast<u64>(i) * static_cast<u64>(p);
                const u64 nip = N_ / ip;
                const int idx = (nip <= static_cast<u64>(r)) ? static_cast<int>(nip) : -static_cast<int>(ip);
                at(-i) -= (at(idx) - sp);
            }

            for (int i = r; i >= 1; --i) {
                if (static_cast<u64>(i) < p2) {
                    break;
                }
                S[static_cast<std::size_t>(i)] -= (S[static_cast<std::size_t>(i / p)] - sp);
            }
        }

        return S;
    }

    u64 rec(const std::vector<int>& cp, int i, u64 n, u64 cn, int k) const {
        if (k == 1) {
            const int idx = (cn < static_cast<u64>(rN_)) ? -static_cast<int>(cn) : static_cast<int>(N_ / cn);
            i64 ans = pc_at(idx) - i;
            if (ans < 0) {
                ans = 0;
            }
            const int start_prime =
                (i < static_cast<int>(primes_.size())) ? primes_[static_cast<std::size_t>(i)] : std::numeric_limits<int>::max();
            for (int p : cp) {
                if (static_cast<u64>(p) <= n && p >= start_prime) {
                    --ans;
                }
            }
            return (ans > 0) ? static_cast<u64>(ans) : 0;
        }

        u64 ans = 0;
        while (true) {
            if (i >= static_cast<int>(primes_.size())) {
                break;
            }
            const int p = primes_[static_cast<std::size_t>(i)];
            if (contains_prime(cp, p)) {
                ++i;
                continue;
            }
            const u64 nn = n / static_cast<u64>(p);
            if (nn < static_cast<u64>(p)) {
                break;
            }
            ans += rec(cp, i + 1, nn, cn * static_cast<u64>(p), k - 1);
            ++i;
        }
        return ans;
    }

    template <class F>
    void dfs(int i, u64 n, u64 c, u64 cb, std::vector<int>& cp, int two_k, F&& emit) const {
        emit(c, cb, cp);

        while (i < static_cast<int>(primes_.size())) {
            const u64 p = static_cast<u64>(primes_[static_cast<std::size_t>(i)]);
            if (p > std::numeric_limits<u64>::max() / p) {
                break;
            }

            u64 nn = n / (p * p);
            if (nn == 0) {
                break;
            }

            u64 cc = c * p * p;
            int e = 2;
            cp.push_back(static_cast<int>(p));

            while (true) {
                const u64 mult = static_cast<u64>(e + (e & 1));
                if (cb <= static_cast<u64>(two_k) / mult) {
                    dfs(i + 1, nn, cc, cb * mult, cp, two_k, emit);
                }

                if (nn < p) {
                    break;
                }
                nn /= p;
                if (cc > std::numeric_limits<u64>::max() / p) {
                    break;
                }
                cc *= p;
                ++e;
            }

            cp.pop_back();
            ++i;
        }
    }
};

int main() {
    {
        Solver s(100);
        assert(s.Q(2) == 51);
    }
    {
        Solver s(1'000'000);
        assert(s.Q(6) == 6'189);
    }

    Solver s(1'000'000'000'000ULL);
    u64 ans = 0;
    for (int k = 2; k <= 10; ++k) {
        ans += s.Q(k);
    }
    std::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

import java.util.ArrayList;
import java.util.HashMap;

public class Euler861 {

    static class Solver {
        long N;
        int rN;
        ArrayList<Integer> primes = new ArrayList<>();
        long[] pc;
        HashMap<Integer, Long> cumulativeCache = new HashMap<>();
        long currentAns = 0;

        Solver(long n) {
            this.N = n;
            this.rN = (int) Math.sqrt(N);
            while ((long) (rN + 1) * (rN + 1) <= N)
                rN++;
            while ((long) rN * rN > N)
                rN--;

            primes.add(2);
            getPrimes(2 * rN);
            pc = lucy();
        }

        long q(int k) {
            if (k < 2)
                return 0;
            return cumulative(k) - cumulative(k - 1);
        }

        int powerOfTwoIndex(long x) {
            if (x == 1)
                return 0;
            if (x == 2)
                return 1;
            if (x == 4)
                return 2;
            if (x == 8)
                return 3;
            if (x == 16)
                return 4;
            return -1;
        }

        long cumulative(int k) {
            if (k < 2)
                return 0;
            if (cumulativeCache.containsKey(k)) {
                return cumulativeCache.get(k);
            }

            int twoK = 2 * k;
            currentAns = 0;
            ArrayList<Integer> cp = new ArrayList<>();

            dfs(0, N, 1, 1, cp, twoK, k);

            cumulativeCache.put(k, currentAns);
            return currentAns;
        }

        void getPrimes(int limit) {
            if (limit < 2)
                return;
            int[] sieve = new int[limit + 1];
            for (int i = 2; i <= limit; i += 2)
                sieve[i] = 2;

            for (int i = 3; i <= limit; i += 2) {
                if (sieve[i] == 0) {
                    primes.add(i);
                    sieve[i] = i;
                    if ((long) i * i <= limit) {
                        for (int j = i * i; j <= limit; j += 2 * i) {
                            if (sieve[j] == 0)
                                sieve[j] = i;
                        }
                    }
                }
            }
        }

        long pcAt(int idx) {
            if (idx >= 0) {
                return pc[idx];
            }
            return pc[pc.length - (-idx)];
        }

        long[] lucy() {
            int r = rN;
            long[] S = new long[2 * r + 1];

            for (int i = 0; i <= r; ++i) {
                S[i] = i - 1;
            }
            for (int i = r; i >= 1; --i) {
                S[S.length - i] = (N / i) - 1;
            }

            for (int p = 2; p <= r; ++p) {
                if (S[p] <= S[p - 1])
                    continue;

                long sp = S[p - 1];
                long p2 = (long) p * p;

                for (int i = 1; i <= r; ++i) {
                    if (N / i < p2)
                        break;
                    long ip = (long) i * p;
                    long nip = N / ip;
                    int idx = (nip <= (long) r) ? (int) nip : -(int) ip;

                    int atI = S.length - i;
                    int atIdx = (idx >= 0) ? idx : (S.length - (-idx));

                    S[atI] -= (S[atIdx] - sp);
                }

                for (int i = r; i >= 1; --i) {
                    if ((long) i < p2)
                        break;
                    S[i] -= (S[i / p] - sp);
                }
            }

            return S;
        }

        long rec(ArrayList<Integer> cp, int i, long n, long cn, int k2) {
            if (k2 == 1) {
                int idx = (cn < (long) rN) ? -(int) cn : (int) (N / cn);
                long ans = pcAt(idx) - i;
                if (ans < 0)
                    ans = 0;

                int startPrime = (i < primes.size()) ? primes.get(i) : Integer.MAX_VALUE;
                for (int p : cp) {
                    if ((long) p <= n && p >= startPrime) {
                        ans--;
                    }
                }
                return ans > 0 ? ans : 0;
            }

            long ans = 0;
            while (true) {
                if (i >= primes.size())
                    break;
                int p = primes.get(i);
                if (cp.contains(p)) {
                    i++;
                    continue;
                }
                long nn = n / p;
                if (nn < p)
                    break;

                ans += rec(cp, i + 1, nn, cn * p, k2 - 1);
                i++;
            }
            return ans;
        }

        void dfs(int i, long n, long c, long cb, ArrayList<Integer> cp, int twoK, int kMax) {
            for (int k0 = 2; k0 <= kMax; ++k0) {
                long twoK0 = 2 * k0;
                if (twoK0 % cb != 0)
                    continue;

                int k2 = powerOfTwoIndex(twoK0 / cb);
                if (k2 < 0)
                    continue;

                if (k2 == 0) {
                    currentAns++;
                } else {
                    currentAns += rec(cp, 0, N / c, c, k2);
                }
            }

            while (i < primes.size()) {
                long p = primes.get(i);
                long nn = n / (p * p);
                if (nn == 0)
                    break;

                long cc = c * p * p;
                int e = 2;
                cp.add((int) p);

                while (true) {
                    long mult = e + (e & 1);
                    if (cb <= twoK / mult) {
                        dfs(i + 1, nn, cc, cb * mult, cp, twoK, kMax);
                    }

                    if (nn < p)
                        break;
                    nn /= p;
                    cc *= p;
                    e++;
                }

                cp.remove(cp.size() - 1);
                i++;
            }
        }
    }

    public static String solve() {
        Solver s = new Solver(1000000000000L);
        long ans = 0;
        for (int k = 2; k <= 10; ++k) {
            ans += s.q(k);
        }
        return Long.toString(ans);
    }

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