欧拉计划问题# 276 - 原始三角形答案

【问题标题】：Project Euler Problem# 276 - Primitive Triangles欧拉计划问题# 276 - 原始三角形
【发布时间】：2010-02-14 11:53:11
【问题描述】：

我试图解决problem# 276 from Project Euler，但到目前为止没有成功。

考虑整数的三角形边 a、b 和 c，a ≤ b ≤ c。一个整数边三角形 (a,b,c) 是如果 gcd(a,b,c)=1，则称为原语。如何许多原始整数边三角形存在周长不超过 10 000 000？

我的代码中的瓶颈是GCD 函数。我的样本输入几乎需要 90% 的执行时间。我很想听听有关如何改进解决方案的提示和 cmets……我的解决方案是用 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;
}

【问题讨论】：

标签： c performance algorithm math

【解决方案1】：

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

您目前仅单独计算原始三角形。相反，问问自己这个问题：你能有效地计算周长为 a+b+c=n 的所有个三角形（不一定是原始三角形）吗？（运行时间 O(n) 就可以了——你当前的算法是 Ω(n³)。）这样做之后，你多算了哪些三角形？例如，有多少个三角形用 p 分割边？（提示：a+b+c=n (a/p) + (b/p) + (c/p) = n/p。）还有so on。

编辑：在解决了这个问题并检查了 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 且所有三个具有相同的奇偶性（与 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) 相关联，以获得递归关系。

（当然，也有一些非常简单的公式，你用谷歌就能找到，但那有什么好玩的呢？）

【讨论】：

【解决方案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。

不幸的是，即使有这些更改以及迄今为止其他答案中提到的更改，这可能也太慢了。我相信真正的解决方案不会枚举所有可能的三角形，而是以某种方式将它们分组计数。

【讨论】：

【解决方案3】：

蛮力解决方案往往真的很慢:)

减少计算的提示：

预先计算所有素数范围 2..10⁷ 并存储它们在查找表中。
在bi_gcd(a, b) 检查 a 或 b 是否是素数并且返回1 否则计算他们的 gcd。
在tri_gcd(a, b, c) 中执行相同操作。

编辑：考虑到@interjay 的评论：

如果例如a 是质数，我们仍然需要检查b 是否是a 的倍数，
像这样(b % a == 0)。在那种情况下，a 是 gcd。

【讨论】：

5 是素数，但 gcd(5,10) 不是 1。您仍然需要检查另一个数是否是素数的倍数。
b%a==0 比(b/a)*a==b 好得多。但随后它由欧几里得步骤处理。

【解决方案4】：

除了 Nick D 的回答（当一个或两个输入为素数时，它会阻止您计算 bi_gcd）之外，我发现自己想知道您必须调用多少次 bi_gcd 与 相同（复合) 数字。无论您计算多少次，GCD(12, 18) 始终为 6。我怀疑，存储结果的机制可能会有所改善。

【讨论】：

顺便说一句，这种技术称为“记忆化”。

【解决方案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-s，我可以将花费的时间从 1.2 秒减少到 1.0 秒（使用 gcc -O3）。

【讨论】：