问题
五角数是由公式\(P_n = \frac{n(3n−1)}{2}\)产生。头十个五角数如下：
1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...
可以看到，\(P_4 + P_7 = 22 + 70 = 92 = P_8\)，即它们的和也是五角数。而它们的差，70 − 22 = 48则不是五角数。
请找出一对五角数，\(P_j 与 P_k\)，使得它们的和与差都是五角数，且它们的差值\(D = |P_k − P_j|\)是最小的。请问D的值是多少？
* 传送门
分析
这个问题有点厉害，因为用暴力解题貌似无解，因为没有明显的上界，直接循环去找答案似乎是没有办法的。于是我先一段一段去找，于是让我蒙对了答案。但是暴力解题还得一段段蒙，这显然不是本题的意图。
于是阅读本题的论坛，发现一位法国老师 Francky 给出了完整的证明和解题代码。他原话说得很强力：你们的解题方式都不对！于是我也只能乖乖受教。
首先要从五角数的通项公式推导出相邻两项的差：
\(\begin{aligned}P_{n+1} – P_n &= \frac{(n+1)[3(n+1) − 1]}{2} – \frac{n(3n−1)}{2} \\
&= \frac{(n+1)[3(n+1) − 1] – n(3n − 1)}{2} \\
&= \frac{3n^2 + 6n + 3 – n – 1 – 3n^2 + n}{2} \\
&= \frac{6n + 2}{2} \\
&= 3n + 1
\end{aligned}\)
也就是说，第n项和第n+1项之间有这个关系：\(P_{n+1} = P_n + 3n + 1\)。
这个关系有什么用处？我们可以用它推广出第n项和第n+i项之间的关系（i > 0）：
\(\begin{aligned}P_{n+i} = P_{n+i-1 + 1} &= P_{n+i-1} + 3(n+i-1) + 1 \\
&= P_{n+i-2 + 1} + 3(n+i-1) + 1 \\
&= [P_{n+i-2} + 3(n+i-2) + 1] + 3(n+i-1) + 1 \\
&= P_{n+i-2} + 3[(n+i-2) + (n+i-1)] + 2 \\
&= …
\end{aligned}\)
不断地展开，会发生什么？这个等式右边的第1项不断地展开，到\(P_n\)就不需要再继续了，因为我们已经到达要求第n项和第n+i项之间的关系。而展开的期间，第2项中括号内部则不断地累积变成从n至n+i-1的这i项的和，即\(\frac{i[n+(n+i-1)]}{2} = \frac{i(2n+i-1)}{2}\)，第3项则会累积成i。
于是便会有：
\(\begin{aligned}
P_{n+i} &= P_n + 3\frac{i(2n+i-1)}{2} + i \\
&= P_n + \frac{6in + 3i^2 – 3i + 2i}{2} \\
&= P_n + 3in + \frac{i(3i – 1)}{2}
\end{aligned}\)
观察最后一项，形式与\(P_n\)的通项公式一致！也就是说可以将它换成\(P_i\)，即\(P_{n+i} = P_n + 3in + P_i\)。也就是说，第n项与第n+i项之间的差，即\(P_{n+i} – P_n = 3in + P_i\)。
由于五角数是单调递增的，所以当k > j，\(P_k > P_j\)，即求它们两个五角数之间差值的绝对值，我们只需要考虑k > j的情况。
设\(i > 0, j > 0, k > 0, k > j, k = j + i, P_k = P_{j+i}, D = P_k – P_j = P_{j+i} – P_j\)。利用上面的结论，便会有\(D = 3ij + P_i\)。
又假设满足题目条件的是第x项和第y项，有\(P_x = P_k + P_j, P_y = P_k – P_j = D, x > y\)，那么便会有\(x > j+i, P_x = (P_j + D) + P_j = 2P_j + D = j(3j – 1) + D, P_y = D = 3ij + P_i, y > i\)。
由于\(P_x = \frac{x(3x−1)}{2}\)移项后可以得到这个关于x的一元二次方程：\(3x^2 – x – 2P_x = 0\)，且由于x是正整数，x的解即为\(x = \frac{1+\sqrt{1+24P_x}}{6}\)。
整个问题就转换为搜索y和i的问题。y不断增大，同时i小于y，找出能满足这些结论的\(P_x\)，看看x是否有正整数解，就知道当前的y和i是否能满足题目要求。由于\(D = P_y\)，y从小到大搜索，所以找出来的第一个符合条件的D值肯定就是最小的D值。
不过即便这个法国老师的方法能快一点，但实际上也是一个\(O(n^2)\)时间复杂度的解题方法。这个问题碰巧是要最小的D，如果说要第二小的，那这个方案恐怕也不会有多快。
还是应该说，这个题目本身就有点莫名其妙呢？
答案
极慢，慎点！！
// answer1: brute force
const start = Date.now()
const P = (n) => n*(3*n-1)/2
const can_solve = (p) => (1+Math.sqrt(1+24*p)) % 6 == 0
let result = Infinity
// 这个搜索范围是碰运气、每隔1000试出来的……
for (let j = 1; j <= 2000; ++j)
for (let k = 3000; k > j; --k) {
const P_j = P(j)
const P_k = P(k)
if (can_solve(P_j + P_k)) {
const D = P_k - P_j
if (can_solve(D) && result > D)
result = D
}
}
alert(result + ', ' + (Date.now() - start) + 'ms')
经过分析后，可以在比较短的时间内找到答案：
// answer2: after analysis
const start = Date.now()
const P = (n) => n*(3*n-1)/2
function* solve_D() {
for (let y = 1; true; ++y) {
const D = P(y)
for (let i = 1; i < y; ++i) {
const _3ij = D - P(i)
const _3i = 3*i
if (_3ij % _3i == 0) { // j有整数解
const j = _3ij / _3i
const P_x = j * (3 * j - 1) + D
const _6x = 1 + Math.sqrt(1 + 24 * P_x)
if (_6x % 6 == 0) { // x有整数解
const k = j+i
const x = _6x / 6
yield [D, j, k, x, y]
}
}
}
}
}
for (const [D, j, k, x, y] of solve_D()) {
result = D // 第一个D肯定是满足条件的最小的D
console.log('P_' + x + ' = P_' + k + ' + P_' + j, P(x) == P(k) + P(j))
console.log('D = P_' + y + ' = P_' + k + ' - P_' + j, D == P(y), P(y) == P(k) - P(j))
break
}
alert(result + ', ' + (Date.now() - start) + 'ms')
