Problem 572: Idempotent Matrices

Project Euler Problem 572 Solution

EulerSolve provides an optimized solution for Project Euler Problem 572, Idempotent Matrices, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count integer $3\times3$ matrices $M$ whose entries lie in $[-n,n]$ and satisfy $$M^2=M.$$ Such matrices are idempotent. A direct search over all $ (2n+1)^9 $ matrices is hopeless, so the solution instead uses the structure theorem for idempotents: every such matrix is a projection, and in dimension $3$ that reduces the problem to counting lattice points on a plane. Mathematical Approach The implementations classify idempotent matrices by rank, express the nontrivial cases through one primitive integer direction vector and one integer linear functional, and then count the admissible integer solutions of a single equation. Step 1: Split the problem by rank If $M^2=M$, then every eigenvalue $\lambda$ satisfies $\lambda^2=\lambda$, so $\lambda\in\{0,1\}$. Therefore an idempotent $3\times3$ matrix can only have rank $0,1,2,$ or $3$. The extreme cases are immediate: $$\operatorname{rank}(M)=0 \Rightarrow M=0,\qquad \operatorname{rank}(M)=3 \Rightarrow M=I.$$ So for positive $n$, two matrices are already known, and the real work is to count the rank-$1$ and rank-$2$ cases. Step 2: Parametrize every rank-$1$ idempotent Let $M$ have rank $1$. Its image is a one-dimensional sublattice of $\mathbb{Z}^3$, so it is generated by a primitive integer vector $\mathbf{p}\in\mathbb{Z}^3$....

Detailed mathematical approach

Problem Summary

We must count integer $3\times3$ matrices $M$ whose entries lie in $[-n,n]$ and satisfy

$$M^2=M.$$

Such matrices are idempotent. A direct search over all $ (2n+1)^9 $ matrices is hopeless, so the solution instead uses the structure theorem for idempotents: every such matrix is a projection, and in dimension $3$ that reduces the problem to counting lattice points on a plane.

Mathematical Approach

The implementations classify idempotent matrices by rank, express the nontrivial cases through one primitive integer direction vector and one integer linear functional, and then count the admissible integer solutions of a single equation.

Step 1: Split the problem by rank

If $M^2=M$, then every eigenvalue $\lambda$ satisfies $\lambda^2=\lambda$, so $\lambda\in\{0,1\}$. Therefore an idempotent $3\times3$ matrix can only have rank $0,1,2,$ or $3$.

The extreme cases are immediate:

$$\operatorname{rank}(M)=0 \Rightarrow M=0,\qquad \operatorname{rank}(M)=3 \Rightarrow M=I.$$

So for positive $n$, two matrices are already known, and the real work is to count the rank-$1$ and rank-$2$ cases.

Step 2: Parametrize every rank-$1$ idempotent

Let $M$ have rank $1$. Its image is a one-dimensional sublattice of $\mathbb{Z}^3$, so it is generated by a primitive integer vector $\mathbf{p}\in\mathbb{Z}^3$. Because $M$ is idempotent, it acts as the identity on its image, hence there exists an integer row vector $\mathbf{w}^{T}$ such that

$$M=\mathbf{p}\mathbf{w}^{T},\qquad \mathbf{w}\cdot\mathbf{p}=1.$$

The primitive condition on $\mathbf{p}$ is forced: if $\gcd(p_1,p_2,p_3)=d>1$, then $\mathbf{w}\cdot\mathbf{p}$ would also be divisible by $d$, contradicting $\mathbf{w}\cdot\mathbf{p}=1$.

Conversely, any integer pair $(\mathbf{p},\mathbf{w})$ with $\mathbf{w}\cdot\mathbf{p}=1$ gives an idempotent, because

$$\left(\mathbf{p}\mathbf{w}^{T}\right)^2=\mathbf{p}\left(\mathbf{w}\cdot\mathbf{p}\right)\mathbf{w}^{T}=\mathbf{p}\mathbf{w}^{T}.$$

The pair $(\mathbf{p},\mathbf{w})$ and $(-\mathbf{p},-\mathbf{w})$ produces the same matrix, so we fix a canonical orientation by requiring the first nonzero component of $\mathbf{p}$ to be positive.

Step 3: Parametrize every rank-$2$ idempotent

If $M$ has rank $2$, then $I-M$ is idempotent of rank $1$. Applying the previous description to $I-M$ yields

$$M=I-\mathbf{p}\mathbf{w}^{T},\qquad \mathbf{w}\cdot\mathbf{p}=1.$$

This is again a complete characterization, since

$$\left(I-\mathbf{p}\mathbf{w}^{T}\right)^2=I-2\mathbf{p}\mathbf{w}^{T}+\mathbf{p}\left(\mathbf{w}\cdot\mathbf{p}\right)\mathbf{w}^{T}=I-\mathbf{p}\mathbf{w}^{T}.$$

Thus both nontrivial ranks are described by the same arithmetic data; only the entry bounds are different.

Step 4: Turn the entry bounds into intervals for $\mathbf{w}$

Write $\mathbf{p}=(p_1,p_2,p_3)^{T}$ and $\mathbf{w}=(w_1,w_2,w_3)^{T}$.

For rank $1$, every entry is $M_{ij}=p_iw_j$, so

$$|p_iw_j|\le n \qquad (1\le i,j\le 3).$$

If $P=\max(|p_1|,|p_2|,|p_3|)$, then each coordinate of $\mathbf{w}$ must satisfy

$$|w_j|\le \left\lfloor\frac{n}{P}\right\rfloor.$$

So rank $1$ is reduced to counting integer solutions of

$$p_1w_1+p_2w_2+p_3w_3=1$$

inside a symmetric cube.

For rank $2$, the off-diagonal entries are $-p_iw_j$ for $i\ne j$. Therefore, for each column $j$,

$$|w_j|\le \left\lfloor\frac{n}{\max_{i\ne j}|p_i|}\right\rfloor$$

whenever at least one of the other two coefficients is nonzero. The diagonal entries are $1-p_jw_j$, so

$$|1-p_jw_j|\le n \iff 1-n\le p_jw_j\le 1+n.$$

When $p_j\ne 0$, this becomes an explicit integer interval via floor and ceiling division. The admissible values of $w_j$ are obtained by intersecting the diagonal and off-diagonal restrictions for each coordinate.

Step 5: Count lattice points on the plane $\mathbf{w}\cdot\mathbf{p}=1$

Once the three coordinate ranges are known, the problem becomes: count integer points in a rectangular box that lie on

$$p_1w_1+p_2w_2+p_3w_3=1.$$

When one or two coefficients vanish, the equation collapses to an easy one-variable or two-variable count. In the generic case, one coordinate is solved from the equation, the other two are enumerated inside their intervals, and the remaining coordinate is checked for integrality and for being within range.

This turns a nine-variable matrix search into a much smaller geometric counting problem on one affine plane.

Worked Example: $n=1$ and $\mathbf{p}=(1,0,0)^{T}$

For rank $1$, the equation $\mathbf{w}\cdot\mathbf{p}=1$ becomes $w_1=1$. The box constraints give $|w_2|\le 1$ and $|w_3|\le 1$, so this one primitive direction contributes

$$3\cdot 3=9$$

rank-$1$ matrices.

For rank $2$, we have $M=I-\mathbf{p}\mathbf{w}^{T}$. The same equation still forces $w_1=1$, the off-diagonal bounds again give $|w_2|,|w_3|\le 1$, and the diagonal entry becomes $1-w_1=0$, which is allowed. So the same direction contributes another $9$ rank-$2$ matrices.

Summing over all canonical primitive directions and then adding $0$ and $I$ gives $C(1)=164$, which matches the small test value used by the implementation.

How the Code Works

The C++, Python, and Java implementations all use this rank decomposition. They begin with the trivial idempotents, enumerate all nonzero primitive direction vectors in $[-n,n]^3$ with a canonical sign choice, and for each direction count two families: rank $1$ matrices in a symmetric cube and rank $2$ matrices in the intersection of three coordinate intervals.

The shared core is a bounded linear Diophantine solver for

$$p_1w_1+p_2w_2+p_3w_3=1.$$

Special cases with zero coefficients are treated separately, while the generic case fixes one coordinate from the equation and loops over the other two. The outer enumeration is split across multiple workers, so independent slices of the primitive direction set can be counted in parallel and added safely at the end.

Complexity Analysis

The outer search box contains $O(n^3)$ candidate direction vectors. For each surviving primitive canonical vector, the counting stage is a bounded two-dimensional scan of integer pairs, so a conservative worst-case bound is $O(n^5)$ time. That upper bound is crude, because large coefficients make the admissible intervals much smaller and many rank-$2$ boxes are empty or very thin. The auxiliary memory is $O(1)$ per worker, plus a small amount of storage for partial sums.

Footnotes and References

Problem page: https://projecteuler.net/problem=572
Idempotent matrix: Wikipedia — Idempotent matrix
Projection (linear algebra): Wikipedia — Projection (linear algebra)
Linear Diophantine equation: Wikipedia — Linear Diophantine equation
Rank-nullity theorem: Wikipedia — Rank-nullity theorem

Problem 572 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

// Project Euler 572: Idempotent Matrices
//
// An integer 3x3 matrix M is idempotent iff it is a projection of the free
// abelian group Z^3, hence its rank over Q is 0,1,2,3. Rank 0 and 3 give
// the zero and identity matrices.
//
// Rank 1 idempotents have the form M = u v^T with u,v in Z^3 and v·u = 1.
// Rank 2 idempotents have the form M = I - u v^T with the same constraint.
// The pair (u,v) and (-u,-v) give the same matrix, so we count only vectors u
// with a canonical sign (first nonzero component positive) and count matching v.
//
// For rank 1, the bound |M_ij| <= n is equivalent to |u_i v_j| <= n for all i,j.
// For rank 2, bounds are |u_i v_j| <= n for i != j and |1 - u_i v_i| <= n.

using i64 = std::int64_t;
using u64 = std::uint64_t;

static inline int iabs(int x) { return x < 0 ? -x : x; }

static inline int gcd3(int a, int b, int c) {
    a = iabs(a);
    b = iabs(b);
    c = iabs(c);
    return std::gcd(a, std::gcd(b, c));
}

static inline bool canonical_u(int a, int b, int c) {
    if (a != 0) return a > 0;
    if (b != 0) return b > 0;
    return c > 0;  // (0,0,0) excluded by caller
}

static inline i64 floor_div(i64 a, i64 b) {
    // floor(a/b) for b != 0
    i64 q = a / b;
    i64 r = a % b;
    if (r != 0 && ((r > 0) != (b > 0))) --q;
    return q;
}

static inline i64 ceil_div(i64 a, i64 b) {
    // ceil(a/b) for b != 0
    return -floor_div(-a, b);
}

static inline u64 count_solutions_symmetric(int a, int b, int c, int V) {
    // Count integer (x,y,z) with -V<=x,y,z<=V and a*x + b*y + c*z = 1.
    if (V < 0) return 0;
    const int len = 2 * V + 1;

    const bool za = (a == 0), zb = (b == 0), zc = (c == 0);
    const int zcnt = int(za) + int(zb) + int(zc);

    if (zcnt == 2) {
        // One variable only.
        const int coef = !za ? a : (!zb ? b : c);
        if (coef == 0) return 0;
        if (1 % coef != 0) return 0;
        const int x = 1 / coef;
        if (iabs(x) > V) return 0;
        return static_cast<u64>(len) * static_cast<u64>(len);
    }
    if (zcnt == 1) {
        // Two variables.
        if (c == 0) {
            // a*x + b*y = 1, z free.
            u64 cnt = 0;
            for (int x = -V; x <= V; ++x) {
                const int rem = 1 - a * x;
                if (rem % b != 0) continue;
                const int y = rem / b;
                if (iabs(y) <= V) cnt += static_cast<u64>(len);
            }
            return cnt;
        }
        if (b == 0) {
            u64 cnt = 0;
            for (int x = -V; x <= V; ++x) {
                const int rem = 1 - a * x;
                if (rem % c != 0) continue;
                const int z = rem / c;
                if (iabs(z) <= V) cnt += static_cast<u64>(len);
            }
            return cnt;
        }
        // a==0
        u64 cnt = 0;
        for (int y = -V; y <= V; ++y) {
            const int rem = 1 - b * y;
            if (rem % c != 0) continue;
            const int z = rem / c;
            if (iabs(z) <= V) cnt += static_cast<u64>(len);
        }
        return cnt;
    }

    // Three-variable case. Pick pivot to minimize divisor checks: prefer +/-1, else largest |coef|.
    int coef[3] = {a, b, c};
    int pivot = 0;
    for (int k = 0; k < 3; ++k) {
        if (iabs(coef[k]) == 1) {
            pivot = k;
            break;
        }
        if (iabs(coef[k]) > iabs(coef[pivot])) pivot = k;
    }
    const int i = (pivot + 1) % 3;
    const int j = (pivot + 2) % 3;

    u64 cnt = 0;
    for (int xi = -V; xi <= V; ++xi) {
        for (int xj = -V; xj <= V; ++xj) {
            const int rem = 1 - coef[i] * xi - coef[j] * xj;
            const int den = coef[pivot];
            if (rem % den != 0) continue;
            const int xp = rem / den;
            if (iabs(xp) <= V) ++cnt;
        }
    }
    return cnt;
}

struct Range {
    int lo = 0;
    int hi = -1;
    int len() const { return hi >= lo ? (hi - lo + 1) : 0; }
};

static inline u64 count_solutions_box(int a, int b, int c, Range r0, Range r1, Range r2) {
    // Count (x,y,z) in r0 x r1 x r2 with a*x + b*y + c*z = 1.
    if (r0.len() == 0 || r1.len() == 0 || r2.len() == 0) return 0;

    int coef[3] = {a, b, c};
    Range rr[3] = {r0, r1, r2};

    // Handle degenerate cases with zero coefficients.
    const bool za = (a == 0), zb = (b == 0), zc = (c == 0);
    const int zcnt = int(za) + int(zb) + int(zc);
    if (zcnt == 3) return 0;

    if (zcnt == 2) {
        // One variable only.
        int k = !za ? 0 : (!zb ? 1 : 2);
        const int den = coef[k];
        if (1 % den != 0) return 0;
        const int x = 1 / den;
        return (x < rr[k].lo || x > rr[k].hi) ? 0ULL
                                              : static_cast<u64>(rr[(k + 1) % 3].len()) *
                                                    static_cast<u64>(rr[(k + 2) % 3].len());
    }

    if (zcnt == 1) {
        // Two variables + one free.
        int free_k = (za ? 0 : (zb ? 1 : 2));
        int i = (free_k + 1) % 3;
        int j = (free_k + 2) % 3;
        // Solve coef[i]*xi + coef[j]*xj = 1, x_free arbitrary.
        u64 cnt2 = 0;
        for (int xi = rr[i].lo; xi <= rr[i].hi; ++xi) {
            const int rem = 1 - coef[i] * xi;
            const int den = coef[j];
            if (rem % den != 0) continue;
            const int xj = rem / den;
            if (xj >= rr[j].lo && xj <= rr[j].hi) ++cnt2;
        }
        return cnt2 * static_cast<u64>(rr[free_k].len());
    }

    // Choose pivot with nonzero coefficient to minimize iterations on the other two ranges.
    int best_pivot = -1;
    i64 best_cost = (i64)1 << 62;
    for (int k = 0; k < 3; ++k) {
        if (coef[k] == 0) continue;
        int i = (k + 1) % 3;
        int j = (k + 2) % 3;
        i64 cost = static_cast<i64>(rr[i].len()) * static_cast<i64>(rr[j].len());
        // Prefer +/-1 pivots when tied.
        if (cost < best_cost || (cost == best_cost && iabs(coef[k]) == 1)) {
            best_cost = cost;
            best_pivot = k;
        }
    }
    const int k = best_pivot;
    const int i = (k + 1) % 3;
    const int j = (k + 2) % 3;

    u64 cnt = 0;
    for (int xi = rr[i].lo; xi <= rr[i].hi; ++xi) {
        for (int xj = rr[j].lo; xj <= rr[j].hi; ++xj) {
            const int rem = 1 - coef[i] * xi - coef[j] * xj;
            const int den = coef[k];
            if (rem % den != 0) continue;
            const int xk = rem / den;
            if (xk >= rr[k].lo && xk <= rr[k].hi) ++cnt;
        }
    }
    return cnt;
}

static inline Range intersect(Range a, Range b) {
    Range r;
    r.lo = std::max(a.lo, b.lo);
    r.hi = std::min(a.hi, b.hi);
    return r;
}

static inline u64 compute_C(int n) {
    // base: zero matrix always; identity only if n>=1
    const u64 base = 1 + (n >= 1 ? 1 : 0);

    const int threads = std::max(1u, std::min(8u, std::thread::hardware_concurrency() ? std::thread::hardware_concurrency() : 4u));
    std::vector<std::thread> pool;
    std::vector<u64> part_r1(threads, 0), part_r2(threads, 0);

    const int a_min = -n;
    const int a_max = n;
    const int total_a = a_max - a_min + 1;
    const int block = (total_a + threads - 1) / threads;

    for (int t = 0; t < threads; ++t) {
        const int alo = a_min + t * block;
        const int ahi = std::min(a_max, alo + block - 1);
        pool.emplace_back([&, t, alo, ahi] {
            u64 r1 = 0, r2 = 0;
            for (int a = alo; a <= ahi; ++a) {
                for (int b = -n; b <= n; ++b) {
                    for (int c = -n; c <= n; ++c) {
                        if (a == 0 && b == 0 && c == 0) continue;
                        if (!canonical_u(a, b, c)) continue;
                        if (gcd3(a, b, c) != 1) continue;

                        const int Umax = std::max({iabs(a), iabs(b), iabs(c)});
                        const int V = n / Umax;
                        if (V > 0) {
                            r1 += count_solutions_symmetric(a, b, c, V);
                        }

                        // Rank 2: bounds for v0,v1,v2.
                        Range rv[3] = {{-n - 1, n + 1}, {-n - 1, n + 1}, {-n - 1, n + 1}};
                        const int u[3] = {a, b, c};
                        for (int j = 0; j < 3; ++j) {
                            int max_other = 0;
                            for (int i = 0; i < 3; ++i) {
                                if (i == j) continue;
                                max_other = std::max(max_other, iabs(u[i]));
                            }
                            if (max_other > 0) {
                                const int off = n / max_other;
                                rv[j] = intersect(rv[j], Range{-off, off});
                            }
                            if (u[j] != 0) {
                                // Diagonal constraint: 1-n <= u[j]*v[j] <= 1+n (note sign flip if u[j] < 0).
                                i64 lo, hi;
                                if (u[j] > 0) {
                                    lo = ceil_div(1LL - n, u[j]);
                                    hi = floor_div(1LL + n, u[j]);
                                } else {
                                    lo = ceil_div(1LL + n, u[j]);
                                    hi = floor_div(1LL - n, u[j]);
                                }
                                rv[j] = intersect(rv[j], Range{static_cast<int>(lo), static_cast<int>(hi)});
                            }
                        }
                        r2 += count_solutions_box(a, b, c, rv[0], rv[1], rv[2]);
                    }
                }
            }
            part_r1[t] = r1;
            part_r2[t] = r2;
        });
    }
    for (auto& th : pool) th.join();

    const u64 rank1 = std::accumulate(part_r1.begin(), part_r1.end(), 0ULL);
    const u64 rank2 = std::accumulate(part_r2.begin(), part_r2.end(), 0ULL);
    return base + rank1 + rank2;
}

int main() {
    // Validation from the statement.
    if (compute_C(1) != 164ULL) {
        std::cerr << "Validation failed: C(1)\n";
        return 1;
    }
    if (compute_C(2) != 848ULL) {
        std::cerr << "Validation failed: C(2)\n";
        return 1;
    }

    std::cout << compute_C(200) << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler572 {
    static int gcd3(int a, int b, int c) {
        a = Math.abs(a);
        b = Math.abs(b);
        c = Math.abs(c);
        long g = a;
        long h = b;
        while (h != 0) {
            long t = g % h;
            g = h;
            h = t;
        }
        long r = c;
        while (r != 0) {
            long t = g % r;
            g = r;
            r = t;
        }
        return (int) g;
    }

    static boolean canonicalU(int a, int b, int c) {
        if (a != 0)
            return a > 0;
        if (b != 0)
            return b > 0;
        return c > 0;
    }

    static long floorDiv(long a, long b) {
        long q = a / b;
        long r = a % b;
        if (r != 0 && ((r > 0) != (b > 0)))
            q--;
        return q;
    }

    static long ceilDiv(long a, long b) {
        return -floorDiv(-a, b);
    }

    static long countSolutionsSymmetric(int a, int b, int c, int V) {
        if (V < 0)
            return 0;
        int len = 2 * V + 1;

        boolean za = (a == 0), zb = (b == 0), zc = (c == 0);
        int zcnt = (za ? 1 : 0) + (zb ? 1 : 0) + (zc ? 1 : 0);

        if (zcnt == 2) {
            int coef = !za ? a : (!zb ? b : c);
            if (coef == 0)
                return 0;
            if (1 % coef != 0)
                return 0;
            int x = 1 / coef;
            if (Math.abs(x) > V)
                return 0;
            return (long) len * len;
        }

        if (zcnt == 1) {
            if (c == 0) {
                long cnt = 0;
                for (int x = -V; x <= V; x++) {
                    int rem = 1 - a * x;
                    if (rem % b != 0)
                        continue;
                    int y = rem / b;
                    if (Math.abs(y) <= V)
                        cnt += len;
                }
                return cnt;
            }
            if (b == 0) {
                long cnt = 0;
                for (int x = -V; x <= V; x++) {
                    int rem = 1 - a * x;
                    if (rem % c != 0)
                        continue;
                    int z = rem / c;
                    if (Math.abs(z) <= V)
                        cnt += len;
                }
                return cnt;
            }
            long cnt = 0;
            for (int y = -V; y <= V; y++) {
                int rem = 1 - b * y;
                if (rem % c != 0)
                    continue;
                int z = rem / c;
                if (Math.abs(z) <= V)
                    cnt += len;
            }
            return cnt;
        }

        int[] coef = { a, b, c };
        int pivot = 0;
        for (int k = 0; k < 3; k++) {
            if (Math.abs(coef[k]) == 1) {
                pivot = k;
                break;
            }
            if (Math.abs(coef[k]) > Math.abs(coef[pivot]))
                pivot = k;
        }
        int i = (pivot + 1) % 3;
        int j = (pivot + 2) % 3;

        long cnt = 0;
        for (int xi = -V; xi <= V; xi++) {
            for (int xj = -V; xj <= V; xj++) {
                int rem = 1 - coef[i] * xi - coef[j] * xj;
                int den = coef[pivot];
                if (rem % den != 0)
                    continue;
                int xp = rem / den;
                if (Math.abs(xp) <= V)
                    cnt++;
            }
        }
        return cnt;
    }

    static class Range {
        int lo, hi;

        Range(int lo, int hi) {
            this.lo = lo;
            this.hi = hi;
        }

        int len() {
            return hi >= lo ? (hi - lo + 1) : 0;
        }
    }

    static Range intersect(Range a, Range b) {
        return new Range(Math.max(a.lo, b.lo), Math.min(a.hi, b.hi));
    }

    static long countSolutionsBox(int a, int b, int c, Range r0, Range r1, Range r2) {
        if (r0.len() == 0 || r1.len() == 0 || r2.len() == 0)
            return 0;

        int[] coef = { a, b, c };
        Range[] rr = { r0, r1, r2 };

        boolean za = (a == 0), zb = (b == 0), zc = (c == 0);
        int zcnt = (za ? 1 : 0) + (zb ? 1 : 0) + (zc ? 1 : 0);
        if (zcnt == 3)
            return 0;

        if (zcnt == 2) {
            int k = !za ? 0 : (!zb ? 1 : 2);
            int den = coef[k];
            if (1 % den != 0)
                return 0;
            int x = 1 / den;
            return (x < rr[k].lo || x > rr[k].hi) ? 0 : (long) rr[(k + 1) % 3].len() * rr[(k + 2) % 3].len();
        }

        if (zcnt == 1) {
            int freeK = za ? 0 : (zb ? 1 : 2);
            int i = (freeK + 1) % 3;
            int j = (freeK + 2) % 3;
            long cnt2 = 0;
            for (int xi = rr[i].lo; xi <= rr[i].hi; xi++) {
                int rem = 1 - coef[i] * xi;
                int den = coef[j];
                if (rem % den != 0)
                    continue;
                int xj = rem / den;
                if (xj >= rr[j].lo && xj <= rr[j].hi)
                    cnt2++;
            }
            return cnt2 * rr[freeK].len();
        }

        int bestPivot = -1;
        long bestCost = Long.MAX_VALUE;
        for (int k = 0; k < 3; k++) {
            if (coef[k] == 0)
                continue;
            int i = (k + 1) % 3;
            int j = (k + 2) % 3;
            long cost = (long) rr[i].len() * rr[j].len();
            if (cost < bestCost || (cost == bestCost && Math.abs(coef[k]) == 1)) {
                bestCost = cost;
                bestPivot = k;
            }
        }

        int k = bestPivot;
        int i = (k + 1) % 3;
        int j = (k + 2) % 3;

        long cnt = 0;
        for (int xi = rr[i].lo; xi <= rr[i].hi; xi++) {
            for (int xj = rr[j].lo; xj <= rr[j].hi; xj++) {
                int rem = 1 - coef[i] * xi - coef[j] * xj;
                int den = coef[k];
                if (rem % den != 0)
                    continue;
                int xk = rem / den;
                if (xk >= rr[k].lo && xk <= rr[k].hi)
                    cnt++;
            }
        }
        return cnt;
    }

    static long computeC(int n) {
        long base = 1 + (n >= 1 ? 1 : 0);

        int threads = Runtime.getRuntime().availableProcessors();
        if (threads <= 0)
            threads = 1;
        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Callable<long[]>> tasks = new ArrayList<>();

        int aMin = -n;
        int aMax = n;
        int totalA = aMax - aMin + 1;
        int block = (totalA + threads - 1) / threads;

        for (int t = 0; t < threads; t++) {
            final int alo = aMin + t * block;
            final int ahi = Math.min(aMax, alo + block - 1);
            if (alo > ahi)
                continue;
            tasks.add(() -> {
                long r1 = 0, r2 = 0;
                for (int a = alo; a <= ahi; a++) {
                    for (int b = -n; b <= n; b++) {
                        for (int c = -n; c <= n; c++) {
                            if (a == 0 && b == 0 && c == 0)
                                continue;
                            if (!canonicalU(a, b, c))
                                continue;
                            if (gcd3(a, b, c) != 1)
                                continue;

                            int Umax = Math.max(Math.abs(a), Math.max(Math.abs(b), Math.abs(c)));
                            int V = n / Umax;
                            if (V > 0) {
                                r1 += countSolutionsSymmetric(a, b, c, V);
                            }

                            Range[] rv = {
                                    new Range(-n - 1, n + 1),
                                    new Range(-n - 1, n + 1),
                                    new Range(-n - 1, n + 1)
                            };

                            int[] u = { a, b, c };
                            for (int j = 0; j < 3; j++) {
                                int maxOther = 0;
                                for (int i = 0; i < 3; i++) {
                                    if (i == j)
                                        continue;
                                    maxOther = Math.max(maxOther, Math.abs(u[i]));
                                }
                                if (maxOther > 0) {
                                    int off = n / maxOther;
                                    rv[j] = intersect(rv[j], new Range(-off, off));
                                }
                                if (u[j] != 0) {
                                    long lo, hi;
                                    if (u[j] > 0) {
                                        lo = ceilDiv(1L - n, u[j]);
                                        hi = floorDiv(1L + n, u[j]);
                                    } else {
                                        lo = ceilDiv(1L + n, u[j]);
                                        hi = floorDiv(1L - n, u[j]);
                                    }
                                    rv[j] = intersect(rv[j], new Range((int) lo, (int) hi));
                                }
                            }
                            r2 += countSolutionsBox(a, b, c, rv[0], rv[1], rv[2]);
                        }
                    }
                }
                return new long[] { r1, r2 };
            });
        }

        long rank1 = 0;
        long rank2 = 0;
        try {
            for (Future<long[]> res : executor.invokeAll(tasks)) {
                long[] part = res.get();
                rank1 += part[0];
                rank2 += part[1];
            }
        } catch (Exception e) {
        }
        executor.shutdown();

        return base + rank1 + rank2;
    }

    public static String solve() {
        return String.valueOf(computeC(200));
    }

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