UVA 106 - Fermat vs. Pythagoras

題目敘述

輸入 N (1 ≤ N ≤ 1000000)，求X² + Y² = Z²（0 < X < Y < Z ≤ N）的解中

三數互質 gcd(X, Y, Z) = 1 的解有幾組
１ ~ Ｎ中有幾個正整數沒被任何一組解包含

解法

原本想枚舉 Y 和 Z，但時間上會是 10¹¹ 等級，顯然不是一個好方法。

這時，突然想到了「畢式三元數」(Pythagorean triple)，畢竟每次數學培訓上數論都會講到，前幾天翻數學競賽的書也剛好看到，不過似乎只針對「互質畢式三元數」討論而已。

「互質畢式三元數」性質如下：

X, Y, Z ∈ 𝚴，X² + Y² ＝ Z² 且 gcd(X, Y, Z) = 1，則可表示為

X = m² - n²

Y = 2mn

Z = m² + n²

其中 m > n，gcd(m, n) = 1，m, n 為一奇一偶。

這題的範圍不僅僅限於 gcd(X, Y, Z) = 1，但只要將每組「畢式三元數」除以 gcd(X, Y, Z)，就都能變成「互質畢式三元數」了。換句話說，「畢式三元數」一定是由某一組「互質畢式三元數」乘ｋ倍而來。表示為 X = k (m² - n²)，Y = 2kmn，Z = k (m² + n²)。

考慮 m, n 的範圍，1000000 ≥ N ≥ Z ＝ m² + n²，得 1000 ≥ √N ≥ m > n ≥ 1。枚舉 m, n的時間為 10⁶ 等級，明顯能快很多。

我的方法是跑一層 for 枚舉 m，裡面再跑一層 for 枚舉 n，if(gcd(m, n) == 1)，則跑一層 k。

寫完後，丟到 UVA 上排名 21，時間 0.050。（第一名 0.008，我超好奇怎麼辦到？）天啊！這真的怪爆了，我到底哪裡比別人快啊？於是 google 了別人的解法。

看到大部分的人都是對每筆輸入 N 跑一次 for m in range(1,sqrt(N))，這方法我沒想到……因為我寫 DP 寫中毒了，於是動不動就想直接給他一次找完，也就是直接「建表」。

這題的數論難關解決後，建表的過程確實有點卡，於是就來記一下重點部分。

大部分人是直接開個 bool 紀錄數字出現過了沒，還有開個變數紀錄有幾組互質的。

首先，來談互質的組數 p。用 DP 的角度切入，p[i] 和 p[i-1] 有關嗎？

p[i] = p[i - 1] ＋ (z == i 的互質組數)。所以就開個 cnt 陣列紀錄，每找到一組互質畢式三元數，cnt[Z]++。p[N]就是cnt[N]的前綴和，O(N)一次做完即可。

接著，是算有幾個沒用過的數。1 ~ i + 1 會比 1 ~ i 多一個數能用，而 1 ~ i 用過的數 1 ~ i + 1 也同樣被用過。num[i + 1] ＝ num[i] + 1 - (有幾個數在 Z == i + 1 時才出現)。所以重點是計算有幾個數在 Z = i + 1 時第一次出現，也就是紀錄每個數搭配的 Z 的最小值，最後在掃過一遍就行了。

因為都是由互質三元數去找其他三元數，所以不會找到重複的，所以效率佳。寫這篇時發現自己沒有應用到一奇一偶的性質，改成 n += 2 後，丟到UVA上排名 13，時間 0.040。

Code

 1#include<bits/stdc++.h>
 2using namespace std;
 3#define MAX 1000001
 4int f[MAX], num[MAX], p[MAX], ans[MAX] = {0, 1};
 5inline bool coprime(int x, int y, int z) {
 6    return __gcd(__gcd(x, y), z) == 1;
 7}
 8int main() {
 9    fill(f, f + MAX, MAX);
10    for (int m = 2; m < 1000; m++)
11        for (int n = (m & 1) + 1; n < m; n += 2) {
12            if (__gcd(m, n) != 1) continue;
13            int x = m * m - n * n, y = 2 * m * n, z = m * m + n * n;
14            if (z >= MAX) break;
15            if (coprime(x, y, z)) {
16                p[z]++;
17                for(int k = MAX / z; k >= 1; k--) {
18                    int xx = x * k, yy = y * k, zz = z * k;
19                    f[xx] = min(f[xx], zz);
20                    f[yy] = min(f[yy], zz);
21                    f[zz] = zz;
22                }
23            }
24        }
25    for(int i = 1; i <= MAX; i++)
26        p[i] += p[i - 1], num[f[i]]++;
27    for (int i = 2; i <= MAX; i++)
28        ans[i] = ans[i - 1] + 1 - num[i];
29    int N;
30    while (scanf("%d", &N) == 1) printf("%d %d\n", p[N], ans[N]);
31    return 0;
32}