假设x为奇数，y为偶数，则z为奇数，2z与2x的最大公因数为2，2z和2x可分别写作
- 2z = (z + x) + (z - x)
- 2x = (z + x) - (z - x)
那么跟据最大公因数性质，z + x和z - x的最大公因数也为2，又因为：
- (z + x)(z - x) = y2，两边同除以4得：
((z + x) / 2)((z - x) / 2) = (y / 2)2
故可令：
- z + x = 2m2, z - x = 2n2
其中z = m + n, x = m - n（m与n互质）
则有：
- y2 = z2 - x2 = 2m22n2 = 4m2n2
即y = 2mn。
综上所述，可得到下式：
- x = m2 - n2, y = 2mn, z = m2 + n2. (m, n为任意自然数)
这里还有一个问题：题目要求统计(x, y, z)三元组的数量时只统计x，y和z两两互质的的情况，这个问题用上面的算法就可以解决了。但对于统计p的数量，题目并不限定三元组是两两互质的。但是上式不能生成所有x, y, z并不是两两互质的情况。然而假设x与y最大公因数w不为1，则z也必能被w整除，因此w为x, y, z三个数的公因数。归纳总结可知，所有非两两互质的x0, y0, z0都可由一组互质的x, y, z乘以系数得到。根据以上理论就可以快速的求解了。
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define MAX 1000010
using namespace std;
bool vis[MAX];
int prime[MAX][],n;
int gcd(int a, int b)
{
return b == ? a : gcd(b, a%b);
}
int upper(int l, int r, int v)
{
int m;
while (l < r)
{
m = l + (r - l) / ;
if (prime[m][] <= v) l = m + ;
else r = m;
}
return r;
}
int cmp(const void*a, const void*b)
{
return ((int*)a)[] - ((int*)b)[];
}
int main()
{
int i,j,z,x,y,k=;
int imax = int(sqrt(MAX >> )+0.5),jmax;
for (i = ; i <= imax; i++)
{
jmax = int(sqrt(MAX - i*i) + 0.5);
for (j = i + ; j <= jmax; j++)
if ((i & ) + (j & )== && gcd(i, j) == )//(i&1)+(j&1)==1 一奇一偶
{
y = * i*j;
z = i*i + j*j;
x = j*j - i*i;
if (x*x+y*y==z*z&&z<=)
{
prime[k][] = x;
prime[k][] = y;
prime[k++][] = z;
}
}
}
qsort(prime,k,sizeof(prime[]),cmp);
while (scanf("%d", &n) == )
{
int a = upper(, k, n);
memset(vis, , n + );
for (i = ; i < a; i++)
for (j = ; j*prime[i][] <= n; j++)
{
vis[j*prime[i][]] = ;
vis[j*prime[i][]] = ;
vis[j*prime[i][]] = ;
}
int count = ;
for (i = ; i <= n; i++)
if (!vis[i]) count++;
printf("%d %d\n", a, count);
}
return ;
}