# A double series with divisibility restrictions

How may I restrict $i,j$ such that they run over $\mathbb{N}$ excepting the numbers divisible by $2$ or $3$ or $7$?

Sum[1/( i j (i + j + 1)), {i, 1, Infinity}, {j, 1, Infinity}]

• Could you please add a bit of context just to say how you got into this problem ? Commented Sep 13, 2013 at 15:20
• @b.gatessucks I'm doing some personal research on multiseries. Commented Sep 13, 2013 at 15:32

Could do explicit inclusion/exclusion as below (I forget what the correct term is but you'll get the idea). First do one index.

s1 = Sum[1/(i j (i + j + 1)), {j, 1, Infinity}];

s2 = Sum[1/(i j (i + j + 1)), {j, 2, Infinity, 2}];
s3 = Sum[1/(i j (i + j + 1)), {j, 3, Infinity, 3}];
s7 = Sum[1/(i j (i + j + 1)), {j, 7, Infinity, 7}];

s23 = Sum[1/(i j (i + j + 1)), {j, 2*3, Infinity, 2*3}];
s27 = Sum[1/(i j (i + j + 1)), {j, 2*7, Infinity, 2*7}];
s37 = Sum[1/(i j (i + j + 1)), {j, 3*7, Infinity, 3*7}];
s237 = Sum[1/(i j (i + j + 1)), {j, 2*3*7, Infinity, 2*3*7}];

stot = s1 - (s2 + s3 + s7) + (s23 + s27 + s37) - s237;


Can now repeat the process with the other index. Below would be the first step.

t1 = Sum[stot, {i, 1, Infinity}]


Unfortunately I cannot get a closed form for this.

If there is a way to rewrite as some form of convolution sum that might be amenable to Fourier series type methods. Offhand I don't see such a rewrite but I could easily be missing something of that sort.

• OK. Thank you for answer! Commented Sep 13, 2013 at 15:11

I tried to go the usual tricky route in order to use symmetries :

$$\sum_{i,j} \frac{1}{i\ j (i+j+1)} = \int_{0}^{\infty} dx\ e^{-x} \left(\sum_{i} \frac{e^{-i x}}{i} \right)^2$$

The sum can be performed in closed form :

s = Sum[Exp[-i x]/i Boole[Mod[i, 2] != 0] Boole[Mod[i, 3] != 0] Boole[
Mod[i, 7] != 0] , {i, 1, Infinity}]


but the integral cannot. Numerically :

NIntegrate[Exp[-x] s s, {x, 0, Infinity}]
(* 0.480345 *)


which looks in line with @tintin.

You can also use Boole as follows:

Sum[1/(i j (i + j + 1)) Boole[! Divisible[i, 2] && !
Divisible[i, 3] && ! Divisible[i, 7] && ! Divisible[j, 2] && !
Divisible[j, 3] && ! Divisible[j, 7]], {i, 1, 1000}, {j, 1,
1000}] // N


0.478819

Which gives the same numerical value as tintin's Mod version.

To speed things up, you can use the new iteration form in Mathematica

First define lis as follows:

lis = Select[Range[1, 1000, 2], Not@Divisible[#, 3] && Not@Divisible[#, 7] &] // N;


Then:

Sum[1/(i j (i + j + 1)), {i, lis}, {j, lis}]


0.478819

• OK. Thank you for answer! Commented Sep 13, 2013 at 15:42
• 2x speedup: add N@ before Select Commented Sep 13, 2013 at 17:52
• @ybeltukov. Yes, good point. And this makes the // N after Sum redundant. Thanks. Commented Sep 13, 2013 at 18:10
f[i_, j_] := 0;
f[i_, j_] :=
1/(i j (i + j + 1)) /;
Mod[i, 2] != 0 && Mod[i, 3] != 0 && Mod[i, 7] != 0 &&
Mod[j, 2] != 0 && Mod[j, 3] != 0 && Mod[j, 7] != 0


you can sum it with

Sum[f[i,j],{i,1,1000},{j,1,1000}]


which I get the number of 0.478819 I write 1000 because I cant deal with the infinity thing

• note you can do Sum[f[i, j], {i, 1, 1000, 2}, {j, 1, 1000, 2}] and omit the Mod[,2] check.. speed it up a tad. Commented Sep 13, 2013 at 15:13
• @george2079,that is definately a big improvement Commented Sep 13, 2013 at 15:19

I think you can achieve this way too,

Total[Total[
Table[If[(Divisible[i, #] && Divisible[j, #]) & /@ {2, 3, 7}  /.
List -> Or, 1/(i j (i + j + 1)), 0], {i, 1, 1000}, {j, 1,
1000}]]] // N


0.344854

Sum will achieve the same result, its gives same result as a Table with two iterators as in your case. If you do NSum[1/(i j (i + j + 1)), {i, 1, 1000}, {j, 1, 1000}] you get

SequenceLimit::seqlim: The general form of the sequence could not be determined, and the result may be incorrect. 1.97715

But if you use Total[Total[Table[1/(i j (i + j + 1)), {i, 1, 1000}, {j, 1, 1000}]] // N]

1.98443 (No warnings)

Now to see that same results are generated by Sum and custom loop check these,

Sum[1/(i j (i + j + 1)), {i, 1, 50}, {j, 1, 50}] // N


1.81095

Total[Total[Table[1/(i j (i + j + 1)), {i, 1, 50}, {j, 1, 50}]] // N]


1.81095

'NSum' probably uses approximations hence its result is pretty close but might be wrong if high precision needed.