Problem 962: Angular Bisector and Tangent 2
View on Project EulerProject Euler Problem 962 Solution
EulerSolve provides an optimized solution for Project Euler Problem 962, Angular Bisector and Tangent 2, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Problem 962 asks for all integer triples \((a,b,c)\) with \(a \le b \le c \lt a+b\) and \(a+b+c \le N\) for which the geometric construction in the statement leads to \[ Q=\frac{a^3(a+b-c)(a+b+c)}{b(a+b)^2} \] being an integer perfect square. For the Project Euler instance, \(N=10^6\). A direct sweep over all triangles is far too slow, so the solution replaces the \((a,b,c)\)-search by a number-theoretic count over a much smaller parameter space. Mathematical Approach The key step is to separate the common scale of \(a\) and \(b\), prove that every valid triangle has a rigid difference-of-squares form, and then count the remaining possibilities through squarefree data and interval arithmetic. Separating the common scale of \(a\) and \(b\) Let \[ g=\gcd(a,b), \qquad a=ug, \qquad b=vg, \] with \(\gcd(u,v)=1\) and \(u \le v\). Write \[ k=u+v. \] Then \(a+b=kg\), so the square test becomes \[ Q=\frac{u^3\bigl(k^2g^2-c^2\bigr)}{vk^2}. \] This already shows why the reduced ratio \(u:v\) is the real state space: once \(u\) and \(v\) are fixed, the remaining freedom sits inside \(g\) and \(c\). Why \(c\) must be a multiple of \(k\) If \(Q\) is an integer, then \(k^2\) divides the numerator \(u^3(k^2g^2-c^2)\). Because \(\gcd(u,k)=1\), the factor \(u^3\) is irrelevant for divisibility by \(k^2\), so \[ k^2 \mid \bigl(k^2g^2-c^2\bigr)....
Detailed mathematical approach
Problem Summary
Problem 962 asks for all integer triples \((a,b,c)\) with \(a \le b \le c \lt a+b\) and \(a+b+c \le N\) for which the geometric construction in the statement leads to
\[ Q=\frac{a^3(a+b-c)(a+b+c)}{b(a+b)^2} \]
being an integer perfect square. For the Project Euler instance, \(N=10^6\). A direct sweep over all triangles is far too slow, so the solution replaces the \((a,b,c)\)-search by a number-theoretic count over a much smaller parameter space.
Mathematical Approach
The key step is to separate the common scale of \(a\) and \(b\), prove that every valid triangle has a rigid difference-of-squares form, and then count the remaining possibilities through squarefree data and interval arithmetic.
Separating the common scale of \(a\) and \(b\)
Let
\[ g=\gcd(a,b), \qquad a=ug, \qquad b=vg, \]
with \(\gcd(u,v)=1\) and \(u \le v\). Write
\[ k=u+v. \]
Then \(a+b=kg\), so the square test becomes
\[ Q=\frac{u^3\bigl(k^2g^2-c^2\bigr)}{vk^2}. \]
This already shows why the reduced ratio \(u:v\) is the real state space: once \(u\) and \(v\) are fixed, the remaining freedom sits inside \(g\) and \(c\).
Why \(c\) must be a multiple of \(k\)
If \(Q\) is an integer, then \(k^2\) divides the numerator \(u^3(k^2g^2-c^2)\). Because \(\gcd(u,k)=1\), the factor \(u^3\) is irrelevant for divisibility by \(k^2\), so
\[ k^2 \mid \bigl(k^2g^2-c^2\bigr). \]
The term \(k^2g^2\) is already a multiple of \(k^2\), hence \(k^2 \mid c^2\), which implies
\[ k \mid c. \]
So every valid triangle can be written as
\[ c=kt \]
for some integer \(t\) with \(0 \lt t \lt g\). Substituting this into the previous formula collapses the \(k^2\)-denominator:
\[ Q=\frac{u^3(g^2-t^2)}{v}. \]
From the triangle to the pair \((r,s)\)
Now factor the difference of squares by setting
\[ r=g-t, \qquad s=g+t. \]
Then
\[ g=\frac{r+s}{2}, \qquad t=\frac{s-r}{2}, \]
so \(r\) and \(s\) are positive integers with \(r \lt s\) and
\[ r \equiv s \pmod 2. \]
The triangle sides become
\[ a=\frac{u(r+s)}{2}, \qquad b=\frac{v(r+s)}{2}, \qquad c=\frac{k(s-r)}{2}, \]
and the perimeter is especially simple:
\[ a+b+c=ks. \]
The square test is now
\[ Q=\frac{u^3rs}{v}. \]
This is the form used by the optimized algorithm. Instead of iterating over \(c\), it counts admissible pairs \((r,s)\) for each reduced ratio \(u:v\).
Turning the square condition into squarefree data
Write
\[ u=d\,h^2, \]
where \(d=\operatorname{sf}(u)\) is the squarefree part of \(u\). Then
\[ u^3=d^3h^6=d^2(dh^3)^2, \]
so \(u^3\) contributes one squarefree copy of \(d\) and everything else is already a square. Therefore
\[ Q \text{ is a perfect square } \Longleftrightarrow \frac{drs}{v} \text{ is a perfect square.} \]
That is why the implementations precompute the squarefree kernel of \(u\). For primes dividing \(v\), the product \(rs\) must supply enough exponent to cancel the denominator and leave an even valuation. For primes dividing \(d\), the product \(rs\) must carry an odd valuation. These prime-by-prime requirements are compressed into a smallest admissible base value \(s_0\), and every valid \(s\) for fixed \((u,v,r)\) has the form
\[ s=s_0n^2 \]
with \(n \ge 1\).
Bounds coming from the triangle ordering and the perimeter cap
Let
\[ L=\left\lfloor\frac{N}{k}\right\rfloor. \]
Since the perimeter equals \(ks\), we must have \(s \le L\). The remaining geometric restriction is \(c \ge b\), which becomes
\[ \frac{k(s-r)}{2} \ge \frac{v(r+s)}{2}. \]
After simplifying with \(k=u+v\), this becomes
\[ us \ge (u+2v)r. \]
So for fixed \(u,v,r\), the first feasible value of \(s\) is
\[ s_{\min}=\left\lceil\frac{(u+2v)r}{u}\right\rceil. \]
Because \(s=s_0n^2\), the allowed range for \(n\) is
\[ n_{\min}=\left\lceil\sqrt{\left\lceil\frac{s_{\min}}{s_0}\right\rceil}\right\rceil, \qquad n_{\max}=\left\lfloor\sqrt{\frac{L}{s_0}}\right\rfloor. \]
The parity condition \(r \equiv s \pmod 2\) now becomes a filter on \(n\). If \(s_0\) is even, then every admissible \(s\) is even, so only even \(r\) can contribute. If \(s_0\) is odd, then \(s \equiv n^2 \equiv n \pmod 2\), so we need
\[ n \equiv r \pmod 2. \]
Each branch therefore contributes only the number of integers \(n\) in a short interval with the correct parity.
Why the outer loop stops at \(k \le \lfloor\sqrt[3]{2N^2}\rfloor+2\)
The solutions use a sharp global bound on \(k\). Since \(v \ge k/2\) and \(d \ge 1\), any valid branch must satisfy
\[ vd \le rs \le s^2 \le L^2 \le \left(\frac{N}{k}\right)^2. \]
Therefore
\[ \frac{k}{2} \le \left(\frac{N}{k}\right)^2, \qquad\text{so}\qquad k^3 \le 2N^2. \]
This is why only \(k\) up to roughly \((2N^2)^{1/3}\) need to be examined.
Worked example
Consider the valid triangle \((a,b,c)=(12,15,18)\). Then
\[ g=\gcd(12,15)=3,\qquad u=4,\qquad v=5,\qquad k=9. \]
Since \(c=18=9\cdot 2\), we have \(t=2\). Hence
\[ r=g-t=1,\qquad s=g+t=5. \]
Recovering the sides from \((u,v,r,s)\) gives
\[ a=\frac{4(1+5)}{2}=12,\qquad b=\frac{5(1+5)}{2}=15,\qquad c=\frac{9(5-1)}{2}=18. \]
The square test becomes
\[ Q=\frac{4^3 \cdot 1 \cdot 5}{5}=64=8^2. \]
Here \(\operatorname{sf}(4)=1\), so the square condition reduces to \(s/5\) being a square when \(r=1\). The smallest choice is \(s_0=5\), and \(s=5=s_0\cdot 1^2\) is indeed the first admissible value above
\[ s_{\min}=\left\lceil\frac{(4+2\cdot 5)\cdot 1}{4}\right\rceil=4. \]
This single example mirrors the whole algorithm: fix a reduced ratio \(u:v\), turn the triangle into \(r\) and \(s\), build the minimal base \(s_0\), then count square multiples inside the allowed interval.
How the Code Works
Prime tables and squarefree kernels
The C++, Python, and Java implementations begin by computing a smallest-prime-factor table. From it they derive two reusable data sets: the squarefree part \(\operatorname{sf}(u)\) for every possible \(u\), and a compact prime-exponent factorization for every integer that can appear as \(v\), \(r\), or \(\operatorname{sf}(u)\). This removes repeated factorization work from the main count.
Processing one reduced ratio
For each \(k\), the implementation iterates over \(v\) from \(\lceil k/2 \rceil\) to \(k-1\), keeps only \(\gcd(v,k)=1\), and sets \(u=k-v\). It then computes \(L=\lfloor N/k \rfloor\) and discards the branch immediately when \(v\operatorname{sf}(u) \gt L^2\), because then even the largest possible product \(rs\) cannot satisfy the square condition.
Next it loops over
\[ 1 \le r \le \left\lfloor\frac{uL}{u+2v}\right\rfloor. \]
For each such \(r\), the implementation merges the prime data of \(r\), \(v\), and \(\operatorname{sf}(u)\) to construct the smallest \(s_0\) that can make \(drs/v\) a square. If \(s_0 \gt L\), the branch is impossible. Otherwise the code converts the conditions into the interval \([n_{\min},n_{\max}]\), applies the parity rule, and adds the number of admissible \(n\).
Final accumulation
The total answer is the sum of these contributions over all \(k\). The Python implementation evaluates the outer \(k\)-loop serially. The C++ and Java implementations use parallel execution for that outer sweep, because each \(k\) contributes independently once the factor tables have been built.
Complexity Analysis
The outer bound \(k \le O(N^{2/3})\) already cuts the search space drastically. For a fixed \(k\), there are \(O(k)\) candidates for \(v\), and for each surviving \((u,v)\) the parameter \(r\) ranges up to \(O(N/k)\). A coarse upper bound for the arithmetic branches is therefore
\[ \sum_{k \le O(N^{2/3})} O(k)\,O\!\left(\frac{N}{k}\right)=O(N^{5/3}). \]
Each branch only performs a short merge of prime-exponent lists together with constant-time integer-root and interval calculations, so the practical runtime is much better than a triangle-by-triangle scan. The early coprimality filter and the \(v\operatorname{sf}(u) \le L^2\) pruning remove many branches before any expensive work happens.
Memory usage is dominated by the prime tables: the smallest-prime-factor sieve, the squarefree-part array, and the factor table up to the largest needed integer. This is essentially linear in the sieve bound, with only a small logarithmic factor from storing prime exponents.
Footnotes and References
- Problem page: https://projecteuler.net/problem=962
- Angle bisector theorem: Wikipedia - Angle bisector theorem
- Greatest common divisor: Wikipedia - Greatest common divisor
- Difference of two squares: Wikipedia - Difference of two squares
- Square-free integer: Wikipedia - Square-free integer
- Perfect square: Wikipedia - Perfect square
- Triangle inequality: Wikipedia - Triangle inequality
Problem 962 source code
C++
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <thread>
#include <vector>
namespace {
using u64 = std::uint64_t;
using i64 = std::int64_t;
struct Options {
int limit = 1000000;
unsigned threads = std::thread::hardware_concurrency();
bool run_checkpoints = true;
};
struct PrimeExp {
int p;
int e;
};
bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
int parsed = 0;
for (char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10 + static_cast<int>(c - '0');
}
value = parsed;
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg, const std::string& prefix,
unsigned& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
unsigned parsed = 0;
for (char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10U + static_cast<unsigned>(c - '0');
}
value = parsed;
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
continue;
}
if (parse_unsigned_after_prefix(arg, "--threads=", options.threads)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.limit < 3) {
return false;
}
if (options.threads == 0) {
options.threads = 1;
}
return true;
}
int iroot_floor(i64 x) {
if (x <= 0) {
return 0;
}
i64 r = static_cast<i64>(std::sqrt(static_cast<long double>(x)));
while ((r + 1) * (r + 1) <= x) {
++r;
}
while (r * r > x) {
--r;
}
return static_cast<int>(r);
}
int iroot_ceil(i64 x) {
if (x <= 0) {
return 0;
}
int r = iroot_floor(x);
if (static_cast<i64>(r) * r < x) {
++r;
}
return r;
}
int cube_root_floor(i64 x) {
int lo = 0;
int hi = static_cast<int>(std::cbrt(static_cast<long double>(x))) + 5;
while (static_cast<i64>(hi) * hi * hi <= x) {
++hi;
}
while (lo + 1 < hi) {
const int mid = lo + (hi - lo) / 2;
const i64 cube = static_cast<i64>(mid) * mid * mid;
if (cube <= x) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
std::vector<int> build_spf(int n) {
std::vector<int> spf(n + 1);
for (int i = 0; i <= n; ++i) {
spf[i] = i;
}
if (n >= 0) {
spf[0] = 0;
}
if (n >= 1) {
spf[1] = 1;
}
for (int i = 2; static_cast<i64>(i) * i <= n; ++i) {
if (spf[i] != i) {
continue;
}
for (int j = i * i; j <= n; j += i) {
if (spf[j] == j) {
spf[j] = i;
}
}
}
return spf;
}
std::vector<int> build_squarefree_parts(int n, const std::vector<int>& spf) {
std::vector<int> sf(n + 1, 1);
sf[0] = 0;
for (int x = 1; x <= n; ++x) {
int t = x;
int value = 1;
while (t > 1) {
const int p = spf[t];
int c = 0;
while (t % p == 0) {
t /= p;
c ^= 1;
}
if (c != 0) {
value *= p;
}
}
sf[x] = value;
}
return sf;
}
std::vector<std::vector<PrimeExp>> build_factor_table(int n, const std::vector<int>& spf) {
std::vector<std::vector<PrimeExp>> table(n + 1);
for (int x = 2; x <= n; ++x) {
int t = x;
while (t > 1) {
const int p = spf[t];
int e = 0;
while (t % p == 0) {
t /= p;
++e;
}
table[x].push_back({p, e});
}
}
return table;
}
bool multiply_power_limited(i64& value, int p, int e, int limit) {
for (int i = 0; i < e; ++i) {
value *= p;
if (value > limit) {
return false;
}
}
return true;
}
u64 count_for_k(int limit, int k, const std::vector<int>& squarefree_u,
const std::vector<std::vector<PrimeExp>>& factors) {
const int L = limit / k;
if (L < 2) {
return 0;
}
u64 subtotal = 0;
const int v_min = (k + 1) / 2;
for (int v = v_min; v < k; ++v) {
if (std::gcd(v, k) != 1) {
continue;
}
const int u = k - v;
const int sf_u = squarefree_u[u];
const i64 W = static_cast<i64>(v) * sf_u;
if (W > static_cast<i64>(L) * L) {
continue;
}
std::vector<PrimeExp> wf = factors[v];
if (sf_u != 1) {
for (const PrimeExp pe : factors[sf_u]) {
wf.push_back({pe.p, 1});
}
}
std::sort(wf.begin(), wf.end(), [](const PrimeExp& a, const PrimeExp& b) {
return a.p < b.p;
});
const int r_max = static_cast<int>((static_cast<i64>(u) * L) / (u + 2 * v));
for (int r = 1; r <= r_max; ++r) {
const std::vector<PrimeExp>& rf = factors[r];
i64 s0 = 1;
std::size_t i = 0;
std::size_t j = 0;
bool valid = true;
while (i < wf.size() || j < rf.size()) {
int p = 0;
int e = 0;
int a = 0;
if (j == rf.size() || (i < wf.size() && wf[i].p < rf[j].p)) {
p = wf[i].p;
e = wf[i].e;
a = 0;
++i;
} else if (i == wf.size() || rf[j].p < wf[i].p) {
p = rf[j].p;
e = 0;
a = rf[j].e;
++j;
} else {
p = wf[i].p;
e = wf[i].e;
a = rf[j].e;
++i;
++j;
}
int b0 = 0;
if (a < e) {
b0 = e - a;
} else {
b0 = (a - e) & 1;
}
if (b0 > 0) {
if (!multiply_power_limited(s0, p, b0, L)) {
valid = false;
break;
}
}
}
if (!valid) {
continue;
}
const int s_min = static_cast<int>(
(static_cast<i64>(u + 2 * v) * r + (u - 1)) / u);
if (s_min > L) {
continue;
}
const i64 ratio_min = (static_cast<i64>(s_min) + s0 - 1) / s0;
int n_min = iroot_ceil(ratio_min);
const int n_max = iroot_floor(L / s0);
if (n_min > n_max) {
continue;
}
if ((s0 & 1LL) == 0) {
if ((r & 1) == 0) {
subtotal += static_cast<u64>(n_max - n_min + 1);
}
continue;
}
const int want_parity = r & 1;
if ((n_min & 1) != want_parity) {
++n_min;
}
if (n_min <= n_max) {
subtotal += static_cast<u64>((n_max - n_min) / 2 + 1);
}
}
}
return subtotal;
}
u64 solve_fast(int limit, unsigned threads) {
const i64 n2_twice = 2LL * limit * limit;
const int k_max = cube_root_floor(n2_twice) + 2;
const int max_needed = std::max(k_max, limit / 6 + 5);
const std::vector<int> spf = build_spf(max_needed);
const std::vector<int> squarefree_u = build_squarefree_parts(k_max, spf);
const std::vector<std::vector<PrimeExp>> factors = build_factor_table(max_needed, spf);
if (threads <= 1 || k_max < 128) {
u64 total = 0;
for (int k = 2; k <= k_max; ++k) {
total += count_for_k(limit, k, squarefree_u, factors);
}
return total;
}
const unsigned use_threads = std::min<unsigned>(threads, static_cast<unsigned>(k_max - 1));
std::atomic<int> next_k(2);
std::vector<std::thread> pool;
std::vector<u64> partial(use_threads, 0);
pool.reserve(use_threads);
for (unsigned t = 0; t < use_threads; ++t) {
pool.emplace_back([&, t]() {
u64 local = 0;
while (true) {
const int k = next_k.fetch_add(1, std::memory_order_relaxed);
if (k > k_max) {
break;
}
local += count_for_k(limit, k, squarefree_u, factors);
}
partial[t] = local;
});
}
for (auto& th : pool) {
th.join();
}
u64 total = 0;
for (u64 part : partial) {
total += part;
}
return total;
}
u64 solve_bruteforce(int limit) {
u64 total = 0;
for (int a = 1; a <= limit / 3; ++a) {
for (int b = a; b <= (limit - a) / 2; ++b) {
const int s = a + b;
const int c_max = std::min(s - 1, limit - s);
for (int c = b; c <= c_max; ++c) {
const i64 num = static_cast<i64>(a) * a * a * (s - c) * (s + c);
const i64 den = static_cast<i64>(b) * s * s;
if (num % den != 0) {
continue;
}
const i64 q = num / den;
const i64 rt = static_cast<i64>(std::sqrt(static_cast<long double>(q)));
if (rt * rt == q || (rt + 1) * (rt + 1) == q) {
++total;
}
}
}
}
return total;
}
bool run_checkpoints() {
for (int limit : {50, 80, 120}) {
const u64 fast = solve_fast(limit, 1);
const u64 brute = solve_bruteforce(limit);
if (fast != brute) {
std::cerr << "Bruteforce checkpoint failed for limit=" << limit << ": fast=" << fast
<< ", brute=" << brute << '\n';
return false;
}
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 2;
}
const u64 answer = solve_fast(options.limit, options.threads);
std::cout << answer << '\n';
return 0;
}
Python
import math
def solve():
limit = 1000000
def isqrt_f(x):
if x <= 0: return 0
r = int(math.isqrt(x))
while (r+1)*(r+1)<=x: r+=1
while r*r>x: r-=1
return r
def isqrt_c(x):
if x<=0: return 0
r = isqrt_f(x)
if r*r < x: r += 1
return r
def cbrt_f(x):
lo = 0; hi = int(x**(1/3))+5
while hi*hi*hi <= x: hi += 1
while lo+1 < hi:
mid = (lo+hi)//2
if mid*mid*mid <= x: lo = mid
else: hi = mid
return lo
# SPF sieve
n2 = 2*limit*limit; k_max = cbrt_f(n2)+2
mx = max(k_max, limit//6+5)
spf = list(range(mx+1)); spf[0] = 0; spf[1] = 1
for i in range(2, int(mx**0.5)+1):
if spf[i] != i: continue
for j in range(i*i, mx+1, i):
if spf[j] == j: spf[j] = i
# Squarefree parts
sf = [1]*(k_max+1); sf[0] = 0
for x in range(1, k_max+1):
t = x; val = 1
while t > 1:
p = spf[t]; c = 0
while t%p==0: t//=p; c ^= 1
if c: val *= p
sf[x] = val
# Factor tables
def factorize(x):
fs = []; t = x
while t > 1:
p = spf[t]; e = 0
while t%p==0: t//=p; e += 1
fs.append((p, e))
return fs
factors = [[] for _ in range(mx+1)]
for x in range(2, mx+1): factors[x] = factorize(x)
def count_for_k(k):
L = limit // k
if L < 2: return 0
sub = 0; v_min = (k+1)//2
for v in range(v_min, k):
if math.gcd(v, k) != 1: continue
u = k - v; sf_u = sf[u]; W = v * sf_u
if W > L*L: continue
wf = list(factors[v])
if sf_u != 1: wf += [(pe[0], 1) for pe in factors[sf_u]]
wf.sort()
r_max = (u * L) // (u + 2*v)
for r in range(1, r_max+1):
rf = factors[r]; s0 = 1; i = j = 0; ok = True
while i < len(wf) or j < len(rf):
if j == len(rf) or (i < len(wf) and wf[i][0] < rf[j][0]):
p, e = wf[i]; a = 0; i += 1
elif i == len(wf) or rf[j][0] < wf[i][0]:
p, a = rf[j]; e = 0; j += 1
else:
p, e = wf[i]; a = rf[j][1]; i += 1; j += 1
b0 = e-a if a < e else (a-e)&1
if b0 > 0:
s0 *= p**b0
if s0 > L: ok = False; break
if not ok: continue
s_min = ((u+2*v)*r + u-1)//u
if s_min > L: continue
ratio_min = (s_min+s0-1)//s0
n_min = isqrt_c(ratio_min); n_max = isqrt_f(L//s0)
if n_min > n_max: continue
if s0 & 1 == 0:
if r & 1 == 0: sub += n_max-n_min+1
continue
wp = r & 1
if (n_min & 1) != wp: n_min += 1
if n_min <= n_max: sub += (n_max-n_min)//2+1
return sub
total = sum(count_for_k(k) for k in range(2, k_max+1))
return str(total)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.stream.IntStream;
public class Euler962 {
static int irootFloor(long x) {
if (x <= 0)
return 0;
long r = (long) Math.sqrt(x);
while ((r + 1) * (r + 1) <= x) {
++r;
}
while (r * r > x) {
--r;
}
return (int) r;
}
static int irootCeil(long x) {
if (x <= 0)
return 0;
int r = irootFloor(x);
if ((long) r * r < x) {
++r;
}
return r;
}
static int cubeRootFloor(long x) {
int lo = 0;
int hi = (int) Math.cbrt(x) + 5;
while ((long) hi * hi * hi <= x) {
++hi;
}
while (lo + 1 < hi) {
int mid = lo + (hi - lo) / 2;
long cube = (long) mid * mid * mid;
if (cube <= x) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
static int[] buildSpf(int n) {
int[] spf = new int[n + 1];
for (int i = 0; i <= n; ++i) {
spf[i] = i;
}
for (int i = 2; (long) i * i <= n; ++i) {
if (spf[i] != i)
continue;
for (int j = i * i; j <= n; j += i) {
if (spf[j] == j) {
spf[j] = i;
}
}
}
return spf;
}
static int[] buildSquarefreeParts(int n, int[] spf) {
int[] sf = new int[n + 1];
Arrays.fill(sf, 1);
sf[0] = 0;
for (int x = 1; x <= n; ++x) {
int t = x;
int value = 1;
while (t > 1) {
int p = spf[t];
int c = 0;
while (t % p == 0) {
t /= p;
c ^= 1;
}
if (c != 0) {
value *= p;
}
}
sf[x] = value;
}
return sf;
}
static class PrimeExp implements Comparable<PrimeExp> {
int p;
int e;
PrimeExp(int p, int e) {
this.p = p;
this.e = e;
}
@Override
public int compareTo(PrimeExp o) {
return Integer.compare(this.p, o.p);
}
}
static List<List<PrimeExp>> buildFactorTable(int n, int[] spf) {
List<List<PrimeExp>> table = new ArrayList<>(n + 1);
for (int i = 0; i <= n; i++) {
table.add(new ArrayList<>());
}
for (int x = 2; x <= n; ++x) {
int t = x;
while (t > 1) {
int p = spf[t];
int e = 0;
while (t % p == 0) {
t /= p;
++e;
}
table.get(x).add(new PrimeExp(p, e));
}
}
return table;
}
static long gcd(long a, long b) {
while (b != 0) {
long t = b;
b = a % b;
a = t;
}
return a;
}
static long countForK(int limit, int k, int[] squarefreeU, List<List<PrimeExp>> factors) {
int L = limit / k;
if (L < 2)
return 0;
long subtotal = 0;
int vMin = (k + 1) / 2;
for (int v = vMin; v < k; ++v) {
if (gcd(v, k) != 1)
continue;
int u = k - v;
int sfU = squarefreeU[u];
long W = (long) v * sfU;
if (W > (long) L * L)
continue;
List<PrimeExp> wf = new ArrayList<>(factors.get(v));
if (sfU != 1) {
for (PrimeExp pe : factors.get(sfU)) {
boolean found = false;
for (PrimeExp existing : wf) {
if (existing.p == pe.p) {
existing.e += 1;
found = true;
break;
}
}
if (!found) {
wf.add(new PrimeExp(pe.p, 1));
}
}
}
Collections.sort(wf);
int rMax = (int) (((long) u * L) / (u + 2 * v));
for (int r = 1; r <= rMax; ++r) {
List<PrimeExp> rf = factors.get(r);
long s0 = 1;
int i = 0;
int j = 0;
boolean valid = true;
while (i < wf.size() || j < rf.size()) {
int p = 0, e = 0, a = 0;
if (j == rf.size() || (i < wf.size() && wf.get(i).p < rf.get(j).p)) {
p = wf.get(i).p;
e = wf.get(i).e;
a = 0;
++i;
} else if (i == wf.size() || rf.get(j).p < wf.get(i).p) {
p = rf.get(j).p;
e = 0;
a = rf.get(j).e;
++j;
} else {
p = wf.get(i).p;
e = wf.get(i).e;
a = rf.get(j).e;
++i;
++j;
}
int b0 = (a < e) ? (e - a) : ((a - e) & 1);
if (b0 > 0) {
for (int k_idx = 0; k_idx < b0; ++k_idx) {
s0 *= p;
if (s0 > L) {
valid = false;
break;
}
}
if (!valid)
break;
}
}
if (!valid)
continue;
int sMin = (int) (((long) (u + 2 * v) * r + (u - 1)) / u);
if (sMin > L)
continue;
long ratioMin = (sMin + s0 - 1) / s0;
int nMin = irootCeil(ratioMin);
int nMax = irootFloor(L / s0);
if (nMin > nMax)
continue;
if ((s0 & 1L) == 0) {
if ((r & 1) == 0) {
subtotal += (nMax - nMin + 1);
}
continue;
}
int wantParity = r & 1;
if ((nMin & 1) != wantParity) {
++nMin;
}
if (nMin <= nMax) {
subtotal += (nMax - nMin) / 2 + 1;
}
}
}
return subtotal;
}
static long solveFast(int limit) {
long n2Twice = 2L * limit * limit;
int kMax = cubeRootFloor(n2Twice) + 2;
int maxNeeded = Math.max(kMax, limit / 6 + 5);
int[] spf = buildSpf(maxNeeded);
int[] squarefreeU = buildSquarefreeParts(kMax, spf);
List<List<PrimeExp>> factors = buildFactorTable(maxNeeded, spf);
return IntStream.rangeClosed(2, kMax)
.parallel()
.mapToLong(k -> countForK(limit, k, squarefreeU, factors))
.sum();
}
public static String solve() {
return Long.toString(solveFast(1000000));
}
public static void main(String[] args) {
System.out.println(solve());
}
}