讨论 $n$ 的奇偶性，利用通项 $(n^2-1)^2+(2n)^2=(n^2+1)^2$，枚举 $\gcd(a,b,c)$ 即可。

#include <bits/stdc++.h>
using namespace std;
int c[1000012], l[1000012], r[1000012], d[20000012];
int e[1000012], f = 0, g = 0;
struct A
{
	int a;
	int b;
	int c;
};
A a[10000012];
bool cmp(A x, A y)
{
	if (x.a != y.a) return x.a < y.a;
	if (x.b != y.b) return x.b < y.b;
	return x.c < y.c;
}
bool operator!=(A x, A y)
{
	return (x.a != y.a || x.b != y.b || x.c != y.c);
}
signed main()
{
	ios::sync_with_stdio(0);
	int n;
	cin >> n;
	for (int i = 1; i <= n; i++)
		for (int j = i; j <= n; j += i)
			c[j]++;
	for (int i = 1; i <= n; i++)
	{
		r[i] = r[i - 1] + c[i];
		l[i] = r[i - 1] + 1;
	}
	memset(c, 0, sizeof c);
	for (int i = 1; i <= n; i++)
		for (int j = i; j <= n; j += i)
		{
			c[j]++;
			d[l[j] + c[j] - 1] = i;
		}
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= f; j++)
			e[j] = 0;
		f = 0;
		for (int j = l[i]; j <= r[i]; j++)
			for (int k = j; k <= l[i] + r[i] - j; k++)
			{
				f++;
				e[f] = d[j] * d[k];
			}
		for (int j = 1; j <= f; j++)
		{
			long long u = 1ll * i * i / e[j];
			if (u - e[j] > 2 * i && (u - e[j]) % 2 == 0 && u + e[j] <= 2 * n)
			{
				g++;
				a[g] = (A) {i, (u - e[j]) / 2, (u + e[j]) / 2};
			}
		}
	}
	sort(a + 1, a + g + 1, cmp);
	for (int i = 1; i <= g; i++)
		if (a[i] != a[i - 1]) printf("%d %d %d\n", a[i].a, a[i].b, a[i].c);
	return 0;
}