# An efficient way to evaluate this expression with deeply nested Do loops

Is there an efficient way to evaluate something like this?

Np=10;
xx = 0.;
Do[
Do[
Do[
Do[
Do[
Do[
If[(1. ix - 1. jx)^2 + (1. iy - 1. jy)^2 + (1. iz - 1. jz)^2 ==
0., xx += 0.,
xx += (1. iz - 1. jz)^2/(Sqrt[(1. ix - 1. jx)^2 + (1. iy -
1. jy)^2 + (1. iz - 1. jz)^2])^3
, {jz, Np}],
{iz, Np}],
{jy, Np}],
{iy, Np}],
{jx, Np}],
{ix, Np}]


The problem is not the complexity, but the fact that there are six independent variables and the test in that If. You can imagine what happens if Np is increased. The time it take to evaluate this increases quite a bit.

I tried to evaluate the expression that can become zero separately and put the results in a six way nested Table. The plan was to then Flatten the Table and take out the values equal to 0 and then perform the summation with something like Total on the inverse values. But this way I bump into a memory problem, because if I increase Np then all that huge Table has to be written to memory.

So the question is: Is there a way to efficiently evaluate the sum in my example, for Np=100?

Edit @ciao

It appears that I oversimplified the question too much. i thought that only the denominator had issues.

So the terms look something like this in reality:

$$\frac{(iz-jz)^2}{((ix-jx)^2+(iy-jy)^2+(iz-jz)^2)^{3/2}}$$, where ix, jx, iy, jy, iz, jz go from 1 to N. I didn't think that the nominator should be included in the original question. How does this modify your answer?

• Just a side question - what is the \$ symbol? – jmm Feb 18 '20 at 1:47
• @jmm Probably mistakenly inserted LaTeX delimiters. – Rohit Namjoshi Feb 18 '20 at 2:31
• yes, sorry, I had an attempt to format the code better but I didn't manage at 3 am :D – lucian Feb 18 '20 at 7:14
• "How does this modify your answer?": It breaks it, so I have removed it. I am loath to work on questions that become moving goalposts. You should be able to modify the first function I posted for your needs, assuming that my interpretation of the ambiguous "... terms look something like this..." is correct. I doubt the faster second function I posted can be equivalently modified. – ciao Feb 18 '20 at 17:30
• I am trully sorry. I didn't realize the complexity introduced with the real term. I thought that the only problem was checking for the denominator not being zero. I am working on modifying you first answer to my needs. I really appreciate your effort! You should leave the answer and I will accept it because it does solve the initial problem I posted in a very elegant way! – lucian Feb 18 '20 at 19:14

## 1 Answer

UPDATE: After reposting the below function after modifications from the deleted post to make it work after OP question was changed, I determined a way to utilize ideas from the original deleted post's second, faster function. This is much faster than the second function below for larger arguments (~100 X faster for Np argument of 500) and has less memory pressure. It will comfortably handle Np arguments into four digit range.

sumx[np_] := Module[{t, c2, u, r, v, p, x, c},
t = Tally[(Subtract @@@ Tuples[Range@np, {2}])^2];
p = Tr[x^t[[All, 1]]*t[[All, 2]]];
c = CoefficientList[p, x];
c2 = ListConvolve[c, c, {1, -1}, 0];
r = Pick[Range[0, Length[c2] - 1], u = Unitize[c2], 1];
v = Pick[c2, u, 1];
Tr[Rest[Join @@ Outer[Times, v, t[[All, 2]]]] Rest[Join @@ ConstantArray[t[[All, 1]], Length@r]]/
Rest[Join @@ Outer[Plus, r, t[[All, 1]]]]^N@(3/2)]]


==================================================

Original reposted answer:

Per OP comment request, modified original function to calculate this sum:

sum[np_] :=
Module[{r2 = Tuples[Tally[(Subtract @@@ Tuples[N@Range@np,{2}])^2], {3}]},
Tr[Times @@@ Transpose[{r2[[2 ;;, 3, 1]]/Total[r2[[2 ;;, All, 1]], {2}]^(3/2),
Times @@@ r2[[2 ;;, All, 2]]}]]];


This should be usefully quicker, and reasonable time to execute up to low 100s arguments:

sum//AbsoluteTiming

{2.97555,4.67618*10^7}

• Thank you for your answer! I will try it in a few mins! I do find it a bit hard to follow thou. Can you add a bit more explanations to it? Thank you! – lucian Feb 18 '20 at 7:15
• I forgot to add more info. This is part of a larger expression where this fraction appears two more times. This is the term that creates the problem because I have to make sure i don't end up with 1/0. So the numbers you got are ok. they look weird yes! but they are ok!. I just have trouble understanding the code. I am going through it now. – lucian Feb 18 '20 at 7:31
• @lucian - not much to explain, I think: the intermediate values in the sum are much smaller cardinality when treated as distinct values. I took advantage of the multiplicity of each (which happens to count the number of distinct boxes of integer dimensions X, Y, Z by diagonals - is that part of what you're doing?) and then form the combinations, and applying the multiplicities. This shrinks the number of terms to be operated on by many orders of magnitude, with concomitant speed increase. There's another level of abstraction that should do the same I think, but I've not spent time on it. – ciao Feb 18 '20 at 8:36
• I think there is a problem with you solution, but it does give me some ideas. I tried to break it down and understand it and instead of N@Range@np I replaced it with a generic list {a,b,c} to see the terms better. When that goes through the Tr at the end, I have three types of terms: i) 27/((a - b)^2)^(3/2); ii) 18/((-a + b)^2 + (a - c)^2)^(3/2) and iii) 6/((-a + b)^2 + (-a + c)^2 + (-b + c)^2)^(3/2). Only the last type of term should exist in the sum. if you look in the example I gave, that type of term is the only one that it is added. – lucian Feb 18 '20 at 9:05
• @lucian - Well, using your code, for Np from 1 to 12 (I do not have the patience to run for larger numbers) and only changing your code and mine to use exact results, my method result and your code result match exactly. So I'm not sure what you mean. Please either verify or clarify. In any case, I have a new method that's orders of magnitude faster beyond my current one, but if the results are not what you're after, comment and I'll delete the answer. – ciao Feb 18 '20 at 10:08