UVA106 - Fermat vs. Pythagoras

假设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\n", a, count);
65     }
66     return 0;
67 }
View Code

 

posted @ 2015-09-30 15:55  cdongyang  阅读(228)  评论(0)  编辑  收藏  举报