假设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乘以系数得到。根据以上理论就可以快速的求解了。
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<algorithm> 5 #define MAX 1000010 6 using namespace std; 7 bool vis[MAX]; 8 int prime[MAX][3],n; 9 int gcd(int a, int b) 10 { 11 return b == 0 ? a : gcd(b, a%b); 12 } 13 int upper(int l, int r, int v) 14 { 15 int m; 16 while (l < r) 17 { 18 m = l + (r - l) / 2; 19 if (prime[m][2] <= v) l = m + 1; 20 else r = m; 21 } 22 return r; 23 } 24 int cmp(const void*a, const void*b) 25 { 26 return ((int*)a)[2] - ((int*)b)[2]; 27 } 28 int main() 29 { 30 int i,j,z,x,y,k=0; 31 int imax = int(sqrt(MAX >> 1)+0.5),jmax; 32 for (i = 1; i <= imax; i++) 33 { 34 jmax = int(sqrt(MAX - i*i) + 0.5); 35 for (j = i + 1; j <= jmax; j++) 36 if ((i & 1) + (j & 1 )== 1 && gcd(i, j) == 1)//(i&1)+(j&1)==1 一奇一偶 37 { 38 y = 2 * i*j; 39 z = i*i + j*j; 40 x = j*j - i*i; 41 if (x*x+y*y==z*z&&z<=1000000) 42 { 43 prime[k][0] = x; 44 prime[k][1] = y; 45 prime[k++][2] = z; 46 } 47 } 48 } 49 qsort(prime,k,sizeof(prime[0]),cmp); 50 while (scanf("%d", &n) == 1) 51 { 52 int a = upper(0, k, n); 53 memset(vis, 0, n + 1); 54 for (i = 0; i < a; i++) 55 for (j = 1; j*prime[i][2] <= n; j++) 56 { 57 vis[j*prime[i][0]] = 1; 58 vis[j*prime[i][1]] = 1; 59 vis[j*prime[i][2]] = 1; 60 } 61 int count = 0; 62 for (i = 1; i <= n; i++) 63 if (!vis[i]) count++; 64 printf("%d %d ", a, count); 65 } 66 return 0; 67 }