Problem 450: Hypocycloid and Lattice Points

Project Euler Problem 450 Solution

EulerSolve provides an optimized solution for Project Euler Problem 450, Hypocycloid and Lattice Points, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For integers $R \gt r \gt 0$, the hypocycloid is defined by $$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$ We only keep parameters $t$ for which $\sin t$ and $\cos t$ are rational, and among those points we keep only the ones with integer coordinates. For a fixed pair $(R,r)$, let $S(R,r)$ be the sum of $|x|+|y|$ over the distinct lattice points obtained this way. The global target is $$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$ The implementations verify the checkpoints $S(3,1)=10$, $T(10)=524$, $T(100)=580442$, and $T(1000)=583108600$. Mathematical Approach Step 1: Separate the Common Scale Write $$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$ The condition $r\lt R/2$ becomes $q \lt p/2$, so $e \gt q$. Now set $t=q\theta$ and define $z=e^{i\theta}$. Then $z^q=e^{it}$, and the hypocycloid coordinates combine into one complex expression: $$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$ This identity is the key simplification. For each reduced pair $(p,q)$, the geometry is fixed, and changing $g$ only scales the final lattice point linearly....

Detailed mathematical approach

Problem Summary

For integers $R \gt r \gt 0$, the hypocycloid is defined by

$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$

We only keep parameters $t$ for which $\sin t$ and $\cos t$ are rational, and among those points we keep only the ones with integer coordinates. For a fixed pair $(R,r)$, let $S(R,r)$ be the sum of $|x|+|y|$ over the distinct lattice points obtained this way. The global target is

$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$

The implementations verify the checkpoints $S(3,1)=10$, $T(10)=524$, $T(100)=580442$, and $T(1000)=583108600$.

Mathematical Approach

Step 1: Separate the Common Scale

Write

$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$

The condition $r\lt R/2$ becomes $q \lt p/2$, so $e \gt q$. Now set $t=q\theta$ and define $z=e^{i\theta}$. Then $z^q=e^{it}$, and the hypocycloid coordinates combine into one complex expression:

$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$

This identity is the key simplification. For each reduced pair $(p,q)$, the geometry is fixed, and changing $g$ only scales the final lattice point linearly.

Step 2: Rational Points on the Unit Circle

Every rational point on the unit circle can be written as

$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$

The trivial case $D=1$ gives the four roots of unity $(\pm1,0)$ and $(0,\pm1)$. Every nontrivial rational point comes from a primitive Pythagorean triple, together with sign changes and swapping the two legs.

Substituting this parametrization into the complex form gives

$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$

If we write

$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e,$$

then the entire arithmetic question becomes: after cancelling the common factors of $D^e$, when does the remaining denominator divide $g$?

Step 3: Exact Integrality Criterion

Let

$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$

After reducing the fraction, the point becomes

$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$

Therefore the coordinates are integers if and only if

$$\delta \mid g.$$

If $g=m\delta$, then this one geometric shape contributes

$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$

Summing over all admissible multiples of $\delta$ produces a triangular number:

$$\sum_{m=1}^{m_{\max}} m=\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$

This is why the global computation can be organized by reduced shapes instead of looping over all $(R,r)$ one by one.

Step 4: The $D=1$ Case Produces the Main Term

When $D=1$, the only rational unit-circle points are the four roots of unity. For a fixed reduced pair $(p,q)$ and a scale $g$, their total Manhattan contribution is

$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$

The parity of $p$ determines whether one axis pair has norm $p$ or $p-2q$. A small example is $(R,r)=(3,1)$, where $(g,p,q,e)=(1,3,1,2)$. The four points are

$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$

$$S(3,1)=3+1+3+3=10.$$

Step 5: Sum the Main Term with Totients and Möbius Inversion

For a fixed $p$, define

$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$

and

$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$

If $M_p=\lfloor N/p\rfloor$, then the full $D=1$ contribution is

$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$

with

$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$

The nontrivial part is the coprime sum $Q(p)$. Möbius inversion removes the gcd condition:

$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$

This is exactly the divisor sum accumulated by the implementation before the final outer summation over $p$.

Step 6: Explicit Correction for $D \gt 1$

All remaining lattice points come from nontrivial rational points on the unit circle, so $D \gt 1$. For each reduced pair $(p,q)$, each primitive Pythagorean triple, and each sign or coordinate-swap symmetry, Step 3 gives a required denominator $\delta$. Whenever $\delta\le M_p$, that one family contributes

$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor.$$

The direct single-pair routine removes coincident points explicitly before summing $|x|+|y|$. For the bulk computation of $T(N)$, the implementations combine the closed form for $D=1$ with this explicit $D \gt 1$ correction. At $N=10^6$, the correction is validated to vanish for $p \gt 12$, so the expensive branch can be truncated safely.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They first build tables for the Möbius function and Euler's totient function up to $N$. Those tables are then used to accumulate the divisor sums needed for $Q(p)$, which yields the dominant arithmetic contribution from the $D=1$ case.

After that, the implementation enumerates primitive Pythagorean triples only for the nontrivial $D \gt 1$ branch. For each admissible reduced pair and each rational unit-circle point generated from the triple symmetries, it computes the reduced denominator $\delta$ and adds the corresponding triangular-scale contribution. A separate direct evaluator for one pair $(R,r)$ is used for checkpoint verification and deduplicates geometric coincidences explicitly.

Complexity Analysis

The sieve for $\mu$ and $\varphi$ is $O(N)$ time and $O(N)$ memory. The divisor accumulation that evaluates the Möbius formula for every $p\le N$ is near-linear, with the usual $O(N\log\log N)$ sieve-like behavior. The $D \gt 1$ correction is much smaller in practice because only small reduced $p$ values survive the validation cutoff and most primitive triples are eliminated immediately by the denominator test. Overall the method is dominated by arithmetic precomputation with $O(N)$ storage.

Footnotes and References

Problem page: https://projecteuler.net/problem=450
Hypocycloid: Wikipedia — Hypocycloid
Rational points on the unit circle and primitive triples: Wikipedia — Pythagorean triple
Möbius inversion: Wikipedia — Möbius inversion formula
Euler's totient function: Wikipedia — Euler's totient function

Problem 450 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <set>
#include <string>
#include <thread>
#include <tuple>
#include <utility>
#include <vector>

using namespace std;

using int64 = long long;
using i128 = __int128_t;

static inline i128 tri_i128(int64 n) { return (i128)n * (n + 1) / 2; }

static string to_string_i128(i128 v) {
    if (v == 0) return "0";
    bool neg = v < 0;
    if (neg) v = -v;
    string s;
    while (v > 0) {
        int digit = (int)(v % 10);
        s.push_back(char('0' + digit));
        v /= 10;
    }
    if (neg) s.push_back('-');
    reverse(s.begin(), s.end());
    return s;
}

static inline i128 iabs128(i128 x) { return x < 0 ? -x : x; }

static i128 gcd128(i128 a, i128 b) {
    if (a < 0) a = -a;
    if (b < 0) b = -b;
    while (b != 0) {
        i128 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

static i128 gcd128_3(i128 a, i128 b, i128 c) {
    return gcd128(gcd128(a, b), c);
}

static inline pair<i128, i128> pow_gauss_i128(i128 x, i128 y, int n) {
    i128 rx = 1, ry = 0;
    for (int i = 0; i < n; ++i) {
        i128 nx = rx * x - ry * y;
        i128 ny = rx * y + ry * x;
        rx = nx;
        ry = ny;
    }
    return {rx, ry};
}

static inline i128 ipow_i128(i128 base, int exp) {
    i128 r = 1;
    for (int i = 0; i < exp; ++i) r *= base;
    return r;
}

static inline tuple<int64, i128, i128, i128> denom_required(int p, int q, int A, int B, int D) {
    int e = p - q;
    auto aq = pow_gauss_i128((i128)A, (i128)B, q);
    auto aec = pow_gauss_i128((i128)A, (i128)(-B), e);

    i128 Deq = ipow_i128((i128)D, e - q);
    i128 ux = (i128)e * aq.first * Deq + (i128)q * aec.first;
    i128 uy = (i128)e * aq.second * Deq + (i128)q * aec.second;

    i128 De = ipow_i128((i128)D, e);
    i128 g0 = gcd128_3(De, ux, uy);
    i128 den128 = De / g0;

    int64 den = (int64)den128;
    return {den, ux, uy, g0};
}

static inline array<pair<int, int>, 8> variations(int a, int b) {
    return array<pair<int, int>, 8>{
        pair<int, int>{a, b},
        pair<int, int>{a, -b},
        pair<int, int>{-a, b},
        pair<int, int>{-a, -b},
        pair<int, int>{b, a},
        pair<int, int>{b, -a},
        pair<int, int>{-b, a},
        pair<int, int>{-b, -a},
    };
}

static vector<tuple<int, int, int>> primitive_triples(int Cmax) {
    vector<tuple<int, int, int>> out;
    int mlim = (int)std::sqrt((long double)Cmax) + 2;
    for (int m = 2; m <= mlim; ++m) {
        for (int n = 1; n < m; ++n) {
            if (((m - n) & 1) == 0) continue;
            if (std::gcd(m, n) != 1) continue;
            int a = m * m - n * n;
            int b = 2 * m * n;
            int c = m * m + n * n;
            if (c > Cmax) continue;
            out.emplace_back(a, b, c);
        }
    }
    return out;
}

static void sieve_mu_phi(int N, vector<int> &mu, vector<int> &phi) {
    mu.assign(N + 1, 0);
    phi.assign(N + 1, 0);
    vector<int> primes;
    vector<char> comp(N + 1, 0);
    mu[1] = 1;
    phi[1] = 1;
    for (int i = 2; i <= N; ++i) {
        if (!comp[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for (int p : primes) {
            long long v = 1LL * i * p;
            if (v > N) break;
            comp[(int)v] = 1;
            if (i % p == 0) {
                phi[(int)v] = phi[i] * p;
                mu[(int)v] = 0;
                break;
            } else {
                phi[(int)v] = phi[i] * (p - 1);
                mu[(int)v] = -mu[i];
            }
        }
    }
}

static vector<int64> compute_B_parallel(int N, const vector<int> &mu, int maxThreads = 0) {
    vector<int64> B(N + 1, 0);

    unsigned hw = thread::hardware_concurrency();
    int T = (maxThreads > 0) ? maxThreads : (hw ? (int)hw : 4);
    T = max(1, min(T, 8));

    vector<vector<int64>> local(T, vector<int64>(N + 1, 0));

    auto worker = [&](int tid, int dStart, int dEnd) {
        for (int d = dStart; d <= dEnd; ++d) {
            int md = mu[d];
            if (md == 0) continue;
            int64 coef = (int64)md * d;
            for (int p = d; p <= N; p += d) {
                int64 m = (p - 1) / (2LL * d);
                local[tid][p] += coef * (m * (m + 1) / 2);
            }
        }
    };

    vector<thread> th;
    th.reserve(T);
    int block = N / T;
    int cur = 1;
    for (int t = 0; t < T; ++t) {
        int start = cur;
        int end = (t == T - 1) ? N : (cur + block - 1);
        cur = end + 1;
        th.emplace_back(worker, t, start, end);
    }
    for (auto &x : th) x.join();

    for (int p = 1; p <= N; ++p) {
        int64 s = 0;
        for (int t = 0; t < T; ++t) s += local[t][p];
        B[p] = s;
    }
    return B;
}

static i128 compute_T(int N) {
    vector<int> mu, phi;
    sieve_mu_phi(N, mu, phi);
    vector<int64> B = compute_B_parallel(N, mu);

    i128 total = 0;

    for (int p = 3; p <= N; ++p) {
        int64 A = phi[p] / 2;
        if (A == 0) continue;
        int64 Bsum = B[p];

        int64 inner;
        if (p & 1) {
            inner = 4LL * p * A - 2LL * Bsum;
        } else {
            if (p % 4 == 0) inner = 4LL * p * A;
            else inner = 4LL * p * A - 4LL * Bsum;
        }

        int64 M = N / p;
        total += (i128)inner * tri_i128(M);
    }

    int Dmax = (int)std::sqrt((long double)(N / 3)) + 2;
    auto triples = primitive_triples(Dmax);

    // For N = 1e6, checking D = 5 (the smallest hypotenuse) shows no D>1
    // contributions for p > 12. We cap extra work here accordingly.
    int maxPextra = 12;

    for (int p = 3; p <= min(N, maxPextra); ++p) {
        int64 M = N / p;
        for (int q = 1; 2 * q < p; ++q) {
            if (std::gcd(p, q) != 1) continue;

            for (auto [a, b, D] : triples) {
                auto vars = variations(a, b);
                for (auto [A, Bv] : vars) {
                    auto [den, ux, uy, g0] = denom_required(p, q, A, Bv, D);
                    if (den > M) continue;

                    i128 bx = ux / g0;
                    i128 by = uy / g0;
                    i128 base_abs = iabs128(bx) + iabs128(by);

                    int64 mmax = M / den;
                    total += base_abs * tri_i128(mmax);
                }
            }
        }
    }

    return total;
}

static pair<i128, size_t> compute_S_and_count(int R, int r) {
    int g = std::gcd(R, r);
    int p = R / g;
    int q = r / g;
    int e = p - q;

    int Dmax = 600;
    auto triples = primitive_triples(Dmax);

    set<pair<long long, long long>> pts;

    auto ipow = [](int k) -> pair<int, int> {
        k %= 4;
        if (k < 0) k += 4;
        switch (k) {
            case 0: return {1, 0};
            case 1: return {0, 1};
            case 2: return {-1, 0};
            default: return {0, -1};
        }
    };

    for (int k = 0; k < 4; ++k) {
        auto wq = ipow(k * q);
        auto wne = ipow(-k * e);
        long long zx = (long long)e * wq.first + (long long)q * wne.first;
        long long zy = (long long)e * wq.second + (long long)q * wne.second;
        pts.insert({(long long)g * zx, (long long)g * zy});
    }

    for (auto [a, b, D] : triples) {
        auto vars = variations(a, b);
        for (auto [A, Bv] : vars) {
            auto [den, ux, uy, g0] = denom_required(p, q, A, Bv, D);
            if (den == 0) continue;
            if (g % den != 0) continue;
            long long k = g / den;

            i128 bx = ux / g0;
            i128 by = uy / g0;
            i128 x = (i128)k * bx;
            i128 y = (i128)k * by;

            pts.insert({(long long)x, (long long)y});
        }
    }

    i128 S = 0;
    for (auto [x, y] : pts) {
        S += (i128)llabs(x) + (i128)llabs(y);
    }
    return {S, pts.size()};
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    {
        auto [S31, c31] = compute_S_and_count(3, 1);
        assert((long long)S31 == 10);
        assert(c31 == 4);

        auto [Sbig, cbig] = compute_S_and_count(2500, 1000);
        assert((long long)Sbig == 24944);
        assert(cbig == 12);
    }
    {
        i128 t3 = compute_T(3);
        i128 t10 = compute_T(10);
        i128 t100 = compute_T(100);
        i128 t1000 = compute_T(1000);
        assert((long long)t3 == 10);
        assert((long long)t10 == 524);
        assert((long long)t100 == 580442);
        assert((long long)t1000 == 583108600);
    }

    i128 ans = compute_T(1'000'000);
    cout << to_string_i128(ans) << "\n";
    return 0;
}

Python

import sys
import math
from multiprocessing import Pool, cpu_count

def gcd(a, b):
    a = abs(a)
    b = abs(b)
    while b:
        a, b = b, a % b
    return a

def gcd3(a, b, c):
    return gcd(gcd(a, b), c)

def pow_gauss(x, y, n):
    rx, ry = 1, 0
    for _ in range(n):
        nx = rx * x - ry * y
        ny = rx * y + ry * x
        rx, ry = nx, ny
    return rx, ry

def denom_required(p, q, A, B, D):
    e = p - q
    aq_x, aq_y = pow_gauss(A, B, q)
    aec_x, aec_y = pow_gauss(A, -B, e)
    
    Deq = D ** (e - q)
    ux = e * aq_x * Deq + q * aec_x
    uy = e * aq_y * Deq + q * aec_y
    
    De = D ** e
    g0 = gcd3(De, ux, uy)
    
    den = De // g0
    return den, ux, uy, g0

def variations(a, b):
    return [
        (a, b), (a, -b), (-a, b), (-a, -b),
        (b, a), (b, -a), (-b, a), (-b, -a)
    ]

def primitive_triples(Cmax):
    out = []
    mlim = int(math.sqrt(Cmax)) + 2
    for m in range(2, mlim + 1):
        for n in range(1, m):
            if (m - n) % 2 == 0: continue
            if math.gcd(m, n) != 1: continue
            a = m * m - n * n
            b = 2 * m * n
            c = m * m + n * n
            if c > Cmax: continue
            out.append((a, b, c))
    return out

def sieve_mu_phi(N):
    mu = [0] * (N + 1)
    phi = [0] * (N + 1)
    primes = []
    comp = bytearray(N + 1)
    
    mu[1] = 1
    phi[1] = 1
    
    for i in range(2, N + 1):
        if not comp[i]:
            primes.append(i)
            phi[i] = i - 1
            mu[i] = -1
        for p in primes:
            v = i * p
            if v > N: break
            comp[v] = 1
            if i % p == 0:
                phi[v] = phi[i] * p
                mu[v] = 0
                break
            else:
                phi[v] = phi[i] * (p - 1)
                mu[v] = -mu[i]
                
    return mu, phi

pool_mu = None
def init_worker(shared_mu):
    global pool_mu
    pool_mu = shared_mu

def worker_B(args):
    d_start, d_end, N = args
    mu = pool_mu
    local = [0] * (N + 1)
    for d in range(d_start, d_end + 1):
        md = mu[d]
        if md == 0:
            continue
        coef = md * d
        for p in range(d, N + 1, d):
            m = (p - 1) // (2 * d)
            local[p] += coef * (m * (m + 1) // 2)
    return local

def compute_B_parallel(N, mu):
    threads = max(1, cpu_count())
    threads = min(threads, 8)
    
    tasks = []
    block = N // threads
    cur = 1
    for t in range(threads):
        start = cur
        end = N if t == threads - 1 else cur + block - 1
        tasks.append((start, end, N))
        cur = end + 1
        
    with Pool(threads, initializer=init_worker, initargs=(mu,)) as pool:
        locals_res = pool.map(worker_B, tasks)
        
    B = [0] * (N + 1)
    for p in range(1, N + 1):
        B[p] = sum(loc[p] for loc in locals_res)
        
    return B

def compute_T(N):
    mu, phi = sieve_mu_phi(N)
    B = compute_B_parallel(N, mu)
    
    total = 0
    for p in range(3, N + 1):
        A = phi[p] // 2
        if A == 0: continue
        Bsum = B[p]
        
        if p % 2 != 0:
            inner = 4 * p * A - 2 * Bsum
        else:
            if p % 4 == 0:
                inner = 4 * p * A
            else:
                inner = 4 * p * A - 4 * Bsum
                
        M = N // p
        total += inner * (M * (M + 1) // 2)
        
    Dmax = int(math.sqrt(N // 3)) + 2
    triples = primitive_triples(Dmax)
    
    maxPextra = 12
    for p in range(3, min(N, maxPextra) + 1):
        M = N // p
        for q in range(1, (p + 1) // 2):
            if math.gcd(p, q) != 1: continue
            
            for a, b, D in triples:
                for A, Bv in variations(a, b):
                    den, ux, uy, g0 = denom_required(p, q, A, Bv, D)
                    if den > M: continue
                    
                    bx = ux // g0
                    by = uy // g0
                    base_abs = abs(bx) + abs(by)
                    
                    mmax = M // den
                    total += base_abs * (mmax * (mmax + 1) // 2)
                    
    return total

def solve():
    return str(compute_T(1000000))

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler450 {

    static BigInteger gcdBig(BigInteger a, BigInteger b) {
        a = a.abs();
        b = b.abs();
        while (!b.equals(BigInteger.ZERO)) {
            BigInteger t = a.remainder(b);
            a = b;
            b = t;
        }
        return a;
    }

    static BigInteger[] powGauss(BigInteger x, BigInteger y, int n) {
        BigInteger rx = BigInteger.ONE, ry = BigInteger.ZERO;
        for (int i = 0; i < n; i++) {
            BigInteger nx = rx.multiply(x).subtract(ry.multiply(y));
            BigInteger ny = rx.multiply(y).add(ry.multiply(x));
            rx = nx;
            ry = ny;
        }
        return new BigInteger[] { rx, ry };
    }

    static BigInteger triBig(long n) {
        return BigInteger.valueOf(n).multiply(BigInteger.valueOf(n + 1)).divide(BigInteger.valueOf(2));
    }

    static class DenomResult {
        long den;
        BigInteger ux, uy, g0;

        DenomResult(long den, BigInteger ux, BigInteger uy, BigInteger g0) {
            this.den = den;
            this.ux = ux;
            this.uy = uy;
            this.g0 = g0;
        }
    }

    static DenomResult denomRequired(int p, int q, int A, int B, int D) {
        int e = p - q;
        BigInteger[] aq = powGauss(BigInteger.valueOf(A), BigInteger.valueOf(B), q);
        BigInteger[] aec = powGauss(BigInteger.valueOf(A), BigInteger.valueOf(-B), e);

        BigInteger Deq = BigInteger.valueOf(D).pow(e - q);
        BigInteger ux = BigInteger.valueOf(e).multiply(aq[0]).multiply(Deq).add(BigInteger.valueOf(q).multiply(aec[0]));
        BigInteger uy = BigInteger.valueOf(e).multiply(aq[1]).multiply(Deq).add(BigInteger.valueOf(q).multiply(aec[1]));

        BigInteger De = BigInteger.valueOf(D).pow(e);
        BigInteger g0 = gcdBig(gcdBig(De, ux), uy);

        BigInteger denBig = De.divide(g0);
        return new DenomResult(denBig.longValue(), ux, uy, g0);
    }

    static int[][] variations(int a, int b) {
        return new int[][] {
                { a, b }, { a, -b }, { -a, b }, { -a, -b },
                { b, a }, { b, -a }, { -b, a }, { -b, -a }
        };
    }

    static class Triple {
        int a, b, c;

        Triple(int a, int b, int c) {
            this.a = a;
            this.b = b;
            this.c = c;
        }
    }

    static List<Triple> primitiveTriples(int Cmax) {
        List<Triple> out = new ArrayList<>();
        int mlim = (int) Math.sqrt(Cmax) + 2;
        for (int m = 2; m <= mlim; m++) {
            for (int n = 1; n < m; n++) {
                if (((m - n) & 1) == 0)
                    continue;
                if (gcdInt(m, n) != 1)
                    continue;
                int a = m * m - n * n;
                int b = 2 * m * n;
                int c = m * m + n * n;
                if (c > Cmax)
                    continue;
                out.add(new Triple(a, b, c));
            }
        }
        return out;
    }

    static int gcdInt(int a, int b) {
        while (b != 0) {
            int t = a % b;
            a = b;
            b = t;
        }
        return a;
    }

    static class SieveResult {
        byte[] mu;
        int[] phi;
    }

    static SieveResult sieveMuPhi(int N) {
        SieveResult res = new SieveResult();
        res.mu = new byte[N + 1];
        res.phi = new int[N + 1];
        boolean[] comp = new boolean[N + 1];
        List<Integer> primes = new ArrayList<>();

        res.mu[1] = 1;
        res.phi[1] = 1;

        for (int i = 2; i <= N; i++) {
            if (!comp[i]) {
                primes.add(i);
                res.phi[i] = i - 1;
                res.mu[i] = -1;
            }
            for (int p : primes) {
                long v = (long) i * p;
                if (v > N)
                    break;
                comp[(int) v] = true;
                if (i % p == 0) {
                    res.phi[(int) v] = res.phi[i] * p;
                    res.mu[(int) v] = 0;
                    break;
                } else {
                    res.phi[(int) v] = res.phi[i] * (p - 1);
                    res.mu[(int) v] = (byte) -res.mu[i];
                }
            }
        }
        return res;
    }

    static long[] computeBParallel(int N, byte[] mu) {
        int threads = Math.max(1, Math.min(8, Runtime.getRuntime().availableProcessors()));

        long[][] local = new long[threads][N + 1];
        int block = N / threads;

        IntStream.range(0, threads).parallel().forEach(t -> {
            int start = (t == 0) ? 1 : t * block + 1;
            int end = (t == threads - 1) ? N : (t + 1) * block;

            for (int d = start; d <= end; d++) {
                int md = mu[d];
                if (md == 0)
                    continue;
                long coef = (long) md * d;
                for (int p = d; p <= N; p += d) {
                    long m = (p - 1) / (2L * d);
                    local[t][p] += coef * (m * (m + 1) / 2);
                }
            }
        });

        long[] B = new long[N + 1];
        for (int p = 1; p <= N; p++) {
            long s = 0;
            for (int t = 0; t < threads; t++) {
                s += local[t][p];
            }
            B[p] = s;
        }

        return B;
    }

    public static String solve() {
        int N = 1000000;
        SieveResult sieve = sieveMuPhi(N);
        long[] B = computeBParallel(N, sieve.mu);

        BigInteger total = BigInteger.ZERO;

        for (int p = 3; p <= N; p++) {
            long A = sieve.phi[p] / 2;
            if (A == 0)
                continue;
            long Bsum = B[p];

            long inner;
            if ((p & 1) != 0) {
                inner = 4L * p * A - 2L * Bsum;
            } else {
                if (p % 4 == 0) {
                    inner = 4L * p * A;
                } else {
                    inner = 4L * p * A - 4L * Bsum;
                }
            }

            long M = N / p;
            total = total.add(BigInteger.valueOf(inner).multiply(triBig(M)));
        }

        int Dmax = (int) Math.sqrt(N / 3) + 2;
        List<Triple> triples = primitiveTriples(Dmax);

        int maxPextra = 12;

        for (int p = 3; p <= Math.min(N, maxPextra); p++) {
            long M = N / p;
            for (int q = 1; 2 * q < p; q++) {
                if (gcdInt(p, q) != 1)
                    continue;

                for (Triple t : triples) {
                    int[][] vars = variations(t.a, t.b);
                    for (int[] var : vars) {
                        DenomResult dr = denomRequired(p, q, var[0], var[1], t.c);
                        if (dr.den > M)
                            continue;

                        BigInteger bx = dr.ux.divide(dr.g0);
                        BigInteger by = dr.uy.divide(dr.g0);
                        BigInteger baseAbs = bx.abs().add(by.abs());

                        long mmax = M / dr.den;
                        total = total.add(baseAbs.multiply(triBig(mmax)));
                    }
                }
            }
        }

        return total.toString();
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}