# Can this simple loop with several counters be improved or parallelized

I am running the following calculations. The timing for H = 2000 is about 8 days on my machine. Can this be parallelized or improved upon in some other way. Also I am trying to avoid storing the results of each each step, this example with H = 2000 already has some 32 billion cases.

S2IntegerS2A2ThreeParameterCounts[H_] :=
Timing[Print["Start time:  ", Date[]];
Block[{Ncases, NS2, NA2, Nreducible, Nirreducible},
Ncases = 0;
NS2 = 0;
NA2 = 0;
Nreducible = 0;
Nirreducible = 0;
Monitor[Do[Ncases++; If[GCD[r, s, t] == 1,
If[IntegerQ[Sqrt[s^2 - 4 r t]], NA2++; Nreducible++, NS2++;
Nirreducible++]],
{r, 1, +H}, {s, -H, +H}, {t, -H, +H}],
ProgressIndicator[Ncases, {0, (H (2*H + 1)^2)}]];
Print["Height H =  ", H, "  where  ", -H,
" \[LessEqual] a1, a0 \[LessEqual] +", H];
Print["Number of initial test cases:  ", Ncases];
Print["Number of reducible cases:  ", Nreducible];
Print["Number of irreducible cases:  ", Nirreducible];
Print["Number of unique test cases:  ", Nreducible + Nirreducible]];
Print["End time:  ", Date[]]]

• Would you take a few minutes to time about how long it takes for each increment in r? r==1 won't be fair, but timing each of the next few increments should give a baseline. Then new experiment when you nest two Do loops, the outer one with r and the inner one with s and t. Eliminate the Monitor and print your diagnostic info after the inner Do finishes. Doing this will try to estimate how much overhead your Monitor is taking and still give you a printout, perhaps every 10 minutes or so. It looks like your GCD==1 is only eliminating about 17% of the cases so not a huge amount to save there. – Bill Apr 16 '16 at 6:57

Baseline

As a baseline for computational results and speed, use

H = 100; Ncases = H (2 H + 1)^2; Nreducible = 0; Nirreducible = 0;
Do[If[GCD[r, s, t] == 1, If[IntegerQ[Sqrt[s^2 - 4 r t]], Nreducible++, Nirreducible++]],
{r, 1, +H}, {s, -H, +H}, {t, -H, +H}] // AbsoluteTiming
{Ncases, Nreducible, Nirreducible, Nreducible + Nirreducible}


which yields

(* {4040100, 76728, 3278393, 3355121} *)


in about 29 sec. Increasing H by a factor of two increase run-time by factor of eight, as it should, so the Do overhead is negligible for H this large.

Improved Algorithm

Since s enters the argument of Sqrt, the sign of s does not matter. So, it makes sense to perform the computation for s > 0, double its output, and add to it the output of the computation for s = 0. Additionally, Sqrt need be computed only if its argument is non-negative. With these changes,

H = 100; Ncases = H (2 H + 1)^2; Nreducible = 0; Nirreducible = 0;
(Do[If[GCD[r, s, t] == 1, If[s^2 - 4 r t >= 0 && IntegerQ[Sqrt[s^2 - 4 r t]],
Nreducible++, Nirreducible++]], {r, 1, +H}, {s, 1, +H}, {t, -H, H}];
Nreducible = 2 Nreducible; Nirreducible = 2 Nirreducible;
Do[If[GCD[r, t] == 1, If[IntegerQ[Sqrt[- r t]], Nreducible++, Nirreducible++]],
{r, 1, +H}, {t, -H, +H}]) // AbsoluteTiming
{Ncases, Nreducible, Nirreducible, Nreducible + Nirreducible}


which yields the same result in about 12 sec.

Parallelize Code

There are a variety of ways to run this code in parallel. For instance,

f[h1_, h2_] := Module[{nred = 0, nirred = 0}, Do[If[GCD[r, s, t] == 1,
If[IntegerQ[s^2 - 4 r t >= 0 && Sqrt[s^2 - 4 r t]], nred++, nirred++]],
{r, h1, h2}, {s, 1, +H}, {t, -H, H}]; {nred, nirred}]

H = 100; n = 4; Ncases = H (2 H + 1)^2;
({Nreducible, Nirreducible} = 2 Total[ParallelTable[f[1 + (i - 1) H/n, i H/n],
{i, n}, Method -> "FinestGrained"]];
Do[If[GCD[r, t] == 1, If[IntegerQ[Sqrt[- r t]], Nreducible++, Nirreducible++]],
{r, 1, +H}, {t, -H, +H}]) // AbsoluteTiming
{Ncases, Nreducible, Nirreducible, Nreducible + Nirreducible}


which yields the same result in a bit over 4 sec. Multiple factors seem to prevent reaching the 3 sec theoretical limit. For instance, on my 4-processor PC the CPU utilization per kernel is lower by almost 15% with ParallelTable as compared with Table. Also, the four parallel calculations do not all take the same amount of time. Note that H/n must be an integer for the code above to produce correct results.

In combination these two approaches reduce total run-time by almost a factor of seven.

• This was very informative and helpful my series of calculations. I am getting a performance gain of about 6.5 in my case. I was looking for the method to split the calculation between the processors. – Lorenz H Menke Apr 18 '16 at 21:31