问题
如果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 更新