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

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