#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;
}