Detailed mathematical approach

Problem Summary

Fix a bound \(B\) and look only at primes \(p \le B\). The relevant primes split into the two nonzero residue classes modulo 3:

$\mathcal P_1=\{p \le B : p \text{ prime},\ p \equiv 1 \pmod{3}\}, \qquad \mathcal P_2=\{p \le B : p \text{ prime},\ p \equiv 2 \pmod{3}\}.$

The prime \(3\) is intentionally excluded, because it belongs to neither nonzero class and does not change the congruence that drives the count. For any integer \(n\), define

$\omega_1(n)=\#\{p\in \mathcal P_1 : p\mid n\},\qquad \omega_2(n)=\#\{p\in \mathcal P_2 : p\mid n\}.$

The computation counts

$F(N,B)=\#\{1\le n\le N : \omega_1(n)+2\omega_2(n)\equiv0\pmod3\}.$

Equivalently, it counts the integers for which \(\omega_1(n)\equiv\omega_2(n)\pmod3\). Only the presence or absence of relevant prime divisors matters, so prime exponents never enter the formula. The challenge is to evaluate this count for \(N=10^{18}\) and \(B=120\) without scanning all integers individually.

Mathematical Approach

The solution is a two-dimensional inclusion-exclusion over subsets of the relevant primes. The goal is to turn the modular condition on \(\omega_1(n)\) and \(\omega_2(n)\) into a weighted sum over squarefree products.

The relevant prime signature

Let

$\mathcal P_B=\mathcal P_1\cup\mathcal P_2.$

For a subset \(S\subseteq \mathcal P_B\), write

$d(S)=\prod_{p\in S} p,\qquad a(S)=\#(S\cap\mathcal P_1),\qquad b(S)=\#(S\cap\mathcal P_2).$

If \(d(S)\mid n\), then every prime in \(S\) divides \(n\). Since the definition depends only on which relevant primes divide \(n\), every integer is summarized by the signature \((\omega_1(n),\omega_2(n))\). Factors of \(3\), powers of primes already present, and prime divisors larger than \(B\) are all invisible to this signature.

Detecting the congruence class

Define the indicator

$G(x,y)=\begin{cases}1,& x+2y\equiv0\pmod3,\\0,& \text{otherwise}.\end{cases}$

A clean way to isolate this condition is with a cubic root-of-unity filter. If \(\zeta\neq1\) and \(\zeta^3=1\), then

$G(x,y)=\frac{1+\zeta^{x+2y}+\zeta^{2x+y}}{3}.$

This identity shows that only the counts modulo 3 matter. The implementations never use complex arithmetic directly, but this filter explains why the final weight table depends only on the two residue-class counts.

Binomial inversion turns the indicator into integer weights

Suppose an integer \(n\) has \(x\) relevant prime divisors from \(\mathcal P_1\) and \(y\) from \(\mathcal P_2\). We want weights \(w(a,b)\) such that

$G(x,y)=\sum_{a=0}^{x}\sum_{b=0}^{y}\binom{x}{a}\binom{y}{b}w(a,b).$

The combinatorial meaning is that we choose any \(a\) of the \(x\) class-\(1\) primes and any \(b\) of the \(y\) class-\(2\) primes. Each choice corresponds to one subset of the relevant prime divisors of \(n\).

Applying binomial inversion in both variables gives

$w(a,b)=\sum_{i=0}^{a}\sum_{j=0}^{b}(-1)^{a-i+b-j}\binom{a}{i}\binom{b}{j}G(i,j).$

Because \(G(i,j)\) vanishes unless \(i+2j\equiv0\pmod3\), this becomes

$w(a,b)=(-1)^{a+b}\sum_{\substack{0\le i\le a\\0\le j\le b\\i+2j\equiv0\!\!\!\pmod3}}(-1)^{i+j}\binom{a}{i}\binom{b}{j}.$

This is exactly the formula precomputed by the C++, Python, and Java implementations.

Collapsing the count to squarefree products

For each integer \(n\), let \(R_B(n)\subseteq\mathcal P_B\) be the set of relevant primes dividing \(n\), and let \(\mathbf{1}_B(n)\) be the indicator of \(B\)-trivisibility. The inversion formula implies

$\mathbf{1}_B(n)=\sum_{S\subseteq R_B(n)} w\big(a(S),b(S)\big).$

Now sum over all \(n\le N\) and swap the order of summation. A fixed subset \(S\) contributes once for every multiple of \(d(S)\), so its total contribution is

$\left\lfloor\frac{N}{d(S)}\right\rfloor.$

Therefore

$\boxed{F(N,B)=\sum_{S\subseteq\mathcal P_B} w\big(a(S),b(S)\big)\left\lfloor\frac{N}{d(S)}\right\rfloor.}$

This is the core reduction. The original search over \(1,\dots,N\) is replaced by a sum over subsets of primes up to \(B\), and only squarefree products appear because each prime is either selected once or not selected at all.

Worked example: \(B=10\)

For \(B=10\), the relevant primes are

$\mathcal P_1=\{7\},\qquad \mathcal P_2=\{2,5\}.$

The condition is \(\omega_1(n)+2\omega_2(n)\equiv0\pmod3\). Among the integers \(1\le n\le10\), the numbers \(1,3,9\) have signature \((0,0)\) and are counted. The numbers \(2,4,5,6,8\) have signature \((0,1)\), so they fail the congruence. The number \(7\) has signature \((1,0)\), and \(10\) has signature \((0,2)\); both are rejected as well. Hence

$F(10,10)=3,$

which matches the small checkpoint used by the implementations. The same reasoning explains \(F(10,4)=5\): if the only relevant prime is \(2\), then the counted integers are exactly the odd ones.

How the Code Works

The C++, Python, and Java implementations follow the derivation almost literally. They preprocess the prime classes, build the integer weight table, and then enumerate prime subsets recursively.

Generating the relevant primes

First, the implementation runs a sieve up to \(B\), removes the prime \(3\), and counts how many remaining primes lie in \(\mathcal P_1\) and \(\mathcal P_2\). These counts determine the dimensions of the weight table.

Building Pascal's triangle and the weight table

Next, it constructs binomial coefficients \(\binom{n}{k}\) by Pascal's rule. Using those coefficients, it fills the table \(w(a,b)\) with the inversion formula above. The root-of-unity viewpoint remains only a proof device; the actual computation stays entirely in integer arithmetic.

Recursive enumeration of squarefree products

The main recursion decides, for each relevant prime, whether that prime is excluded from the current subset or included in it. Along the recursion the implementation tracks the current product \(d(S)\), the number of selected primes from \(\mathcal P_1\), and the number selected from \(\mathcal P_2\). At a leaf it adds

$w\big(a(S),b(S)\big)\left\lfloor\frac{N}{d(S)}\right\rfloor$

to the answer. If including the next prime would make the product exceed \(N\), the branch is cut immediately, because every deeper extension would contribute zero. The C++ and Java implementations also verify a few small checkpoint values before printing the final result.

Complexity Analysis

Let \(P=|\mathcal P_B|\), \(M_1=|\mathcal P_1|\), \(M_2=|\mathcal P_2|\), and \(D=\max(M_1,M_2)\). Building the binomial table costs \(O(D^2)\) time and memory. Filling the weight table by the direct inversion formula costs \(O(M_1^2M_2^2)\) time and \(O(M_1M_2)\) memory, which is tiny for \(B=120\).

The recursive subset traversal has worst-case size \(O(2^P)\), since every relevant prime can be absent or present. In practice it is much smaller because branches are pruned as soon as their product exceeds \(N\). For the actual problem, \(P=29\), so the search tree is manageable and the algorithm avoids any dependence on the full scale of \(N=10^{18}\).

Footnotes and References

1. Problem page: Project Euler 967 - B-Trivisible Numbers
2. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
3. Binomial transform and binomial inversion: Wikipedia - Binomial transform
4. Modular arithmetic: Wikipedia - Modular arithmetic
5. Roots of unity: Wikipedia - Root of unity