# Is this Euler-Mascheroni constant calculation from double integrals a true identity?

A prime number is a number that is only divisible by itself and one, that is the number of divisors of a prime number is equal to $2$. One way to illustrate this is to plot a matrix such that if the column index (1,2,3,...) divides the row index (1,2,3,...) then a black square is drawn at that column and row index intersection, like this:

As a matrix this has the definition in Mathematica:

MatrixForm[Table[Table[If[Mod[n, k] == 0, 1, 0 ], {k, 1, 12}], {n, 1, 12}]]

which gives the table below:

$$\small \left( \begin{array}{cccccccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 1 & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{array} \right)$$

The number of divisors of n, also called tau, is the row sums of this matrix and starts: 1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6 ...

So at a row equal to a prime number there are only two black dots and the row sums in the above matrix are equal to 2.

To find the ordinary generating function of the number of divisors of n we consider table $a$ times $b$, or row index times column index:

MatrixForm[Table[Table[a*b, {b, 1, 12}], {a, 1, 12}]]


$$\small \left( \begin{array}{cccccccccccc} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 \\ 2 & 4 & 6 & 8 & 10 & 12 & 14 & 16 & 18 & 20 & 22 & 24 \\ 3 & 6 & 9 & 12 & 15 & 18 & 21 & 24 & 27 & 30 & 33 & 36 \\ 4 & 8 & 12 & 16 & 20 & 24 & 28 & 32 & 36 & 40 & 44 & 48 \\ 5 & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 & 50 & 55 & 60 \\ 6 & 12 & 18 & 24 & 30 & 36 & 42 & 48 & 54 & 60 & 66 & 72 \\ 7 & 14 & 21 & 28 & 35 & 42 & 49 & 56 & 63 & 70 & 77 & 84 \\ 8 & 16 & 24 & 32 & 40 & 48 & 56 & 64 & 72 & 80 & 88 & 96 \\ 9 & 18 & 27 & 36 & 45 & 54 & 63 & 72 & 81 & 90 & 99 & 108 \\ 10 & 20 & 30 & 40 & 50 & 60 & 70 & 80 & 90 & 100 & 110 & 120 \\ 11 & 22 & 33 & 44 & 55 & 66 & 77 & 88 & 99 & 110 & 121 & 132 \\ 12 & 24 & 36 & 48 & 60 & 72 & 84 & 96 & 108 & 120 & 132 & 144 \end{array} \right)$$

By stretching out this table in the downward direction we get this table:

MatrixForm[Table[Table[If[Mod[a, b] == 0, a, 0], {b, 1, 12}], {a, 1, 12}]]


$$\small \left( \begin{array}{cccccccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 3 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 4 & 4 & 0 & 4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 5 & 0 & 0 & 0 & 5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 6 & 6 & 6 & 0 & 0 & 6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 7 & 0 & 0 & 0 & 0 & 0 & 7 & 0 & 0 & 0 & 0 & 0 \\ 8 & 8 & 0 & 8 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 \\ 9 & 0 & 9 & 0 & 0 & 0 & 0 & 0 & 9 & 0 & 0 & 0 \\ 10 & 10 & 0 & 0 & 10 & 0 & 0 & 0 & 0 & 10 & 0 & 0 \\ 11 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 11 & 0 \\ 12 & 12 & 12 & 12 & 0 & 12 & 0 & 0 & 0 & 0 & 0 & 12 \end{array} \right)$$

This looks very much like the divisorplot at the beginning. To find the generating function for the number of divisors of $n$ we then introduce the variable $x$ like this:

Clear[x]
MatrixForm[
Table[Table[If[Mod[a, b] == 0, x^a, 0], {b, 1, 12}], {a, 1, 12}]]


$$\small \left( \begin{array}{cccccccccccc} x & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^2 & x^2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^3 & 0 & x^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^4 & x^4 & 0 & x^4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^5 & 0 & 0 & 0 & x^5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^6 & x^6 & x^6 & 0 & 0 & x^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^7 & 0 & 0 & 0 & 0 & 0 & x^7 & 0 & 0 & 0 & 0 & 0 \\ x^8 & x^8 & 0 & x^8 & 0 & 0 & 0 & x^8 & 0 & 0 & 0 & 0 \\ x^9 & 0 & x^9 & 0 & 0 & 0 & 0 & 0 & x^9 & 0 & 0 & 0 \\ x^{10} & x^{10} & 0 & 0 & x^{10} & 0 & 0 & 0 & 0 & x^{10} & 0 & 0 \\ x^{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & x^{11} & 0 \\ x^{12} & x^{12} & x^{12} & x^{12} & 0 & x^{12} & 0 & 0 & 0 & 0 & 0 & x^{12} \end{array} \right)$$

The generating function for the number of divisors of $n$ is then the infinite double sum which is the sum of the infinite matrix above:

$$\sum\limits_{n=1}^{\infty} d(n)x^{n} = \sum\limits_{a=1}^{\infty} \sum\limits_{b=1}^{\infty} x^{a b}$$ $$= 1\, x^{1} + 2\, x^{2} + 2\, x^{3} + 3\, x^{4} + 2\, x^{5} + 4\, x^{6} + 2\, x^{7} + 4\, x^{8} + 3\, x^{9} + 4\, x^{10} + 2\, x^{11} + 6\, x^{12} + ...$$

where $d(n)$ = the number of divisors of $n$

A sum $\sum$ is a sum of distances - usually under a function, while an integral $\int$ is a sum of areas under a function, and a double integral $\int \int$ is a sum of volumes and so on... Often in mathematics one wants to compare a sum to an integral, or at least so I have heard. In this case we have a double sum:

$$\sum\limits_{a=1}^{\infty} \sum\limits_{b=1}^{\infty} x^{a b}$$

Considering the double integral instead we have:

$$\int\limits_{a=1}^{\infty} \int\limits_{b=1}^{\infty} x^{a b}$$

As suggested by Chris's wise sister in the chat room about two weeks ago, the similar triple integral evaluates faster at $x=\frac{1}{e}$:

According to Mathematica:

Integrate[Integrate[1/Exp[1]^(a*b), {a, 1, Infinity}], {b, 1, Infinity}]


the double integral then evaluates to:

$$\int\limits_{a=1}^{\infty} \int\limits_{b=1}^{\infty} \frac{1}{e^{a b}}=-\text{Ei}(-1)$$

Or as Mathematica calls it, the exponential integral. Reading about the exponential integral in Wikipedia that passing a logarithm as argument to it one gets the logarithmic integral and remembering that the logarithmic integral has to do with the prime number theorem, I did some experimenting in Mathematica with the program:

Monitor[MatrixForm[
Table[Table[
Integrate[
Integrate[Exp[a*b], {a, k, Infinity},
GenerateConditions -> False], {b, 0, n},
GenerateConditions -> False] -
Integrate[
Integrate[Exp[a*b], {a, 0, n}, GenerateConditions -> False], {b,
k, Infinity}, GenerateConditions -> False], {n, 1, 4}], {k, 0,
4}]], k]


and came to the conclusion that for $n \gt 0$ and $k \geq 0$, this appears to be an identity:

$$\gamma = \int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db$$

where $\gamma$ = Euler-Mascheroni constant, approximately equal to 0.57721566490153286060651209...

My question is: Why is this identity true?

I also noted that:

For $k \lt 0$ and $n \lt 0$, it appears to be true that:

$$\gamma = \Re(\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db)$$

$$-\pi = \Im(\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db)$$

But this is probably less interesting.

As a starting point I believe it has something to do with the fact:

-Integrate[Integrate[Exp[a*b], {a, 0, 1}, GenerateConditions -> False], {b, 0,
Infinity}, GenerateConditions -> False]


$$\gamma = -\int_0^{\infty } \left(\int_0^1 e^{a b} \, da\right) \, db$$

and that it somehow can be separated out from both of the double integrals: $$\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db$$ $$\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db$$

depending on the values of $n$ and $k$ since:

Monitor[MatrixForm[
Table[Table[
Integrate[
Integrate[Exp[a*b], {a, k, Infinity},
GenerateConditions -> False], {b, 0, n},
GenerateConditions -> False] -
0*Integrate[
Integrate[Exp[a*b], {a, 0, n},
GenerateConditions -> False], {b, k, Infinity},
GenerateConditions -> False], {n, 1, 4}], {k, 0, 4}]], k]
Monitor[MatrixForm[
Table[Table[
0*Integrate[
Integrate[Exp[a*b], {a, k, Infinity},
GenerateConditions -> False], {b, 0, n},
GenerateConditions -> False] -
Integrate[
Integrate[Exp[a*b], {a, 0, n}, GenerateConditions -> False], {b,
k, Infinity}, GenerateConditions -> False], {n, 1, 4}], {k, 0,
4}]], k]


for $n \geq 1$ and $k \geq 0$ we have for the first double integral:

$$\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db$$

(for $n=1,2,3,4...$ and $k=0,1,2,3,4...$):

$$\small \left( \begin{array}{cccc} 0 & -\log (2) & -\log (3) & -\log (4) \\ \gamma -\text{Ei}(1) & \gamma -\text{Ei}(2) & \gamma -\text{Ei}(3) & \gamma -\text{Ei}(4) \\ -\text{Ei}(2)+\gamma +\log (2) & -\text{Ei}(4)+\gamma +\log (2) & -\text{Ei}(6)+\gamma +\log (2) & -\text{Ei}(8)+\gamma +\log (2) \\ -\text{Ei}(3)+\gamma +\log (3) & -\text{Ei}(6)+\gamma +\log (3) & -\text{Ei}(9)+\gamma +\log (3) & -\text{Ei}(12)+\gamma +\log (3) \\ -\text{Ei}(4)+\gamma +\log (4) & -\text{Ei}(8)+\gamma +\log (4) & -\text{Ei}(12)+\gamma +\log (4) & -\text{Ei}(16)+\gamma +\log (4) \end{array} \right)$$

and for the second double integral:

$$-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db)$$

(for $n=1,2,3,4...$ and $k=0,1,2,3,4...$):

$$\small \left( \begin{array}{cccc} \gamma & \gamma +\log (2) & \gamma +\log (3) & \gamma +\log (4) \\ \text{Ei}(1) & \text{Ei}(2) & \text{Ei}(3) & \text{Ei}(4) \\ \text{Ei}(2)-\log (2) & \text{Ei}(4)-\log (2) & \text{Ei}(6)-\log (2) & \text{Ei}(8)-\log (2) \\ \text{Ei}(3)-\log (3) & \text{Ei}(6)-\log (3) & \text{Ei}(9)-\log (3) & \text{Ei}(12)-\log (3) \\ \text{Ei}(4)-\log (4) & \text{Ei}(8)-\log (4) & \text{Ei}(12)-\log (4) & \text{Ei}(16)-\log (4) \end{array} \right)$$

I was not entirely clear in the use of integration indices $n$ and $k$. Usually $n$ is the row index and $k$ is the column index, while in the matrices above they have been permuted.

Edit 10.6.2013: I should have used the 'PrincipalValue' command instead of the 'GenerateConditions' command, like this:

t = 1;
Table[Integrate[
Integrate[Exp[a*b] + 1/b/t, {a, 0, t}], {b, -Infinity, n},
PrincipalValue -> True], {n, 1, 4}]
Table[Integrate[
Integrate[Exp[a*b] + 1/b/t, {a, 0, t}], {b, -Infinity, Log[n]},
PrincipalValue -> True], {n, 1, 4}]


$$\text{Ei}(n \cdot t) = \int_{-\infty }^n \left(\int_0^t \left(e^{a b}+\frac{1}{b \cdot t}\right) \, da\right) \, db$$

And the relation from Wikipedia I wanted explore is: $$\text{Ei}(\log (n)) = \text{li}(n)$$

Edit 11.6.2013:

This allows us to write:

X=5;
Monitor[Table[
Integrate[
Integrate[Exp[a*b] + n/b, {a, 0, 1/n}], {b, -Infinity, Log[X]},
PrincipalValue -> True], {n, 1, 6}], n]


$$\text{li}(X^{\frac{1}{n}}) = \int_{-\infty }^{\text{Log}[X]} \left(\int_0^{\frac{1}{n}} \left(e^{a b}+\frac{n}{b}\right) \, da\right) \, db$$

Which simplifies to:

X = 5;
Monitor[Table[
Integrate[E^(b/n)/b, {b, -Infinity, Log[X]},
PrincipalValue -> True], {n, 1, 6}], n]


$$\text{li}(X^{\frac{1}{n}}) = \int_{-\infty }^{\text{Log}[X]} \frac{e^{b/n}}{b} \, db$$

• Nice question (+1) Jun 9, 2013 at 15:33
• It looks like Mathematica is evaluating the integrals with different conditions. To me it seems that the integrals should be the same and should not converge. Jun 10, 2013 at 2:35
• Actually when I look at the two double integral now once more, I also think they are the same. Jun 10, 2013 at 3:12
• It was the use of Mathematica's 'GenerateConditions' that mislead me. What I should have used was 'PrincipalValue': t = 2; Table[Integrate[ Integrate[Exp[ab] + 1/b/t, {a, 0, t}], {b, -Infinity, n}, PrincipalValue -> True], {n, 1, 12}] Table[Integrate[ Integrate[Exp[ab] + 1/b/t, {a, 0, t}], {b, -Infinity, Log[n]}, PrincipalValue -> True], {n, 1, 12}] Setting t=1 gives the logarithmic integral. Jun 10, 2013 at 15:41

What justifies diligently the step $\sum\limits_{a=1}^{\infty} \sum\limits_{b=1}^{\infty} x^{a b} \rightarrow \int\limits_{a=1}^{\infty} \int\limits_{b=1}^{\infty} x^{a b}$; this seems likely very vague "what is valid for distances ... valid for areas ... valid for volumes..." so might one argument and extend to higher dimmensions $\int\int\int\int...$ So far cannot follow to justify identity.
$$\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db=-\int^1_0\int^1_0\frac{1-a}{(1-ab)\log (ab)} da\; db = \gamma$$
• There is no need to justify the step between the double sum and the double integral. What I am asking is if the following is true: $$\gamma = \int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db$$ where $\gamma$ is the Euler-Mascheroni constant = 0.57721566490153286060651209... Jun 9, 2013 at 15:17