Computing a double integral with applications to prime numbers I was reading the preprint [1] which contains on p. 7 the following formula (for $4<s\le6$):
$$
f_1(s)=\frac{2e^\gamma}{s}\left\{\log(s-1)+\int_4^s\int_3^t\frac{\log(u-2)}{u-1}du\,dt \right\}
$$
and which concludes, on page 16, with
$$
c=\frac{2e^\gamma}{3f_1(6)}=1.18751\ldots.
$$
But evaluating the function numerically in PARI/GP
f1(s)=2*exp(Euler)/s*(log(s-1)+intnum(t=4,s,intnum(u=3,t,log(u-2)/(u-1))))
\p 200
2*exp(Euler)/3/f1(6)

gives me c = 0.8239599331..., a very different result. I tried it in Maxima
assume(t >= 4);
float(integrate(integrate(log(u-2)/(u-1),u,3,t),t,4,6));

and it gives 0.81786458079115, which corresponds to c = 1.4518071082, not very close to either result. (bfloat doesn't seem to work here -- it has trouble with the polylogarithms it seems.)
How can I reliable compute this double integral?
[1] Pin-Hung Kao, Almost-Prime Polynomials with Prime Arguments, arXiv:1606.03505 [math.NT]
 A: It seems the authors might have used (I admit that the numbers do not match completely)
$$
c=\frac{2e^\gamma}{3F_1(6)}\approx1.18787
$$
(where $F_1$ is also defined on page 7) instead of
$$
c=\frac{2e^\gamma}{3f_1(6)}\approx 0.82396
$$
I don't have time or energy to understand which one is the correct one to put in the end. Mathematica calculates
$$
F_1(s)=\frac{1}{s}\bigl(2e^\gamma(1+\frac{\pi^2}{12}+\log(s-2)\log(s-1)+\text{Li}_2(2-s)\bigr)
$$
and 
$$
\begin{aligned}
f_1(s)&=\frac{2}s e^\gamma\Bigl(s-4+\frac{\pi^2}{12}s+\frac{5}{2}(\log 2)^2
+2\log 2-3\log2\log 3+\log(s-1)\\
&\quad +\log(s-2)(2-s+(s-1)\log(s-1))\\
&\quad+\text{Li}_2(-1/2)+\text{Li}_2(1/4)+(s-1)\text{Li}_2(2-s)\Bigr)
\end{aligned}
$$
This might be a mistake, and I suggest that you contact the author and ask.
A: Hint:
Though it's tedious but the double integral is doable by subbing $u=1+x$ then repeating integration by parts several times. You might use the following relation
\begin{equation}
\operatorname{Li}_2(z)=\int_z^0\frac{\log(1-x)}{x}\ dx
\end{equation}
