欧拉计划问题＃276 - 原始三角形

Question

我尝试解决 Project Euler的问题＃276，但迄今为止没有成功。

考虑具有整数边a、b和c（a≤b≤c）的三角形。如果gcd（a，b，c）= 1，则称为整数边三角形（a，b，c）是原始的。存在多少个原始整数边三角形的周长不超过10,000,000？

我的代码瓶颈在于GCD函数。对于我的样本输入，它几乎占据了执行时间的90％。 我很想听到有关如何改进解决方案的提示和评论...我的解决方案用C编写，但非常简单：

#include <stdio.h>
#include <math.h>

#define sides 3
// This is where my program spends most of its time
size_t bi_gcd (size_t a, size_t b);
size_t tri_gcd (size_t a, size_t b, size_t c);
size_t count_primitive_triangles (size_t maximum_perimeter);

int main()
{
 printf( "primitive_triangles = %lu \n",
            count_primitive_triangles(1000) );
}

size_t bi_gcd (size_t a, size_t b)
{
 size_t t;

 while(b != 0)
 {
  t = b;
  b = a % b;
  a = t;
 }

 return a;
}

size_t tri_gcd (size_t a, size_t b, size_t c)
{
 return bi_gcd(a, bi_gcd(b, c));
}

size_t count_primitive_triangles (size_t max_perimeter)
{
 size_t count = 0; // number of primitive triangles
 size_t a, b, c;   // sides of the triangle
 // the following are the bounds of each side
 size_t 
  // because b >= a && c >= b ==> max of a
  // is when a+b+c > 10,000,000
  // b == a (at least)
  // c == b (at least)
  // ==> a+b+c >= 10,000,000
  // ==> 3*a   >= 10,000,000
  // ==> a= 10,000,000/3
  a_limit = max_perimeter/sides,
  b_limit,
  c_limit;

 for (a = 1; a <= a_limit; ++a)
 {
  // because c >= b && a+b+c <= 10,000,000
  // ==> 2*b + a = 10,000,000
  // ==> 2*b = 10,000,000 - a
  // ==> b = (10,000,000 - a)/2
  for (b = a, b_limit = (max_perimeter-a)/2; b <= b_limit; ++b)
  {
   // The triangle inequality:
   // a+b > c (for a triangle to be valid!)
   // ==> c < a+b
   for (c = b, c_limit = a+b; c < c_limit; ++c)
   {
    if (tri_gcd(a, b, c) == 1)
     ++count;
   }
  }
 }

 return count;
}

- Khaled Alshaya

Answer 1

很好，你在优化之前进行了剖析，但是在gcd函数中花费大量时间并不一定意味着你需要（或可以）使其更快，而是因为你调用它太频繁了。:-) 这里有一个提示——一种算法改进，可以将运行时间提高数个数量级，而不仅仅是实现方面的改进。

你当前只是单独计算原始三角形的数量。相反，问问自己：是否可以有效地计算出周长为a+b+c=n的所有三角形（不一定是原始三角形）？（运行时间O(n)就可以——你当前的算法是Ω(n^3)）。做到这一点之后，你又重复计算了哪些三角形？例如，有多少个以p为边长的三角形？（提示：a+b+c=n <=> (a/p) + (b/p) + (c/p) = n/p.）等等。 编辑：在解决问题并查看Project Euler线程后，我发现还有其他一些不错的方法来解决这个问题，但上述方法是最常见的，并且它有效。对于第一部分，您可以直接计数（有些人已经这样做了；当然可行），或者您可能会发现这个额外的提示/技巧有帮助：

• 假设 1 ≤ a ≤ b ≤ c，且 a+b+c=n。让 p = b+c-a = n-2a，q = c+a-b = n-2b，r = a+b-c = n-2c。那么 1 ≤ r ≤ q ≤ p，p+q+r=a+b+c=n。而且，p、q 和 r 的奇偶性与 n 相同。
• 反之，对于任何满足 1 ≤ r ≤ q ≤ p 且 p+q+r=n（三者的奇偶性相同）的情况， 令 a = (q+r)/2，b = (r+p)/2，c = (p+q)/2。 那么 1 ≤ a ≤ b ≤ c 且 a+b+c=n。此外，它们是整数并形成一个三角形（请检查）。
• 因此，周长为 n 的整数三角形（a,b,c）的数量正好等于n中同奇偶数部分(p,q,r)的组合数。您可能会发现这更容易计算。

此外/或者，尝试将此数字 T(n) 直接与较小的数字 T(n-k) 相关联，以获得递归关系。

当然，你可以通过谷歌搜索找到一些非常简单的公式，但这样做有什么乐趣呢？

- ShreevatsaR

Answer 2

以下是您可以改进的一些内容：

• 在内部循环中，不要计算tri_gcd(a,b,c)，而是在第二个循环中计算g1 = gcd(a,b)，在内部循环中计算g2 = gcd(g1, c)。

• 当gcd(a,b)==1时，可以避免内部循环，并将计数器增加max_c - min_c + 1，因为您知道对于任何c的值，gcd都为1。

• 您的内部循环似乎计数过高，因为它没有检查a+b+c <= 10000000。

不幸的是，即使进行这些更改和其他答案中提到的更改，这可能仍然太慢了。我认为真正的解决方案不是枚举所有可能的三角形，而是以某种方式按组计数。

- interjay

Answer 3

暴力解法通常非常缓慢 :)

减少计算量的提示：

1. 预先计算范围在2..10⁷内的所有质数，并将它们存储在查找表中。
2. 在 bi_gcd(a, b) 中检查 a 或 b 是否为质数，如果是，则返回 1，否则计算它们的gcd。
   在 tri_gcd(a, b, c) 中也要做同样的操作。

编辑：参考 @interjay 的评论:

例如，如果 a 是一个质数，我们仍然需要检查 b 是否是 a 的倍数，
像这样 (b % a == 0)。在这种情况下，a 就是最大公约数。

- Nick Dandoulakis

---

5是质数，但gcd(5,10)不为1。你仍然需要检查另一个数字是否是该质数的倍数。 - interjay

“b%a==0”比“(b/a)*a==b”好得多。但是它会被欧几里得步骤处理。 - kennytm

Answer 4

除了Nick D的答案，当一个或两个输入是质数时，它将阻止您计算bi_gcd。我想知道您必须调用多少次bi_gcd与相同（复合）数字。无论您计算多少次，GCD（12,18）始终为6。存储结果的机制可能会改善情况，我猜测。

- AakashM

---

顺便说一下，这种技术被称为“记忆化”。 - Svante

Answer 5

作为一个小的改进，你可以插入一些特殊情况，例如：

。

size_t bi_gcd (size_t a, size_t b) {
    if (a == 1 || b == 1) return 1;
    ...

size_t tri_gcd (size_t a, size_t b, size_t c) {
    if (a%2==0 && b%2==0 && c%2==0) return 2;
    // of course GCD(a,b,c) can be > 2,
    // but the exact GCD is not needed in this problem.
    ...

有了这两个if，我可以将执行时间从1.2秒减少到1.0秒（使用gcc -O3）。

- kennytm