A double series with divisibility restrictions

Question

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}]

Answer 1

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.

Answer 2

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.

Answer 3

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

Answer 4

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

Answer 5

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.