Sum over n variables

You can write a few helper functions to help you. The following can probably be streamlined...

vars[s_String, n_Integer?Positive] := Table[Symbol[s <> ToString[i]], {i, 1, n}]
vars[sym_Symbol, num_] := vars[SymbolName[sym], num]

nestedRange[vars_List, min_, max_] /; min <= max := 
 Transpose@{vars, ConstantArray[min, Length[vars]], Append[Rest@vars, max]}

nestedSum[f_, vars:{__Symbol}, min_, max_] /; min <= max := 
 With[{r = Sequence @@ Reverse@nestedRange[vars, min, max]}, Sum[f, r]]
nestedSum[f_, {var:(_String|_Symbol), num_Integer?Positive}, 
 min_, max_] /; min <= max := nestedSum[f, vars[var, num], min, max]

Then, for example

nestedSum[f[a, b, c], {a, b, c}, 0, Infinity] // TraditionalForm

produces

A larger sum is

In[]:= nestedSum[Total@vars[i, 4], {i, 4}, 1, 20] // Timing
Out[]= {0.016001, 371910}

which can be compared with evaluating the same thing using Boole

In[]:= v = vars[i, 4];
       With[{r = Sequence@@Table[{n, 1, 20}, {n, v}]}, 
            Sum[Boole[LessEqual@@v] Total@v, r]] // Timing
Out[]= {0.056003, 371910}

Here's another approach. With this solution you would enter the sum in the form sum[f, a0<a1<a2<=...<an], where you can use a mix of < and <= in the specification of the range of summation.

sum[f_, Less[b0_, a__, b1_]] := Sum[f, ##] & @@ 
  Reverse[MapIndexed[{#1[[1]], b0 + #2[[1]], #1[[2]] - 1} &, 
    Partition[{a, b1}, 2, 1]]]

sum[f_, LessEqual[b0_, a__, b1_]] := Sum[f, ##] & @@ 
  Reverse[MapThread[{#1, b0, #2} &, {{a}, Rest[{a, b1}]}]]

sum[f_, (ii : HoldPattern[Inequality[PatternSequence[_, 
   (Less | LessEqual)] .., _]])] := 
  Module[{lst, ops, vars, b0, b1, count},
   lst = List @@ ii;
   ops = lst[[2 ;; -1 ;; 2]] /. {Less -> 1, LessEqual -> 0};
   vars = lst[[;; -1 ;; 2]];
   Sum[f, ##] & @@ 
     Reverse[MapThread[{#1[[1]], lst[[1]] + #2, #1[[2]] - #3} &, 
       {Partition[Rest[vars], 2, 1], Accumulate[Most[ops]], Rest[ops]}
     ]
  ]]

So for example

sum[f[a, b, c], 0 < a <= b < c < d] // TraditionalForm

produces

I will tackle the sum explicitly mentioned in your post.

Suppose you use i[1], i[2] and so on for iteration variables. The product can be formed as follows. -Subsets[Array[i, n], {2}] will form pairs. $i_s-i_r$ for $1\leqslant r<s\leqslant n$.

f[n_Integer?Positive] := Block[{i},
 Sum[Apply[Times, -Subtract @@@ Subsets[Array[i, n], {2}]], ##] & @@ 
  Reverse@Partition[Array[i, n + 1, 1] /. {i[n + 1] -> 2 n + 1}, 2, 1]]

The sequence $f_n$ grows rather fast:

In[79]:= {f[1], f[2], f[3], f[4], f[5], f[6]}

Out[79]= {3, 20, 588, 114048, 194347296, 3634840535040}