Problem 769: Binary Quadratic Form II
View on Project EulerProject Euler Problem 769 Solution
EulerSolve provides an optimized solution for Project Euler Problem 769, Binary Quadratic Form II, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The solution works with a canonical parametrization by primitive integer pairs \((p,q)\) with \(q\ge 0\) and \(\gcd(p,q)=1\). In that parametrization, the size attached to a solution is the discriminant-\(13\) binary quadratic form $$z=\left|3p^2-7pq+3q^2\right|,$$ and we must count all admissible primitive pairs for which \(z\le N\). A second quadratic coordinate must stay positive as well, so the ratio \(p/q\) is not free: only three regions survive after the sign and symmetry reductions. The isolated case \(q=0\) gives \(z=3\). Every other contribution comes from \(q\ge 1\), which is why the implementations loop over \(q\) and count valid \(k\)-intervals instead of scanning the original search space directly. Mathematical Approach Introduce the two quadratic forms $$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$ The counted size is \(|Q(p,q)|\), while admissibility requires \(X(p,q)>0\), \(|p|>q\), and \(\gcd(p,q)=1\). The code turns those conditions into interval counts in a new variable \(k\). Step 1: Replace \(p\) by a positive offset \(k\) Because \(|p|>q\), every admissible pair with \(q>0\) can be written in exactly one of the forms $$p=q+k,\qquad p=-q-k.$$ with \(k\ge 1\)....
Detailed mathematical approach
Problem Summary
The solution works with a canonical parametrization by primitive integer pairs \((p,q)\) with \(q\ge 0\) and \(\gcd(p,q)=1\). In that parametrization, the size attached to a solution is the discriminant-\(13\) binary quadratic form
$$z=\left|3p^2-7pq+3q^2\right|,$$
and we must count all admissible primitive pairs for which \(z\le N\). A second quadratic coordinate must stay positive as well, so the ratio \(p/q\) is not free: only three regions survive after the sign and symmetry reductions.
The isolated case \(q=0\) gives \(z=3\). Every other contribution comes from \(q\ge 1\), which is why the implementations loop over \(q\) and count valid \(k\)-intervals instead of scanning the original search space directly.
Mathematical Approach
Introduce the two quadratic forms
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
The counted size is \(|Q(p,q)|\), while admissibility requires \(X(p,q)>0\), \(|p|>q\), and \(\gcd(p,q)=1\). The code turns those conditions into interval counts in a new variable \(k\).
Step 1: Replace \(p\) by a positive offset \(k\)
Because \(|p|>q\), every admissible pair with \(q>0\) can be written in exactly one of the forms
$$p=q+k,\qquad p=-q-k.$$
with \(k\ge 1\). Coprimality is preserved under this substitution:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
So after fixing \(q\), the arithmetic condition becomes simply \(\gcd(q,k)=1\).
Step 2: Split the search into the three surviving regions
For \(p=q+k\), the auxiliary form becomes
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
If we write
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$$
then \(X>0\) forces either \(k<\alpha q\) or \(k>\beta q\). These are the two branches called region A and region D:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
For \(p=-q-k\), we get
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
so positivity is automatic, and the size becomes
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
Thus the entire count is the sum of region A, region C, region D, and the single base solution \(z=3\).
Step 3: Count coprime \(k\) by Möbius inversion
For fixed \(q\), define
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
Using inclusion-exclusion over the prime divisors of \(q\),
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
Only squarefree divisors matter, so after factoring \(q\) once, the implementations generate all \(2^{\omega(q)}\) squarefree divisors together with their Möbius signs. Any interval \([L,U]\) is then counted by
$$C_q(U)-C_q(L-1).$$
Step 4: Remove the non-canonical branch modulo \(13\)
The discriminant-\(13\) form has a simple congruence:
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
Therefore
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
The implementations keep only the canonical branch with \(13\nmid Q(p,q)\), so one residue class must be removed:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
When \(q\equiv 0\pmod{13}\), coprimality already forbids \(k\equiv 0\pmod{13}\), so there is no extra subtraction. Otherwise the forbidden class is counted with the same Möbius table, but now on an arithmetic progression modulo \(13\).
Step 5: Turn the size bound \(z\le N\) into explicit endpoints
Each region contributes an interval in \(k\).
Region A starts with \(1\le k\le \lfloor \alpha q\rfloor\). The polynomial \(z_A(q,k)=q^2+qk-3k^2\) is concave, so the inequality \(z_A\le N\) can fail only on a middle interval. Solving \(z_A=N\) gives
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
If \(13q^2\le 12N\), the whole region survives; otherwise the integers between the two roots must be removed.
Region C is monotone increasing, so its upper bound comes from
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$$
Region D is also monotone once \(k>\beta q\), so it uses
$$k\ge \lfloor \beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
The implementations compute these bounds from integer square roots and then adjust the endpoints by a few exact checks, eliminating rounding errors.
Worked Example: \(N=100\) and \(q=1\)
Here \(\lfloor \alpha\rfloor=0\), so region A is empty.
Region C satisfies
$$13+13k+3k^2\le 100,$$
which gives \(k=1,2,3\), hence the contributions \(29,51,79\).
Region D starts at \(k=\lfloor \beta\rfloor+1=4\). There we need
$$3k^2-k-1\le 100,$$
so \(k=4,5\), giving \(43\) and \(69\).
For \(q=1\), the forbidden classes are \(k\equiv 11\pmod{13}\) in regions A and D, and \(k\equiv 0\pmod{13}\) in region C, so no subtraction occurs in this layer. Together with the isolated case \(q=0\), this shows exactly how the counting proceeds.
How the Code Works
The C++, Python, and Java implementations first set \(Q_{\max}=\lfloor\sqrt N\rfloor\) and build a smallest-prime-factor table up to that limit. This lets them factor every \(q\) quickly and generate all squarefree divisors needed for the Möbius sums.
For each \(q\), the implementation evaluates the three regions separately. Region A is handled as a prefix count up to \(\lfloor\alpha q\rfloor\) minus the middle interval where \(z_A>N\). Regions C and D are monotone, so they are counted with direct prefix differences. In every case the coprime filter and, when necessary, the forbidden residue class modulo \(13\) are applied through the same divisor table.
The C++ implementation additionally splits the outer \(q\)-range across worker threads and accumulates partial sums in parallel. The Python and Java implementations perform the same arithmetic serially. After all \(q\ge 1\) contributions are finished, the isolated primitive case \(z=3\) is added once \(N\ge 3\).
Complexity Analysis
Let \(Q=\lfloor\sqrt N\rfloor\). Building the smallest-prime-factor sieve costs \(O(Q\log\log Q)\) time and \(O(Q)\) memory. For a fixed \(q\), the number of squarefree divisors is \(2^{\omega(q)}\), and each region uses only a constant number of sums over that table.
So the total running time is
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
which is close to \(O(Q\log Q)\) on average and behaves near-linearly in practice. Memory usage stays \(O(Q)\), with only small per-\(q\) temporary arrays besides the sieve.
Footnotes and References
- Problem page: https://projecteuler.net/problem=769
- Binary quadratic form: Wikipedia — Binary quadratic form
- Möbius inversion formula: Wikipedia — Möbius inversion formula
- Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
- Greatest common divisor: Wikipedia — Greatest common divisor
Problem 769 source code
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
using namespace std;
namespace {
uint64_t isqrt_u64(uint64_t x) {
long double r = sqrtl(static_cast<long double>(x));
uint64_t y = static_cast<uint64_t>(r);
while ((y + 1) * (y + 1) <= x) ++y;
while (y * y > x) --y;
return y;
}
vector<int> build_spf(int limit) {
vector<int> spf(limit + 1, 0);
if (limit >= 1) spf[1] = 1;
for (int i = 2; i <= limit; ++i) {
if (spf[i] == 0) {
spf[i] = i;
if (1LL * i * i <= limit) {
for (int j = i * i; j <= limit; j += i) {
if (spf[j] == 0) spf[j] = i;
}
}
}
}
return spf;
}
struct Counter {
const vector<int>& spf;
long double alpha;
long double beta;
int inv13[13];
explicit Counter(const vector<int>& spf_ref) : spf(spf_ref) {
alpha = 2.0L - sqrtl(3.0L);
beta = 2.0L + sqrtl(3.0L);
inv13[0] = 0;
for (int i = 1; i < 13; ++i) {
for (int j = 1; j < 13; ++j) {
if ((i * j) % 13 == 1) {
inv13[i] = j;
break;
}
}
}
}
void build_squarefree_divs(int q, vector<int>& divs, vector<int>& mus) const {
divs.clear();
mus.clear();
divs.push_back(1);
mus.push_back(1);
int x = q;
while (x > 1) {
int p = spf[x];
while (x % p == 0) x /= p;
int current = static_cast<int>(divs.size());
for (int i = 0; i < current; ++i) {
divs.push_back(divs[i] * p);
mus.push_back(-mus[i]);
}
}
}
long long coprime_count(const vector<int>& divs,
const vector<int>& mus,
long long k_max) const {
if (k_max <= 0) return 0;
long long total = 0;
for (size_t i = 0; i < divs.size(); ++i) {
total += static_cast<long long>(mus[i]) * (k_max / divs[i]);
}
return total;
}
long long coprime_count_residue(const vector<int>& divs,
const vector<int>& mus,
long long k_max,
int residue) const {
if (k_max <= 0) return 0;
long long total = 0;
for (size_t i = 0; i < divs.size(); ++i) {
int d = divs[i];
int d_mod = d % 13;
int inv = inv13[d_mod];
int t0 = (residue * inv) % 13;
if (t0 == 0) continue;
long long first = 1LL * d * t0;
if (first > k_max) continue;
total += static_cast<long long>(mus[i]) *
(1 + (k_max - first) / (13LL * d));
}
return total;
}
long long compute(uint64_t N, int threads) const {
uint64_t q_limit = isqrt_u64(N);
if (q_limit == 0) return 0;
if (threads <= 0) {
unsigned hw = thread::hardware_concurrency();
threads = hw == 0 ? 1 : static_cast<int>(hw);
}
if (q_limit < static_cast<uint64_t>(threads)) {
threads = static_cast<int>(q_limit);
}
vector<long long> partial(threads, 0);
uint64_t chunk = (q_limit + threads - 1) / threads;
auto worker = [&](int idx, uint64_t start, uint64_t end) {
vector<int> divs;
vector<int> mus;
divs.reserve(64);
mus.reserve(64);
long long local = 0;
for (uint64_t q = start; q <= end; ++q) {
build_squarefree_divs(static_cast<int>(q), divs, mus);
long long q_ll = static_cast<long long>(q);
uint64_t q2 = static_cast<uint64_t>(q_ll) * q_ll;
int q_mod13 = static_cast<int>(q % 13);
int residue = q_mod13 == 0 ? 0 : (13 - (2 * q_mod13) % 13) % 13;
auto count_with_residue = [&](long long k_max) -> long long {
if (k_max <= 0) return 0;
long long base = coprime_count(divs, mus, k_max);
if (q_mod13 != 0) {
base -= coprime_count_residue(divs, mus, k_max, residue);
}
return base;
};
auto count_exclude13 = [&](long long k_max) -> long long {
if (k_max <= 0) return 0;
long long base = coprime_count(divs, mus, k_max);
if (q_mod13 != 0) {
base -= coprime_count(divs, mus, k_max / 13);
}
return base;
};
auto zA = [&](long long k) -> long long {
__int128 val = static_cast<__int128>(q2) +
static_cast<__int128>(q_ll) * k -
3 * static_cast<__int128>(k) * k;
return static_cast<long long>(val);
};
auto zC = [&](long long k) -> long long {
__int128 val = 13 * static_cast<__int128>(q2) +
13 * static_cast<__int128>(q_ll) * k +
3 * static_cast<__int128>(k) * k;
return static_cast<long long>(val);
};
auto zD = [&](long long k) -> long long {
__int128 val = 3 * static_cast<__int128>(k) * k -
static_cast<__int128>(q_ll) * k -
static_cast<__int128>(q2);
return static_cast<long long>(val);
};
// Region A: 1 < p/q < 3 - sqrt(3)
long long k_max = static_cast<long long>(floorl(alpha * q));
if (k_max > 0) {
while (k_max > 0) {
__int128 x_val = static_cast<__int128>(q2) -
4 * static_cast<__int128>(q_ll) * k_max +
static_cast<__int128>(k_max) * k_max;
if (x_val > 0) break;
--k_max;
}
if (k_max > 0) {
long long base = count_with_residue(k_max);
__int128 disc128 = 13 * static_cast<__int128>(q2) -
12 * static_cast<__int128>(N);
if (disc128 > 0) {
uint64_t disc = static_cast<uint64_t>(disc128);
uint64_t s = isqrt_u64(disc);
long long L = (q_ll - static_cast<long long>(s) + 5) / 6;
long long R = (q_ll + static_cast<long long>(s)) / 6;
while (L <= k_max && zA(L) <= static_cast<long long>(N)) ++L;
while (R >= 1 && zA(R) <= static_cast<long long>(N)) --R;
if (L <= R && L <= k_max && R >= 1) {
long long L2 = max(1LL, L);
long long R2 = min(k_max, R);
if (L2 <= R2) {
long long forbidden = count_with_residue(R2) -
count_with_residue(L2 - 1);
base -= forbidden;
}
}
}
local += base;
}
}
// Common sqrt for regions C and D.
__int128 disc_plus128 = 13 * static_cast<__int128>(q2) +
12 * static_cast<__int128>(N);
uint64_t disc_plus = static_cast<uint64_t>(disc_plus128);
uint64_t s_plus = isqrt_u64(disc_plus);
// Region C: p/q < -1
long long k_max_c = (static_cast<long long>(s_plus) - 13LL * q_ll) / 6;
if (k_max_c > 0) {
while (k_max_c > 0 && zC(k_max_c) > static_cast<long long>(N)) {
--k_max_c;
}
while (zC(k_max_c + 1) <= static_cast<long long>(N)) {
++k_max_c;
}
if (k_max_c > 0) {
local += count_exclude13(k_max_c);
}
}
// Region D: p/q > 3 + sqrt(3)
long long k_min_d = static_cast<long long>(floorl(beta * q)) + 1;
if (k_min_d < 1) k_min_d = 1;
while (true) {
__int128 x_val = static_cast<__int128>(q2) -
4 * static_cast<__int128>(q_ll) * k_min_d +
static_cast<__int128>(k_min_d) * k_min_d;
if (x_val > 0) break;
++k_min_d;
}
long long k_max_d = (static_cast<long long>(s_plus) + q_ll) / 6;
while (k_max_d > 0 && zD(k_max_d) > static_cast<long long>(N)) {
--k_max_d;
}
while (zD(k_max_d + 1) <= static_cast<long long>(N)) {
++k_max_d;
}
if (k_min_d <= k_max_d) {
local += count_with_residue(k_max_d) -
count_with_residue(k_min_d - 1);
}
}
partial[idx] = local;
};
vector<thread> workers;
workers.reserve(threads);
for (int t = 0; t < threads; ++t) {
uint64_t start = t * chunk + 1;
uint64_t end = min(q_limit, (t + 1) * chunk);
if (start > end) {
partial[t] = 0;
continue;
}
workers.emplace_back(worker, t, start, end);
}
for (auto& th : workers) th.join();
long long total = 0;
for (long long v : partial) total += v;
if (N >= 3) total += 1; // (1,1,3)
return total;
}
};
} // namespace
int main() {
const uint64_t N = 100000000000000ULL;
uint64_t q_limit = isqrt_u64(N);
vector<int> spf = build_spf(static_cast<int>(q_limit));
Counter counter(spf);
// Validation checkpoints.
long long check1 = counter.compute(1000, 1);
long long check2 = counter.compute(1000000, 1);
if (check1 != 142 || check2 != 142463) {
cerr << "Validation failed: C(1e3)=" << check1
<< " C(1e6)=" << check2 << "\n";
return 1;
}
int threads = static_cast<int>(thread::hardware_concurrency());
long long answer = counter.compute(N, threads);
cout << answer << "\n";
return 0;
}
Python
import math
def isqrt(x):
if x < 0:
raise ValueError("isqrt of negative")
if x == 0:
return 0
return int(math.isqrt(x))
def build_spf(limit):
spf = [0] * (limit + 1)
if limit >= 1:
spf[1] = 1
for i in range(2, limit + 1):
if spf[i] == 0:
spf[i] = i
if i * i <= limit:
for j in range(i * i, limit + 1, i):
if spf[j] == 0:
spf[j] = i
return spf
class Counter:
def __init__(self, limit):
self.spf = build_spf(limit)
self.alpha = 2.0 - math.sqrt(3.0)
self.beta = 2.0 + math.sqrt(3.0)
self.inv13 = [0] * 13
for i in range(1, 13):
for j in range(1, 13):
if (i * j) % 13 == 1:
self.inv13[i] = j
break
def build_squarefree_divs(self, q):
divs = [1]
mus = [1]
x = q
while x > 1:
p = self.spf[x]
while x % p == 0:
x //= p
current = len(divs)
for i in range(current):
divs.append(divs[i] * p)
mus.append(-mus[i])
return divs, mus
def coprime_count(self, divs, mus, k_max):
if k_max <= 0: return 0
total = 0
for d, mu in zip(divs, mus):
total += mu * (k_max // d)
return total
def coprime_count_residue(self, divs, mus, k_max, residue):
if k_max <= 0: return 0
total = 0
for d, mu in zip(divs, mus):
d_mod = d % 13
inv = self.inv13[d_mod]
t0 = (residue * inv) % 13
if t0 == 0: continue
first = d * t0
if first > k_max: continue
total += mu * (1 + (k_max - first) // (13 * d))
return total
def compute(self, N):
q_limit = isqrt(N)
if q_limit == 0: return 0
total = 0
for q in range(1, q_limit + 1):
divs, mus = self.build_squarefree_divs(q)
q2 = q * q
q_mod13 = q % 13
residue = 0 if q_mod13 == 0 else (13 - (2 * q_mod13) % 13) % 13
def count_with_residue(k_max):
if k_max <= 0: return 0
base = self.coprime_count(divs, mus, k_max)
if q_mod13 != 0:
base -= self.coprime_count_residue(divs, mus, k_max, residue)
return base
def count_exclude13(k_max):
if k_max <= 0: return 0
base = self.coprime_count(divs, mus, k_max)
if q_mod13 != 0:
base -= self.coprime_count(divs, mus, k_max // 13)
return base
def zA(k):
return q2 + q * k - 3 * k * k
def zC(k):
return 13 * q2 + 13 * q * k + 3 * k * k
def zD(k):
return 3 * k * k - q * k - q2
# Region A
k_max = int(math.floor(self.alpha * q))
if k_max > 0:
while k_max > 0:
x_val = q2 - 4 * q * k_max + k_max * k_max
if x_val > 0: break
k_max -= 1
if k_max > 0:
base = count_with_residue(k_max)
disc = 13 * q2 - 12 * N
if disc > 0:
s = isqrt(disc)
L = (q - s + 5) // 6
R = (q + s) // 6
while L <= k_max and zA(L) <= N: L += 1
while R >= 1 and zA(R) <= N: R -= 1
if L <= R and L <= k_max and R >= 1:
L2 = max(1, L)
R2 = min(k_max, R)
if L2 <= R2:
forbidden = count_with_residue(R2) - count_with_residue(L2 - 1)
base -= forbidden
total += base
disc_plus = 13 * q2 + 12 * N
s_plus = isqrt(disc_plus)
# Region C
k_max_c = (s_plus - 13 * q) // 6
if k_max_c > 0:
while k_max_c > 0 and zC(k_max_c) > N:
k_max_c -= 1
while zC(k_max_c + 1) <= N:
k_max_c += 1
if k_max_c > 0:
total += count_exclude13(k_max_c)
# Region D
k_min_d = int(math.floor(self.beta * q)) + 1
if k_min_d < 1: k_min_d = 1
while True:
x_val = q2 - 4 * q * k_min_d + k_min_d * k_min_d
if x_val > 0: break
k_min_d += 1
k_max_d = (s_plus + q) // 6
while k_max_d > 0 and zD(k_max_d) > N:
k_max_d -= 1
while zD(k_max_d + 1) <= N:
k_max_d += 1
if k_min_d <= k_max_d:
total += count_with_residue(k_max_d) - count_with_residue(k_min_d - 1)
if N >= 3:
total += 1
return total
def solve():
N = 100000000000000
q_limit = isqrt(N)
c = Counter(q_limit)
ans = c.compute(N)
return str(ans)
if __name__ == "__main__":
print(solve())
Java
public class Euler769 {
static long isqrt(long x) {
if (x < 0)
throw new IllegalArgumentException();
if (x == 0)
return 0;
long r = (long) Math.sqrt((double) x);
long y = r;
while ((y + 1) * (y + 1) <= x)
++y;
while (y * y > x)
--y;
return y;
}
static int[] buildSpf(int limit) {
int[] spf = new int[limit + 1];
if (limit >= 1)
spf[1] = 1;
for (int i = 2; i <= limit; ++i) {
if (spf[i] == 0) {
spf[i] = i;
if ((long) i * i <= limit) {
for (int j = i * i; j <= limit; j += i) {
if (spf[j] == 0)
spf[j] = i;
}
}
}
}
return spf;
}
static class Counter {
int[] spf;
double alpha;
double beta;
int[] inv13 = new int[13];
Counter(int[] spfRef) {
spf = spfRef;
alpha = 2.0 - Math.sqrt(3.0);
beta = 2.0 + Math.sqrt(3.0);
for (int i = 1; i < 13; ++i) {
for (int j = 1; j < 13; ++j) {
if ((i * j) % 13 == 1) {
inv13[i] = j;
break;
}
}
}
}
void buildSquarefreeDivs(int q, int[] divs, int[] mus, int[] countOut) {
divs[0] = 1;
mus[0] = 1;
int count = 1;
int x = q;
while (x > 1) {
int p = spf[x];
while (x % p == 0)
x /= p;
int current = count;
for (int i = 0; i < current; ++i) {
divs[count] = divs[i] * p;
mus[count] = -mus[i];
count++;
}
}
countOut[0] = count;
}
long coprimeCount(int[] divs, int[] mus, int count, long kMax) {
if (kMax <= 0)
return 0;
long total = 0;
for (int i = 0; i < count; ++i) {
total += mus[i] * (kMax / divs[i]);
}
return total;
}
long coprimeCountResidue(int[] divs, int[] mus, int count, long kMax, int residue) {
if (kMax <= 0)
return 0;
long total = 0;
for (int i = 0; i < count; ++i) {
int d = divs[i];
int dMod = d % 13;
int inv = inv13[dMod];
int t0 = (residue * inv) % 13;
if (t0 == 0)
continue;
long first = (long) d * t0;
if (first > kMax)
continue;
total += mus[i] * (1 + (kMax - first) / (13L * d));
}
return total;
}
long countWithResidue(int[] divs, int[] mus, int count, long kMax, int qMod13, int residue) {
if (kMax <= 0)
return 0;
long base = coprimeCount(divs, mus, count, kMax);
if (qMod13 != 0) {
base -= coprimeCountResidue(divs, mus, count, kMax, residue);
}
return base;
}
long countExclude13(int[] divs, int[] mus, int count, long kMax, int qMod13) {
if (kMax <= 0)
return 0;
long base = coprimeCount(divs, mus, count, kMax);
if (qMod13 != 0) {
base -= coprimeCount(divs, mus, count, kMax / 13);
}
return base;
}
long compute(long N) {
long qLimit = isqrt(N);
if (qLimit == 0)
return 0;
long total = 0;
int[] divs = new int[256];
int[] mus = new int[256];
int[] countOut = new int[1];
for (long q = 1; q <= qLimit; ++q) {
buildSquarefreeDivs((int) q, divs, mus, countOut);
int cnt = countOut[0];
long q2 = q * q;
int qMod13 = (int) (q % 13);
int residue = qMod13 == 0 ? 0 : (13 - (2 * qMod13) % 13) % 13;
// Region A
long kMax = (long) Math.floor(alpha * q);
if (kMax > 0) {
while (kMax > 0) {
long xVal = q2 - 4 * q * kMax + kMax * kMax;
if (xVal > 0)
break;
kMax--;
}
if (kMax > 0) {
long base = countWithResidue(divs, mus, cnt, kMax, qMod13, residue);
long disc = 13 * q2 - 12 * N;
if (disc > 0) {
long s = isqrt(disc);
long L = (q - s + 5) / 6;
long R = (q + s) / 6;
while (L <= kMax && (q2 + q * L - 3 * L * L) <= N)
L++;
while (R >= 1 && (q2 + q * R - 3 * R * R) <= N)
R--;
if (L <= R && L <= kMax && R >= 1) {
long L2 = Math.max(1, L);
long R2 = Math.min(kMax, R);
if (L2 <= R2) {
long forbidden = countWithResidue(divs, mus, cnt, R2, qMod13, residue) -
countWithResidue(divs, mus, cnt, L2 - 1, qMod13, residue);
base -= forbidden;
}
}
}
total += base;
}
}
// discPlus
// Note: 13 * q^2 + 12 * N can be up to 13*10^14 + 12*10^14 = 2.5 * 10^15 (fits
// in long)
long discPlus = 13 * q2 + 12 * N;
long sPlus = isqrt(discPlus);
// Region C
long kMaxC = (sPlus - 13 * q) / 6;
if (kMaxC > 0) {
while (kMaxC > 0 && (13 * q2 + 13 * q * kMaxC + 3 * kMaxC * kMaxC) > N)
kMaxC--;
while ((13 * q2 + 13 * q * (kMaxC + 1) + 3 * (kMaxC + 1) * (kMaxC + 1)) <= N)
kMaxC++;
if (kMaxC > 0) {
total += countExclude13(divs, mus, cnt, kMaxC, qMod13);
}
}
// Region D
long kMinD = (long) Math.floor(beta * q) + 1;
if (kMinD < 1)
kMinD = 1;
while (true) {
long xVal = q2 - 4 * q * kMinD + kMinD * kMinD;
if (xVal > 0)
break;
kMinD++;
}
long kMaxD = (sPlus + q) / 6;
while (kMaxD > 0 && (3 * kMaxD * kMaxD - q * kMaxD - q2) > N)
kMaxD--;
while ((3 * (kMaxD + 1) * (kMaxD + 1) - q * (kMaxD + 1) - q2) <= N)
kMaxD++;
if (kMinD <= kMaxD) {
total += countWithResidue(divs, mus, cnt, kMaxD, qMod13, residue) -
countWithResidue(divs, mus, cnt, kMinD - 1, qMod13, residue);
}
}
if (N >= 3) {
total += 1;
}
return total;
}
}
public static String solve() {
long n = 100000000000000L;
int qLimit = (int) isqrt(n);
int[] spf = buildSpf(qLimit);
Counter c = new Counter(spf);
long ans = c.compute(n);
return Long.toString(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}