Как оптимизировать мой код для нахождения счета всех интегральных медианов для всех возможных интегральных треугольников с <= b <= c <= 100000?

Я делаю решение этой проблемы из Задача Эйлера Project 513, Интегральная медиана:

ABC является интегральным односторонним треугольником со сторонами a≤b≤c. mc - медиана соединяющий C и середину AB. F (n) - это число таких треугольники с c≤n, для которых mc также имеет целую длину. F (10) = 3 и F (50) = 165.

Найти F (100000).

Анализ:

a <= b <= c <= n == 100000
ABC - это треугольник, поэтому он должен: abs(a-b) < c < a+b
Mc = sqrt(2 * a^2+ 2 * b^2 - c^2) / 2 wikipedia
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) памяти. Я еще не тестировал его.

Ответ 1

Так как это проблема Project Euler, вы должны иметь возможность сделать это примерно за минуту вычислительного времени на современном компьютере. Они не всегда придерживаются этого, но это указывает на то, что время работы 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 в гауссовских целых числах и точно так же, как наибольшее число факторов целого числа не растет слишком быстро, то же самое верно и в гауссовских целых числах. Этот шаг занимает O(n^epsilon) время для любого epsilon>0. Таким образом, общее время работы O(n^(2+epsilon)) для любого epsilon>0.

На практике, если у вас закончилась нехватка памяти, вы можете создать частичные словари, где ключи ограничены в определенных диапазонах. Это тоже хорошо распараллеливается.