Problem 428: Necklace of Circles
View on Project EulerProject Euler Problem 428 Solution
EulerSolve provides an optimized solution for Project Euler Problem 428, Necklace of Circles, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each positive integer \(n\), let \(T(n)\) be the number of positive integer triples \((a,b,c)\) with \(b \le n\) that satisfy the necklace condition from the problem statement. The direct geometric search is far too large for \(n=10^9\), so the solution used by the implementations turns the geometry into a divisor-counting problem and then evaluates that arithmetic description with cached summatory functions. Mathematical Approach The three implementations all start from the same arithmetic reduction of the necklace condition. The relevant chain lengths are \(k=3\), \(k=4\), and \(k=6\), leading to the Diophantine families $$k=3:\ (a-3b)(c-3b)=12b^2,$$ $$k=4:\ (a-b)(c-b)=2b^2,$$ $$k=6:\ (3a-b)(3c-b)=4b^2.$$ So the geometric configuration is reduced to counting admissible factorizations of \(2b^2\), \(4b^2\), \(12b^2\), and \(36b^2\), together with the local congruence restrictions forced by the primes \(2\) and \(3\). Step 1: Separate the Squarefree Part Away from \(2\) and \(3\) For primes \(p \ge 5\), only the squarefree divisor pattern matters. Define $$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$ the number of squarefree divisors of \(m\) that are coprime to \(6\)....
Detailed mathematical approach
Problem Summary
For each positive integer \(n\), let \(T(n)\) be the number of positive integer triples \((a,b,c)\) with \(b \le n\) that satisfy the necklace condition from the problem statement. The direct geometric search is far too large for \(n=10^9\), so the solution used by the implementations turns the geometry into a divisor-counting problem and then evaluates that arithmetic description with cached summatory functions.
Mathematical Approach
The three implementations all start from the same arithmetic reduction of the necklace condition. The relevant chain lengths are \(k=3\), \(k=4\), and \(k=6\), leading to the Diophantine families
$$k=3:\ (a-3b)(c-3b)=12b^2,$$
$$k=4:\ (a-b)(c-b)=2b^2,$$
$$k=6:\ (3a-b)(3c-b)=4b^2.$$
So the geometric configuration is reduced to counting admissible factorizations of \(2b^2\), \(4b^2\), \(12b^2\), and \(36b^2\), together with the local congruence restrictions forced by the primes \(2\) and \(3\).
Step 1: Separate the Squarefree Part Away from \(2\) and \(3\)
For primes \(p \ge 5\), only the squarefree divisor pattern matters. Define
$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$
the number of squarefree divisors of \(m\) that are coprime to \(6\). Its summatory function is
$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$
This function is the basic building block of the whole computation. Intuitively, it captures the contribution of all primes other than \(2\) and \(3\), while the local behavior at those two exceptional primes is handled separately.
Step 2: Count Squarefree Integers Coprime to \(6\)
To evaluate \(B(x)\), the implementations first count squarefree numbers by Möbius inversion:
$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$
Now let \(Q_6(x)\) denote the number of squarefree integers \(\le x\) that are coprime to \(6\). Every squarefree integer can be written uniquely as \(e m\) with \(e \in \{1,2,3,6\}\) and \((m,6)=1\), so
$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$
Rearranging gives the recursion used in the code:
$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$
Once \(Q_6\) is memoized, \(B(x)\) follows from the floor-sum identity above.
Step 3: Local Factors at the Primes \(2\) and \(3\)
After the squarefree part away from \(2\) and \(3\) has been isolated, the remaining arithmetic is encoded in four summatory functions. If \(B\) is the common base term from the previous step, the local combinations used by the implementations are
$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$
$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$
$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$
$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$
These are the exact linear combinations implemented in C++, Python, and Java. They represent the different local weights coming from the three Diophantine families once the odd squarefree part has been factored out.
Step 4: Turn the Outer Sums into Harmonic Floor Sums
For each \(D \in \{2,4,12,36\}\), write
$$P_D(x)=\sum_{m \le x} p_D(m).$$
The corresponding outer contribution is then
$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$
The program does not evaluate this sum term by term. Instead, it groups intervals on which \(\left\lfloor n/m \right\rfloor\) is constant. If \(\left\lfloor n/m \right\rfloor=v\) for all \(m \in [L,R]\), then that whole block contributes
$$v\left(P_D(R)-P_D(L-1)\right).$$
Because the quotient \(\left\lfloor n/m \right\rfloor\) takes only \(O(\sqrt{n})\) distinct values, this reduces a linear scan to a near square-root computation.
Step 5: The Modulo-\(3\) Character Correction
The case \(3 \nmid b\) needs one extra correction term. Introduce the nontrivial Dirichlet character modulo \(3\):
$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$
For squarefree numbers, the weighted prefix sum
$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m)$$
is again computed by Möbius inversion:
$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$
The inner character sum is especially simple, because
$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$
Now define the divisor-weighted function
$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$
Since \(\chi(d)=0\) whenever \(3 \mid d\), we have \(c(3m)=c(m)\). Therefore the summatory function restricted to \(3 \nmid m\) is obtained by subtraction, and the final correction term has the form
$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$
This is exactly the extra term that appears in the \(3 \nmid b\) branch of the implementation.
Step 6: Final Assembly by the \(3\)-Adic Valuation of \(b\)
The total is split according to the power of \(3\) dividing \(b\). For the case \(3 \nmid b\), the contribution is
$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$
If \(v_3(b)=t \ge 1\), the local multiplicity becomes \(2t-1\), so the higher \(3\)-adic layers contribute
$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right).$$
Putting everything together gives the exact formula implemented by the program:
$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$
How the Code Works
The C++, Python, and Java implementations follow the same pipeline. They precompute Möbius values up to \(\lfloor\sqrt{n}\rfloor\), memoize the count of squarefree integers coprime to \(6\), build the four local summatory functions above, evaluate each outer divisor sum through floor-division blocks, and then add the character correction together with the \(3\)-adic layers. As a sanity check, the arithmetic is verified on small checkpoints such as \(T(1)=9\), \(T(20)=732\), and \(T(3000)=438106\) before the full target is evaluated.
Complexity Analysis
The Möbius sieve up to \(\lfloor\sqrt{n}\rfloor\) costs \(O(\sqrt{n})\) time and memory. Every cached summatory function is queried only at values of the form \(\left\lfloor n/k \right\rfloor\), and there are only \(O(\sqrt{n})\) distinct quotients. In practice the running time is dominated by a small number of harmonic block sums plus memoized recursive lookups, which is fast enough for \(n=10^9\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=428
- Möbius function: Wikipedia - Möbius function
- Squarefree integer: Wikipedia - Squarefree integer
- Dirichlet character: Wikipedia - Dirichlet character
- Dirichlet hyperbola method: Wikipedia - Dirichlet hyperbola method
Problem 428 source code
C++
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>
using int64 = long long;
using i128 = __int128_t;
namespace {
int64 isqrt_ll(int64 n) {
int64 r = static_cast<int64>(std::sqrt(static_cast<long double>(n)));
while ((r + 1) * (r + 1) <= n) ++r;
while (r * r > n) --r;
return r;
}
class NecklaceCounter {
public:
explicit NecklaceCounter(int64 max_n) {
init_mu(max_n);
reserve_caches();
}
int64 T(int64 n) {
// k=3: (a-3b)(c-3b)=12 b^2, k=4: (a-b)(c-b)=2 b^2, k=6: (3a-b)(3c-b)=4 b^2.
int64 s2 = S_D(n, 2);
int64 s12 = S_D(n, 12);
int64 s4 = S_D(n, 4);
int64 s36_n3 = S_D(n / 3, 36);
i128 c_not3 = (static_cast<i128>(s4) - s36_n3 - S_h(n)) / 2;
i128 c_div3 = 0;
int64 m = n / 3;
int t = 1;
while (m > 0) {
i128 term = static_cast<i128>(S_D(m, 4)) - S_D(m / 3, 36);
c_div3 += static_cast<i128>(2 * t - 1) * term;
++t;
m /= 3;
}
i128 total = static_cast<i128>(s2) + s12 + c_not3 + c_div3;
return static_cast<int64>(total);
}
private:
std::vector<int> mu;
std::unordered_map<int64, int64> cache_c6;
std::unordered_map<int64, int64> cache_a0;
std::unordered_map<int64, int64> cache_chisq;
std::unordered_map<int64, int64> cache_fg;
std::unordered_map<int64, int64> cache_sh;
std::unordered_map<int64, int64> cache_sd2;
std::unordered_map<int64, int64> cache_sd4;
std::unordered_map<int64, int64> cache_sd12;
std::unordered_map<int64, int64> cache_sd36;
void reserve_caches() {
const size_t reserve = 1 << 18;
cache_c6.reserve(reserve);
cache_a0.reserve(reserve);
cache_chisq.reserve(reserve);
cache_fg.reserve(reserve);
cache_sh.reserve(reserve);
cache_sd2.reserve(reserve);
cache_sd4.reserve(reserve);
cache_sd12.reserve(reserve);
cache_sd36.reserve(reserve);
}
void init_mu(int64 max_n) {
int limit = static_cast<int>(isqrt_ll(max_n)) + 5;
mu.assign(limit + 1, 0);
std::vector<int> primes;
std::vector<int> is_comp(limit + 1, 0);
mu[1] = 1;
for (int i = 2; i <= limit; ++i) {
if (!is_comp[i]) {
primes.push_back(i);
mu[i] = -1;
}
for (int p : primes) {
int64 v = static_cast<int64>(i) * p;
if (v > limit) break;
is_comp[static_cast<int>(v)] = 1;
if (i % p == 0) {
mu[static_cast<int>(v)] = 0;
break;
} else {
mu[static_cast<int>(v)] = -mu[i];
}
}
}
}
int64 squarefree_count(int64 n) {
if (n <= 0) return 0;
int64 r = isqrt_ll(n);
i128 sum = 0;
for (int64 k = 1; k <= r; ++k) {
sum += static_cast<i128>(mu[static_cast<size_t>(k)]) * (n / (k * k));
}
return static_cast<int64>(sum);
}
int64 S_c6(int64 n) {
if (n <= 0) return 0;
auto it = cache_c6.find(n);
if (it != cache_c6.end()) return it->second;
int64 res = squarefree_count(n) - S_c6(n / 2) - S_c6(n / 3) - S_c6(n / 6);
cache_c6[n] = res;
return res;
}
int64 A0(int64 n) {
if (n <= 0) return 0;
auto it = cache_a0.find(n);
if (it != cache_a0.end()) return it->second;
// Sum floor(n/d) over squarefree d coprime to 6.
i128 res = 0;
int64 l = 1;
while (l <= n) {
int64 v = n / l;
int64 r = n / v;
res += static_cast<i128>(v) * (S_c6(r) - S_c6(l - 1));
l = r + 1;
}
int64 ans = static_cast<int64>(res);
cache_a0[n] = ans;
return ans;
}
int64 F_D(int64 n, int D) {
if (n <= 0) return 0;
if (D == 1) {
return A0(n) + A0(n / 2) + A0(n / 3) + A0(n / 6);
}
if (D == 2) {
return 2 * A0(n) + 2 * A0(n / 3);
}
if (D == 4) {
return 3 * A0(n) + 3 * A0(n / 3) - A0(n / 2) - A0(n / 6);
}
if (D == 12) {
return 6 * A0(n) - 2 * A0(n / 2);
}
if (D == 36) {
return 9 * A0(n) - 3 * A0(n / 2) - 3 * A0(n / 3) + A0(n / 6);
}
return 0;
}
int64 S_D(int64 n, int D) {
if (n <= 0) return 0;
std::unordered_map<int64, int64>* cache = nullptr;
if (D == 2) cache = &cache_sd2;
else if (D == 4) cache = &cache_sd4;
else if (D == 12) cache = &cache_sd12;
else if (D == 36) cache = &cache_sd36;
else return 0;
auto it = cache->find(n);
if (it != cache->end()) return it->second;
i128 res = 0;
int64 l = 1;
while (l <= n) {
int64 v = n / l;
int64 r = n / v;
res += static_cast<i128>(v) * (F_D(r, D) - F_D(l - 1, D));
l = r + 1;
}
int64 ans = static_cast<int64>(res);
(*cache)[n] = ans;
return ans;
}
int64 sum_chi_prefix(int64 n) {
if (n <= 0) return 0;
return (n % 3 == 1) ? 1 : 0;
}
int64 S_chisq(int64 n) {
if (n <= 0) return 0;
auto it = cache_chisq.find(n);
if (it != cache_chisq.end()) return it->second;
// Sum_{m<=n} mu^2(m) * chi(m) via square divisors.
int64 r = isqrt_ll(n);
i128 sum = 0;
for (int64 k = 1; k <= r; ++k) {
if (k % 3 == 0) continue;
sum += static_cast<i128>(mu[static_cast<size_t>(k)]) *
sum_chi_prefix(n / (k * k));
}
int64 ans = static_cast<int64>(sum);
cache_chisq[n] = ans;
return ans;
}
int64 F_G(int64 n) {
if (n <= 0) return 0;
auto it = cache_fg.find(n);
if (it != cache_fg.end()) return it->second;
i128 res = 0;
int64 l = 1;
while (l <= n) {
int64 v = n / l;
int64 r = n / v;
res += static_cast<i128>(v) * (S_chisq(r) - S_chisq(l - 1));
l = r + 1;
}
int64 ans = static_cast<int64>(res);
cache_fg[n] = ans;
return ans;
}
int64 F1(int64 n) {
if (n <= 0) return 0;
return F_G(n) - F_G(n / 3);
}
int64 S_h(int64 n) {
if (n <= 0) return 0;
auto it = cache_sh.find(n);
if (it != cache_sh.end()) return it->second;
// Convolution with the mod-3 character yields a simple filter on floor(n/d) mod 3.
i128 res = 0;
int64 l = 1;
while (l <= n) {
int64 v = n / l;
int64 r = n / v;
if (v % 3 == 1) {
res += static_cast<i128>(F1(r) - F1(l - 1));
}
l = r + 1;
}
int64 ans = static_cast<int64>(res);
cache_sh[n] = ans;
return ans;
}
};
} // namespace
int main() {
const int64 target = 1000000000LL;
NecklaceCounter solver(target);
struct Check { int64 n; int64 expected; };
const Check checks[] = {
{1, 9},
{20, 732},
{3000, 438106},
};
for (const auto& chk : checks) {
int64 got = solver.T(chk.n);
if (got != chk.expected) {
std::cerr << "Validation failed for n=" << chk.n
<< ": got " << got << ", expected " << chk.expected << '\n';
return 1;
}
std::cout << "Validation passed for n=" << chk.n << '\n';
}
const int64 answer = solver.T(target);
std::cout << answer << '\n';
std::cout << "Answer: " << answer << '\n';
return 0;
}
Python
import math
def solve():
target = 1000000000
def isqrt(n):
r = int(math.isqrt(n))
while (r+1)*(r+1) <= n: r += 1
while r*r > n: r -= 1
return r
limit = isqrt(target) + 5
mu = [0] * (limit + 1)
mu[1] = 1
primes = []
is_comp = [0] * (limit + 1)
for i in range(2, limit + 1):
if not is_comp[i]:
primes.append(i); mu[i] = -1
for p in primes:
v = i * p
if v > limit: break
is_comp[v] = 1
if i % p == 0: mu[v] = 0; break
mu[v] = -mu[i]
def sqfree_count(n):
if n <= 0: return 0
r = isqrt(n); s = 0
for k in range(1, r + 1): s += mu[k] * (n // (k*k))
return s
cache_c6 = {}
def S_c6(n):
if n <= 0: return 0
if n in cache_c6: return cache_c6[n]
r = sqfree_count(n) - S_c6(n//2) - S_c6(n//3) - S_c6(n//6)
cache_c6[n] = r; return r
cache_a0 = {}
def A0(n):
if n <= 0: return 0
if n in cache_a0: return cache_a0[n]
res = 0; l = 1
while l <= n:
v = n // l; r = n // v
res += v * (S_c6(r) - S_c6(l - 1))
l = r + 1
cache_a0[n] = res; return res
def F_D(n, D):
if n <= 0: return 0
if D == 2: return 2*A0(n) + 2*A0(n//3)
if D == 4: return 3*A0(n) + 3*A0(n//3) - A0(n//2) - A0(n//6)
if D == 12: return 6*A0(n) - 2*A0(n//2)
if D == 36: return 9*A0(n) - 3*A0(n//2) - 3*A0(n//3) + A0(n//6)
return 0
cache_sd = {}
def S_D(n, D):
if n <= 0: return 0
k = (n, D)
if k in cache_sd: return cache_sd[k]
res = 0; l = 1
while l <= n:
v = n // l; r = n // v
res += v * (F_D(r, D) - F_D(l - 1, D))
l = r + 1
cache_sd[k] = res; return res
def chi_prefix(n):
if n <= 0: return 0
return 1 if n % 3 == 1 else 0
cache_chisq = {}
def S_chisq(n):
if n <= 0: return 0
if n in cache_chisq: return cache_chisq[n]
r = isqrt(n); s = 0
for k in range(1, r + 1):
if k % 3 == 0: continue
s += mu[k] * chi_prefix(n // (k*k))
cache_chisq[n] = s; return s
cache_fg = {}
def F_G(n):
if n <= 0: return 0
if n in cache_fg: return cache_fg[n]
res = 0; l = 1
while l <= n:
v = n // l; r = n // v
res += v * (S_chisq(r) - S_chisq(l - 1))
l = r + 1
cache_fg[n] = res; return res
def F1(n):
if n <= 0: return 0
return F_G(n) - F_G(n // 3)
cache_sh = {}
def S_h(n):
if n <= 0: return 0
if n in cache_sh: return cache_sh[n]
res = 0; l = 1
while l <= n:
v = n // l; r = n // v
if v % 3 == 1:
res += F1(r) - F1(l - 1)
l = r + 1
cache_sh[n] = res; return res
s2 = S_D(target, 2)
s12 = S_D(target, 12)
s4 = S_D(target, 4)
s36_n3 = S_D(target // 3, 36)
c_not3 = (s4 - s36_n3 - S_h(target)) // 2
c_div3 = 0; m = target // 3; t = 1
while m > 0:
term = S_D(m, 4) - S_D(m // 3, 36)
c_div3 += (2 * t - 1) * term
t += 1; m //= 3
return str(s2 + s12 + c_not3 + c_div3)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class Euler428 {
static long isqrt(long n) {
long r = (long) Math.sqrt((double) n);
while ((r + 1) * (r + 1) <= n)
r++;
while (r * r > n)
r--;
return r;
}
static class NecklaceCounter {
long maxN;
int[] mu;
Map<Long, Long> cacheC6 = new HashMap<>();
Map<Long, Long> cacheA0 = new HashMap<>();
Map<Long, Long> cacheChisq = new HashMap<>();
Map<Long, Long> cacheFg = new HashMap<>();
Map<Long, Long> cacheSh = new HashMap<>();
Map<Long, Long> cacheSd2 = new HashMap<>();
Map<Long, Long> cacheSd4 = new HashMap<>();
Map<Long, Long> cacheSd12 = new HashMap<>();
Map<Long, Long> cacheSd36 = new HashMap<>();
NecklaceCounter(long maxN) {
this.maxN = maxN;
int limit = (int) isqrt(maxN) + 5;
mu = new int[limit + 1];
initMu(limit);
}
void initMu(int limit) {
int[] isComp = new int[limit + 1];
List<Integer> primes = new ArrayList<>();
mu[1] = 1;
for (int i = 2; i <= limit; i++) {
if (isComp[i] == 0) {
primes.add(i);
mu[i] = -1;
}
for (int p : primes) {
long v = (long) i * p;
if (v > limit)
break;
isComp[(int) v] = 1;
if (i % p == 0) {
mu[(int) v] = 0;
break;
} else {
mu[(int) v] = -mu[i];
}
}
}
}
long squarefreeCount(long n) {
if (n <= 0)
return 0;
long r = isqrt(n);
long s = 0;
for (long k = 1; k <= r; k++) {
s += (long) mu[(int) k] * (n / (k * k));
}
return s;
}
long Sc6(long n) {
if (n <= 0)
return 0;
if (cacheC6.containsKey(n))
return cacheC6.get(n);
long res = squarefreeCount(n) - Sc6(n / 2) - Sc6(n / 3) - Sc6(n / 6);
cacheC6.put(n, res);
return res;
}
long A0(long n) {
if (n <= 0)
return 0;
if (cacheA0.containsKey(n))
return cacheA0.get(n);
long res = 0;
long l = 1;
while (l <= n) {
long v = n / l;
long r = n / v;
res += v * (Sc6(r) - Sc6(l - 1));
l = r + 1;
}
cacheA0.put(n, res);
return res;
}
long FD(long n, int D) {
if (n <= 0)
return 0;
if (D == 1)
return A0(n) + A0(n / 2) + A0(n / 3) + A0(n / 6);
if (D == 2)
return 2 * A0(n) + 2 * A0(n / 3);
if (D == 4)
return 3 * A0(n) + 3 * A0(n / 3) - A0(n / 2) - A0(n / 6);
if (D == 12)
return 6 * A0(n) - 2 * A0(n / 2);
if (D == 36)
return 9 * A0(n) - 3 * A0(n / 2) - 3 * A0(n / 3) + A0(n / 6);
return 0;
}
long SD(long n, int D) {
if (n <= 0)
return 0;
Map<Long, Long> cache;
if (D == 2)
cache = cacheSd2;
else if (D == 4)
cache = cacheSd4;
else if (D == 12)
cache = cacheSd12;
else if (D == 36)
cache = cacheSd36;
else
return 0;
if (cache.containsKey(n))
return cache.get(n);
long res = 0;
long l = 1;
while (l <= n) {
long v = n / l;
long r = n / v;
res += v * (FD(r, D) - FD(l - 1, D));
l = r + 1;
}
cache.put(n, res);
return res;
}
long sumChiPrefix(long n) {
if (n <= 0)
return 0;
return (n % 3 == 1) ? 1 : 0;
}
long Schisq(long n) {
if (n <= 0)
return 0;
if (cacheChisq.containsKey(n))
return cacheChisq.get(n);
long r = isqrt(n);
long s = 0;
for (long k = 1; k <= r; k++) {
if (k % 3 == 0)
continue;
s += (long) mu[(int) k] * sumChiPrefix(n / (k * k));
}
cacheChisq.put(n, s);
return s;
}
long FG(long n) {
if (n <= 0)
return 0;
if (cacheFg.containsKey(n))
return cacheFg.get(n);
long res = 0;
long l = 1;
while (l <= n) {
long v = n / l;
long r = n / v;
res += v * (Schisq(r) - Schisq(l - 1));
l = r + 1;
}
cacheFg.put(n, res);
return res;
}
long F1(long n) {
if (n <= 0)
return 0;
return FG(n) - FG(n / 3);
}
long Sh(long n) {
if (n <= 0)
return 0;
if (cacheSh.containsKey(n))
return cacheSh.get(n);
long res = 0;
long l = 1;
while (l <= n) {
long v = n / l;
long r = n / v;
if (v % 3 == 1) {
res += F1(r) - F1(l - 1);
}
l = r + 1;
}
cacheSh.put(n, res);
return res;
}
long T(long n) {
long s2 = SD(n, 2);
long s12 = SD(n, 12);
long s4 = SD(n, 4);
long s36n3 = SD(n / 3, 36);
long cNot3 = (s4 - s36n3 - Sh(n)) / 2;
long cDiv3 = 0;
long m = n / 3;
long t = 1;
while (m > 0) {
long term = SD(m, 4) - SD(m / 3, 36);
cDiv3 += (2 * t - 1) * term;
t++;
m /= 3;
}
return s2 + s12 + cNot3 + cDiv3;
}
}
public static String solve() {
long target = 1000000000L;
NecklaceCounter solver = new NecklaceCounter(target);
return Long.toString(solver.T(target));
}
public static void main(String[] args) {
System.out.println(solve());
}
}