#include <cassert> #include <cstdio> #include <cstring> #include <cstdlib> #include <cmath> #include <ctime> #include <cctype> #include <algorithm> #include <random> #include <bitset> #include <queue> #include <functional> #include <set> #include <map> #include <vector> #include <chrono> #include <iostream> #include <iomanip> #include <limits> #include <numeric> typedef long long ll; typedef long double ld; // quaternion by Haoxiang Yu: https://yhx-12243.github.io/OI-transit/records/cf1375I.html struct quaternion { ld s, x, y, z; quaternion (ld s0 = 0, ld x0 = 0, ld y0 = 0, ld z0 = 0) : s(s0), x(x0), y(y0), z(z0) {} inline quaternion operator - () const {return quaternion(-s, -x, -y, -z);} inline quaternion conj() const {return quaternion(s, -x, -y, -z);} inline quaternion operator + (const quaternion &B) const {return quaternion(s + B.s, x + B.x, y + B.y, z + B.z);} inline quaternion operator - (const quaternion &B) const {return quaternion(s - B.s, x - B.x, y - B.y, z - B.z);} inline quaternion operator * (const ld k) const {return quaternion(s * k, x * k, y * k, z * k);} inline quaternion operator * (const quaternion &B) const { return quaternion( s * B.s - x * B.x - y * B.y - z * B.z, s * B.x + x * B.s + y * B.z - z * B.y, s * B.y - x * B.z + y * B.s + z * B.x, s * B.z + x * B.y - y * B.x + z * B.s ); } inline operator void * () const {return s || x || y || z ? (void *)this : 0;} inline ld norm2() const {return s * s + x * x + y * y + z * z;} inline ll intnorm2() const {return llroundl(norm2());} inline quaternion inv() const {return conj() * (1.l / norm2());} inline quaternion round() const { quaternion A(roundl(s), roundl(x), roundl(y), roundl(z)), B(floorl(s) + .5l, floorl(x) + .5l, floorl(y) + .5l, floorl(z) + .5l ); return (*this - A).norm2() <= (*this - B).norm2() ? A : B; } inline friend quaternion Ldiv(const quaternion &A, const quaternion &B) {return (B.inv() * A).round();} inline friend quaternion Rdiv(const quaternion &A, const quaternion &B) {return (A * B.inv()).round();} inline friend quaternion Lmod(const quaternion &A, const quaternion &B) {return A - B * Ldiv(A, B);} inline friend quaternion Rmod(const quaternion &A, const quaternion &B) {return A - Rdiv(A, B) * B;} // find (Hurwitz) quaternion d where A = d u, B = d v. inline friend quaternion Lgcd(quaternion A, quaternion B) { static quaternion tmp; for (; B; tmp = Lmod(A, B), A = B, B = tmp); return A; } inline bool is_int() const {return fabsl(s - roundl(s)) < 1e-10l && fabsl(x - roundl(x)) < 1e-10l && fabsl(y - roundl(y)) < 1e-10l && fabsl(z - roundl(z)) < 1e-10l;} }; const int N = 10010; int R; int sqfree[N]; ll coe[N][4]; void test(quaternion q) { if (q.intnorm2() > R) return; if (q == quaternion()) return; quaternion unit, res; unit = quaternion(0, 1, 0, 0); res = q * unit * q.conj(); int a = std::abs(res.x) + std::abs(res.y) + std::abs(res.z); unit = quaternion(0, 0, 1, 0); res = q * unit * q.conj(); int b = std::abs(res.x) + std::abs(res.y) + std::abs(res.z); unit = quaternion(0, 0, 0, 1); res = q * unit * q.conj(); int c = std::abs(res.x) + std::abs(res.y) + std::abs(res.z); for (int i = 1; i <= (R - 1) / std::max({a, b, c}); ++i) if (sqfree[i]) { int u = a * i, v = b * i, w = c * i, p = std::max({u, v, w}) + 1; coe[p][3] += 1; coe[p][2] -= u + v + w; coe[p][1] += u * (ll)v + v * (ll)w + w * (ll)u; coe[p][0] -= u * (ll)v * w; } } std::ostream& operator<<(std::ostream& os, __int128 x) { if (x / 10) os << (x / 10); return os << int(x % 10); } int main() { R = 10000; // std::cin >> R; std::fill(sqfree + 1, sqfree + R + 1, 1); for (int i = 2; i <= R; i += 2) sqfree[i] = 0; for (int i = 2; i * i <= R; ++i) for (int j = i * i; j <= R; j += i * i) sqfree[j] = 0; int rt = ceil(sqrt(R)) + 1; for (int s = -rt; s <= rt; ++s) for (int x = -rt; x <= rt; ++x) for (int y = -rt; y <= rt; ++y) for (int z = -rt; z <= rt; ++z) { test(quaternion(s, x, y, z)); test(quaternion(s + .5, x + .5, y + .5, z + .5)); } for (int i = 1; i <= R; ++i) { __int128 ans = 0; for (int j = 3; j >= 0; --j) ans = (ans * i + (coe[i][j] += coe[i - 1][j])); std::cout << i << " " << ans / 24 << '\n'; } return 0; }