问题
如果p是一个直角三角形的周长，设这个直角三角形的三条{a,b,c}边都是整数，当p=120时，有三个解：
{20,48,52}, {24,45,51}, {30,40,50}
当 p ≤ 1000 时，解的个数最多的p值是多少？
* 传送门
分析
这个问题如果想暴力解题，借助《问题9》对整数直角三角形的代码，从p=3至p=1000遍历一遍就可以得到结果。
然而，当p继续增大，这个暴力算法就会显得十分慢。接近\(O(n^3)\)的时间复杂度的确算不上什么好算法。
观察\(a^2 + b^2 = c^2, a + b + c = p \leq 1000, c > b, c > a, a+b > c\)，将\(c = p – a – b\)代入平方式，得\(a^2 + b^2 = (p – a – b)^2 = p^2 + a^2 + b^2 – 2ap – 2bp + 2ab\)，两边消项后得\(p^2 – 2ap – 2bp + 2ab = 0\)。
观察a的值，有\(a = \frac{p(p – 2b)}{2(p – b)}\)。由于2b肯定是偶数，如果p是奇数，(p-2b)也肯定是奇数，分子两个奇数相乘，再除以分母，其中分母有2，显然p如果是奇数，a不可能有整数解，所以第一个结论就是p肯定是偶数。
然后设b是三条边中最短的边，那肯定有\(b \leq \lfloor \frac{p}{3} \rfloor\)，因为直角三角形的最短边肯定小于平均边长，即b小于三分之一的周长。只要将b从1到\(\lfloor \frac{p}{3} \rfloor\)的整数遍历一次，便可以知道p(p-2b)是否能被2(p-b)整除，即便知道当前的b对应的a值是否有整数解。a、b同时具有整数解的情况下，便是我们所需要的。时间复杂度会小于\(O(n^2)\)，显然那比暴力算法快一个量级。
最后，这个问题还可以用勾股数公式解。
设\(m > n, a = 2mn \cdot t, b = (m^2 – n^2) \cdot t, c = (m^2 + n^2) \cdot t\)，那么便有\(a + b + c = p = t \cdot (2mn + 2m^2) = t \cdot 2m(n + m) \leq 1000\)。
由于p是偶数，所以\(\frac{p}{2} = t \cdot m \cdot (m + n) \leq 500\)肯定有解。所以\(\frac{p}{2}\)必然分别能被t、m和(m+n)整除。
同样地，如果将上式改写为\(\frac{p}{2} = t \cdot m^2 \cdot (1 + \frac{n}{m}) \leq 500\)，移项可得\(m^2 = \frac{p}{2t(1+\frac{n}{m})}\)，即\(m = \sqrt{\frac{p}{2t(1+\frac{n}{m})}}\)。由于\(m > n, t \geq 1\)，所以\(\frac{n}{m} < 1, 1+\frac{n}{m} > 1, t(1+\frac{n}{m}) > 1, m^2 = \frac{p}{2t(1+\frac{n}{m})} < \frac{p}{2}\)，即\(1 \leq m \leq \lfloor \sqrt{\frac{p}{2}} \rfloor\)。
再考虑(m+n)的性质。由于 m > n，所以\((m+n) < 2m\)。又由于\(t \cdot (m+n) = \frac{p}{2m}, t \geq 1\)，所以\((m+n) \leq \frac{p}{2m}\)。
又由于\(m(m+n) \cdot t = \frac{p}{2}\)，t本来就是为了使勾股数公式能推广的公因子，所以如果等式能够成立，m必然和(m+n)互质（不再有公因子了、一旦有便会被算入t）。它们互质，还可以推导出奇偶性。如果m和(m+n)都是偶数，可以想像得出，它们必不互质，所以不存在这种情况；能满足条件的，就是要么一奇一偶、要么都是奇数，即它们至少有一个是奇数。由于m的范围前面已经推导出来，没有明显的奇偶性，所以我们可以只考虑(m+n)是奇数的这种单边情况。
所以(m+n)是一个大于m的奇数，且小于2m，小于等于\(\frac{p}{2m}\)。
基于勾股数公式的推理实际上并不简洁，执行效率和简单地判断 p(p-2b) 与 2(p-b) 之间的整除性没什么差别，或者说还略显逊色。
答案
// answer1: brute force
const start = Date.now()
function count_tripplet(p) {
let count = 0
for (let c = Math.ceil(p / 2) - 1; c > 0; --c) // 最大的c肯定没有总和的一半大
for (let b = c - 1; b > 0; --b) { // b至少比c少1
const a = p - b - c
if (a < b && (a+b) > c && a*a + b*b == c*c)
count++
}
return count
}
let max_count = 0
let result = 0
for (let p = 1000; p > 3; --p) {
const count = count_tripplet(p)
if (count > max_count) {
max_count = count
result = p
}
}
alert(result + ', ' + (Date.now() - start) + 'ms')
// answer2: after analysis
const start = Date.now()
function count_tripplet(p) {
let count = 0
const max_b = Math.floor(p/3)
for (let b = 1; b < max_b; ++b)
if (p*(p-2*b) % (2*(p-b)) == 0) { // 能整除，代表a有整数解
const a = p*(p-2*b) / (2*(p-b))
if (a > b) // 只计算单向情况，b是最短边，避免重复计算
++count
}
return count
}
let max_count = 0
let result = 0
for (let p = 1000; p > 3; p -= 2) {
const count = count_tripplet(p)
if (count > max_count) {
max_count = count
result = p
}
}
alert(result + ', ' + (Date.now() - start) + 'ms')
// answer3: using formula
const start = Date.now()
function gcd(a, b) { // 尾递归求最大公约数
if (b == 0)
return a
else
return gcd(b, a % b)
}
function count_tripplet(p) {
let count = 0
const p_div_2 = p >> 1
const max_m = Math.floor(Math.sqrt(p_div_2))
for (let m = 1; m <= max_m; ++m)
if (p_div_2 % m == 0) {
const _2m = 2*m
const p_div_2m = p_div_2/m
// (m+n)是大于m的奇数
for (let m_plus_n = (m%2 == 0) ? m+1 : m+2; m_plus_n < _2m && m_plus_n <= p_div_2m; m_plus_n += 2)
if (p_div_2m % m_plus_n == 0 && gcd(m_plus_n, m) == 1) // (m+n), m 互质
count++
}
return count
}
let max_count = 0
let result = 0
for (let p = 1000; p > 3; p -= 2) {
const count = count_tripplet(p)
if (count > max_count) {
max_count = count
result = p
}
}
alert(result + ', ' + (Date.now() - start) + 'ms')
由 bruce 于 2018-04-16 20:43 更新
