关于如何快速生成勾股数三元组

#数学 #Python

问题由来

最近在刷 Project Euler，虽然目前刷的都是简单题，还是发现了我在数学方面的很多不足。

比如这个问题：如何生成周长 $p$ 满足 $p <= 1000$ 的所有勾股数三元组?

暴力求解

最容易想到的肯定是暴力求解法, 设三边为 a, b, c, 且 $a < b < c$，简单确定下a的范围

遍历所有的 a , b 和 p, 挨个判断是不是勾股数

暴力代码如下

for a in range(1, 334):
    for b in range(a+1, 1000-334):
        for p in range(1, 1001):
            c = p - a - b
            if c < a or c < b:
                continue
            if c**2 == a**2 + b**2:
                pass # do something

用时 69.73 s

稍微手算一下的做法

根据 $ a^2 + b^2 = c^2 $ 和 $ p = a + b + c $ , 代入消$c$, 可以得到 p, a, b 的等式

$$ b = \frac{p(p - 2a)}{2(p - a)} $$

代码如下

for p in range(1, 1001):
    for a in range(1, 334):
        if (2 * a) >= p:
            continue
        m = p * (p - 2 * a)
        n = 2 * (p - a)
        if m % n == 0:
            b = m // n
            c = p - a - b
            if a < b and b < c:
                pass # do something

用时 0.12 s

对这道题来说, 从 $O(n^3)$ 改成 $O(n^2)$ 就已经可以接受了

欧几里得公式法

从评论中，又学到了一种古老的经典方法，是欧几里得发现的。

$$ a = 2kmn , b = k (m^2 - n^2) , c = k (m^2 + n^2) $$

其中，$a,b,c$ 是三边，$k > 0, m > n > 0$, m, n 互质，奇偶相反

from fractions import gcd
from math import sqrt
for m in range(2, int(sqrt(1000 / 2))):
    for n in range(1, m):
        if gcd(m, n) == 1 and (m - n) % 2:  # No Duplicates
            for p in range(2 * m * (m + n), 1000, 2 * m * (m + n)):
                pass # do something

用时 0.04 s，更快了

还有其他方法

可以去看看维基上的专门讲这个的页面:

传送门

嘻嘻，又水了一篇博客