c - 如何优化我的代码以查找所有可能的积分三角形的所有积分中位数，其中 a <= b <= c <= 100000？

Question

我正在从 Euler Project Problem 513, Integral median 解决这个问题:

ABC is an integral sided triangle with sides a≤b≤c. mc is the median connecting C and the midpoint of AB. F(n) is the number of such triangles with c≤n for which mc has integral length as well. F(10)=3 and F(50)=165.

Find F(100000).

分析:

1. a <= b <= c <= n == 100000
2. ABC 是一个三角形，所以它应该:abs(a-b) < c < a+b
3. Mc = sqrt(2 * a^2+ 2 * b^2 - c^2) / 2 wikipedia
4. Mc是整数所以 2 * a^2+ 2 * b^2 - c^2应该是一个完美的正方形并且可以被 4 整除。

代码:

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

#define N 100000
#define MAX(a,b) (((a)>(b))?(a):(b))

void main(){
    unsigned long int count = 0;
    unsigned long int a,b,c;
    double mc;

    for (a = 1; a <= N; a++) {
        printf("%lu\n", a);
        for (b = a; b <= N; b++)
            for (c = MAX(b, abs(b-a)); c <=N && c < a+b; c++){
                mc = sqrt(2 *a *a + 2 * b * b - c * c)/2.0;
                if (mc-(unsigned long)mc == 0)
                    count++;
            }
    }
     printf("\ncpt == %lu\n", count);

}

问题:

它适用于小型 n但是解决方案的复杂度太高了，我假设是O(n^3) (我错了吗？)n = 100000 需要几天时间.

我怎样才能通过数学或算法的方式改进它？

更新

我得到了那些建议:

• a 的计算能力外面b/c b 的循环和功率外面c环形。这略微提高了性能。
• c不能奇怪。然后 a和 b必须具有相同的奇偶校验。这将性能提高了 4 倍。
• 使用线程将工作分配到多个核心上。它可能会提高接近核心数量的一个因素。
• 在 math.stackexchange 中发布的数学解决方案.它声称 O(N^5/2)一个基本的解决方案，可以实现O(N^2)通过使用 O(N^2)的内存。我还没有测试。

Answer 1

由于这是一个欧拉计划问题，您应该能够在现代计算机上用大约一分钟的计算时间来完成它。他们并不总是坚持这一点，但这表明如果常量不是太差，运行时间 k*n^2 或 k*n^2*log(n) 可能没问题，但可能不是 k*n^2.5 或 k*n^3 。

正如 SleuthEye 评论的那样，边 c 必须是偶数，否则平方根的内部必须是奇数，因此取平方根除以 2 不能得到整数。

您可以将等式简化为 4(mc^2+(c/2)^2) = 2(a^2+b^2) 。

这是一种方法:创建左右两个字典。对于每个，让键为等式那一侧的可能值，并让值成为产生键的 (mc,c/2) 或 (a,b) 对的列表。对于正确的字典，我们只需要考虑 a 和 b 的奇偶性相同的地方，以及 1<=a<=b<=n 的地方。对于左边，我们只需要考虑 1<=c/2<=n/2 和 1<=mc<=sqrt(3)/2 n 以来的 4mc^2 = 2a^2+2b^2-c^2 <= 3b^2 <=3n^2 。

然后遍历可能的键，并比较每个字典中值的元素，找到兼容的 ((mc,c/2),(a,b)) 对的数量，其中 b <= c < a+b 。这个内部步骤不是常数时间，但列表的最大和平均长度不会太长。将 n 写为两个平方和的方法大致对应(最多单位)高斯整数中将 n 分解的方法，正如最大的 number of factors of an integer 不会增长太快一样，高斯整数也是如此。此步骤对于任何 O(n^epsilon) 都需要 epsilon>0 时间。因此，对于任何 O(n^(2+epsilon)) ，总运行时间为 epsilon>0 。

在实践中，如果内存不足，您可以构建部分字典，其中的键被限制在特定范围内。这也可以很好地并行化。