How to optimize my code for finding the count of all integral medians for all possible integral triangles with a <= b <= c <= 100000?

Question

I am doing the solution for this problem from 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).

Analyse:

1. a <= b <= c <= n == 100000
2. ABC is a triangle so it should: abs(a-b) < c < a+b
3. Mc = sqrt(2 * a^2+ 2 * b^2 - c^2) / 2 wikipedia
4. Mc is integer so 2 * a^2+ 2 * b^2 - c^2 should be a perfect square and divisible by 4.

Code:

#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);

}

Issues:

It works fine for small n but the complexity of the solution is too high, I assume it is O(n^3)(am I wrong?) which will take days for n = 100000.

How could I improve this whether with a mathematical or algorithmic way?

Updates

I got those suggestions:

• Calculating power of a outside the b/c loops and power of b outside c loop. This improved slightly the performance.
• c cannot be odd. then a and b must have same parity. This improved the performance 4 times.
• Using threads to divide the work on many cores. It may improve by a factor close to number of cores.
• A mathematical solution posted in math.stackexchange. It claims O(N^5/2) for a basic solution and can achieve O(N^2) by using O(N^2) of memory. I didn't test it yet.

Answer 1

Since this is a Project Euler problem, you are supposed to be able to do it in about a minute of computing time on a modern computer. They don't always stick to that, but it indicates that a running time of k*n^2 or k*n^2*log(n) is probably fine if the constant isn't too bad, but probably not k*n^2.5 or k*n^3.

As SleuthEye commented, the side c must be even or else the inside of the square root would have to be odd, so taking the square root and dividing by 2 could not make an integer.

You can simplify the equation to 4(mc^2+(c/2)^2) = 2(a^2+b^2).

Here is one approach: Create two dictionaries, left and right. For each, let the keys be possible values of that side of the equation, and let the values be a list of the pairs (mc,c/2) or (a,b) which produce the key. For the right dictionary, we only need to consider where a and b have the same parity, and where 1<=a<=b<=n. For the left, we only need to consider 1<=c/2<=n/2 and 1<=mc<=sqrt(3)/2 n since 4mc^2 = 2a^2+2b^2-c^2 <= 3b^2 <=3n^2.

Then go through the possible keys, and compare the elements of the values from each dictionary, finding the number of compatible ((mc,c/2),(a,b)) pairs where b <= c < a+b. This inner step is not constant time, but the maximum and average lengths of the lists are not too long. The ways to write n as a sum of two squares roughly correspond (up to units) to the ways to factor n in the Gaussian integers, and just as the largest number of factors of an integer does not grow too rapidly, the same is true in the Gaussian integers. This step takes O(n^epsilon) time for any epsilon>0. So, the total running time is O(n^(2+epsilon)) for any epsilon>0.

In practice, if you run out of memory, you can construct partial dictionaries where the keys are restricted to be in particular ranges. This parallelizes well, too.