Splitting the integral
We can write
\begin{align}
I &\equiv \int \limits_0^\infty x \operatorname{arccot}(x) [1-x \operatorname{arccot}(x)][2+x \operatorname{arccot}(x)] \, \mathrm{d} x \\
&= \int \limits_0^\infty x \operatorname{arccot}(x) [1-x \operatorname{arccot}(x)] \, \mathrm{d} x + \int \limits_0^\infty x \operatorname{arccot}(x) [1-x^2 \operatorname{arccot}^2(x)] \, \mathrm{d} x \\
&= f_1 (1) + f_2(1) \, ,
\end{align}
where for $a>0$ we let
\begin{align}
f_1(a) &\equiv \int \limits_0^\infty x \operatorname{arccot}(a x) [1-x \operatorname{arccot}(x)] \, \mathrm{d} x \, , \\
f_2(a) &\equiv \int \limits_0^\infty x \operatorname{arccot}(a x) [1-x^2 \operatorname{arccot}^2(x)] \, \mathrm{d} x \, . \\
\end{align}
Computation of $f_1(1)$
Differentiating under the integral sign yields
\begin{align}
f_1 '(a) &= -\int \limits_0^\infty \frac{x^2}{1+a^2 x^2} [1-x \operatorname{arccot}(x)] \, \mathrm{d} x \\
&= -\frac{1}{a^2}\left[\int \limits_0^\infty [1-x \operatorname{arccot}(x)] \, \mathrm{d} x - \int \limits_0^\infty \frac{1-x \operatorname{arccot}(x)}{1+a^2 x^2} \, \mathrm{d} x \right] \\
&\equiv -\frac{I_1 + g_1 (a)}{a^2} \, .
\end{align}
We integrate by parts to find
\begin{align}
I_1 &= \lim_{r \to \infty} \left[r-\frac{r^2}{2} \operatorname{arccot}(r) - \frac{1}{2} \int \limits_0^r \frac{x^2}{1+x^2} \, \mathrm{d} x\right] \\
&= \lim_{r \to \infty} \left[\frac{r-r^2 \operatorname{arccot}(r)}{2}+ \frac{1}{2} \int \limits_0^r \frac{\mathrm{d} x}{1+x^2} \right] \\
&= \frac{\pi}{4} \, .
\end{align}
Another integration by parts reduces $g_1(a)$ to
$$g_1(a) = -\frac{\pi}{2 a} + \frac{1}{2a^2} \int \limits_0^\infty \frac{\ln(1+a^2 x^2)}{1+x^2} \, \mathrm{d} x \, . $$
The remaining integral can be evaluated using differentiation under the integral sign and a partial fraction decomposition. We find
$$g_1(a) = -\frac{\pi}{2 a} + \frac{\pi \ln(1+a)}{2a^2} \, . $$
Since $f_1(\infty) = 0$, we can integrate to get
\begin{align}
f_1(1) &= \int \limits_\infty^1 f_1'(a) \, \mathrm{d} a = \int \limits_1^\infty \left[\frac{\pi}{4a^2} - \frac{\pi}{2 a^3} + \frac{\pi}{2 a^4} \ln(1+a)\right] \, \mathrm{d} a \\
&= \frac{\pi}{3}\ln(2) - \frac{\pi}{12} \, .
\end{align}
Computation of $f_2(1)$
In the same manner we can compute
\begin{align}
f_2 '(a) &= -\int \limits_0^\infty \frac{x^2}{1+a^2 x^2} [1-x^2 \operatorname{arccot}^2(x)] \, \mathrm{d} x \\
&= -\frac{1}{a^2}\left[\int \limits_0^\infty [1-x^2 \operatorname{arccot}^2(x)] \, \mathrm{d} x - \int \limits_0^\infty \frac{1-x^2 \operatorname{arccot}^2(x)}{1+a^2 x^2} \, \mathrm{d} x \right] \\
&\equiv -\frac{I_2 + g_2 (a)}{a^2} \, .
\end{align}
We now need to integrate by parts twice to evaluate
\begin{align}
I_2 &= \frac{1}{3} \left[\int \limits_0^\infty \frac{1}{1+x^2} \, \mathrm{d} x + \int \limits_0^\infty \frac{\ln(1+x^2)}{1+x^2} \, \mathrm{d} x + \lim_{r\to\infty} r \{2-r \operatorname{arccot}(r)[1+r\operatorname{arccot}(r)]\}\right] \\
&= \frac{\pi[2 \ln(2)+1]}{6} \, .
\end{align}
We can write
$$ g_2 (a) = - \frac{\pi}{2 a} + \int \limits_0^\infty \frac{x^2 \operatorname{arccot}^2 (x)}{1+a^2 x^2} \, \mathrm{d} x \equiv - \frac{\pi}{2a} + h(a,1) $$
with
$$ h(a,b) = \int \limits_0^\infty \frac{x^2 \operatorname{arccot}^2 (b x)}{1+a^2 x^2} \, \mathrm{d} x $$
for $a,b >0$. Differentiation under the integral sign and partial fractions reduce the derivative of $h(a,b)$ with respect to $b$ to known integrals:
\begin{align}
\partial_b h(a,b) &= - 2 \int \limits_0^\infty \frac{x^3 \operatorname{arccot} (b x)}{(1+a^2 x^2)(1+a^2 x^2)} \, \mathrm{d} x \\
&= \frac{2}{a^2-b^2}\left[\int \limits_0^\infty \frac{x \operatorname{arccot} (b x)}{1+a^2 x^2} \, \mathrm{d} x - \int \limits_0^\infty \frac{x \operatorname{arccot} (b x)}{1+b^2 x^2} \, \mathrm{d} x \right] \\
&= \frac{\pi[\ln(a+b)-\ln(2b)]}{a^2(a^2-b^2)} - \frac{\pi \ln(2)}{a^2 b^2} \, .
\end{align}
After integrating with respect to $a$ and $b$ and evaluating all elementary integrals we are left with
$$f_2(1) = \pi \ln(2) - \frac{\pi}{24} - \pi J \, ,$$
where
\begin{align}
J &= \frac{1}{2} \int \limits_1^\infty \int \limits_1^\infty \frac{b^{-4} \ln(a) - a^{-4} \ln(b)}{a^2-b^2} \, \mathrm{d} a \, \mathrm{d} b \\
&= \frac{1}{2} \int \limits_0^1 \int \limits_0^1\frac{v^4 \ln(u) - u^4 \ln(v)}{u^2-v^2} \, \mathrm{d} u \, \mathrm{d} v \, .
\end{align}
This integral can be evaluated by introducing the similar integral
$$ K = \frac{1}{2} \int \limits_0^1 \int \limits_0^1\frac{v^4 \ln(v) - u^4 \ln(u)}{u^2-v^2} \, \mathrm{d} u \, \mathrm{d} v \, . $$
We have
\begin{align}
K &= - \frac{1}{2} \frac{\mathrm{d}}{\mathrm{d} s} \int \limits_0^1 \int \limits_0^1\frac{u^s - v^s}{u^2-v^2} \, \mathrm{d} u \, \mathrm{d} v ~~ \Bigg \rvert_{s=4} = - \frac{\mathrm{d}}{\mathrm{d} s} \left[\frac{1}{s} \int \limits_0^1 \frac{1-t^s}{1-t^2} \, \mathrm{d} t \right] \\
&= \frac{1}{16} \int \limits_0^1 (1+t^2)\, \mathrm{d} t - \frac{1}{4} \int \limits_0^1 \frac{- \ln(t) t^4}{1-t^2} \, \mathrm{d} t \\
&= \frac{1}{12} - \frac{1}{4} \sum_{k=0}^\infty \frac{1}{(2k+5)^2} = \frac{13}{36} - \frac{\pi^2}{32} \, .
\end{align}
Since
$$ J + K = \frac{1}{2} \int \limits_0^1 \int \limits_0^1 - \ln(u v) (u^2+v^2) \, \mathrm{d} u \, \mathrm{d} v = \frac{4}{9} \, , $$
we conclude that
$$ J = J + K - K = \frac{1}{12} + \frac{\pi^2}{32} \, .$$
Therefore we arrive at
$$ f_2(1) = \pi \ln(2) - \frac{\pi}{8} - \frac{\pi^3}{32} \, . $$
Final result
We now obtain the value
$$ I = f_1(1) + f_2(1) = \frac{4\pi}{3} \ln(2) - \frac{5\pi}{24} - \frac{\pi^3}{32} \, . $$
for the original integral. Note that $I \approx 1.2800035$ is very close but not equal to the rational number $\frac{32}{25} = 1.28$ .
I will assume $s \in \mathbb{N}_0$ to avoid complex numbers. We can let $x = \mathrm{e}^{-t}$ and use $\operatorname{li}(\mathrm{e}^{-t}) = \operatorname{Ei}(-t)$, where $\operatorname{Ei}$ is the exponential integral defined by
$$ \operatorname{Ei}(z) = \int \limits_{-\infty}^z \frac{\mathrm{e}^{y}}{y} \, \mathrm{d} y \, , \, z \in \mathbb{R} \, .$$
The principal value is to be used for $z>0$ . Then we can integrate by parts (using $\frac{\mathrm{d}}{\mathrm{d} t} \frac{t^{s+1} [(s+1)\ln(t) - 1]}{(s+1)^2} = t^s \ln(t)$)
\begin{align}
G(s) &= \int \limits_0^\infty (-1)^s \operatorname{Ei}(-t) t^s \ln(t) \, \mathrm{d} t = \frac{(-1)^s}{(s+1)^2} \int \limits_0^\infty \left[t^s \mathrm{e}^{-t} - (s+1) t^s \ln(t) \mathrm{e}^{-t}\right] \, \mathrm{d} t \\
&= \frac{(-1)^s}{(s+1)^2} [\Gamma(s+1) - (s+1) \Gamma'(s+1)] \, .
\end{align}
Using $\Gamma(s+1) = s!$ and $\Gamma'(s+1) = s! (H_s - \gamma)$ with the harmonic number $H_s$ now yields
$$G(s) = \frac{(-1)^s s!}{(s+1)^2} [ 1 + (s+1)(\gamma-H_s)] \, . $$
Best Answer
Here's a way that only requires the identity $\Gamma'(1) = - \gamma$ , since all derivatives of higher order cancel: \begin{align} I &=\int \limits_0^\infty (x^2 - 3x+1) \ln^3 (x) \mathrm{e}^{-x} \, \mathrm{d} x \\ &= \left[\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 3 \frac{\mathrm{d}}{\mathrm{d}t}+1 \right] \int \limits_0^\infty \ln^3 (x) \mathrm{e}^{- t x} \, \mathrm{d} x ~\Bigg\vert_{t=1}\\ &= \left[\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 3 \frac{\mathrm{d}}{\mathrm{d}t}+1 \right] \frac{1}{t} \int \limits_0^\infty \left[\ln^3 (y) - 3 \ln^2(y) \ln(t) + 3 \ln(y) \ln^2(t) - \ln^3 (t)\right] \mathrm{e}^{-y} \, \mathrm{d} y ~\Bigg\vert_{t=1}\\ &= \left[\frac{\mathrm{d}^2}{\mathrm{d} t^2} + 3 \frac{\mathrm{d}}{\mathrm{d}t}+1 \right] \frac{1}{t} \left[\Gamma'''(1) - 3 \Gamma''(1) \ln(t) + 3 \Gamma'(1) \ln^2(t) - \ln^3 (t)\right] ~\Bigg\vert_{t=1} \\ &= 2 \Gamma'''(1) + 6 \Gamma''(1) + 3 \Gamma''(1) + 6 \Gamma'(1) - 3 \Gamma'''(1) - 9 \Gamma''(1) + \Gamma'''(1) \\ &= 6 \Gamma'(1)\\ &= - 6 \gamma \, . \end{align}
In fact, integration by parts yields the following generalisation: \begin{align} \gamma &= \int \limits_0^\infty (-\ln (x)) \mathrm{e}^{-x} \, \mathrm{d} x = \int \limits_0^\infty \frac{-\ln (x)}{x} x \mathrm{e}^{-x} \, \mathrm{d} x \\ &= \int \limits_0^\infty (-\ln (x))^2 \frac{1-x}{2} \mathrm{e}^{-x} \, \mathrm{d} x \\ &= \int \limits_0^\infty (-\ln (x))^3 \frac{x^2 - 3x +1}{6} \mathrm{e}^{-x} \, \mathrm{d} x \\ &= \dots \, \\ &= \int \limits_0^\infty (-\ln (x))^{n+1} \frac{p_n (x)}{(n+1)!} \mathrm{e}^{-x} \, \mathrm{d} x \, . \end{align} The polynomials $p_n$ are defined recursively by $p_0(x) = 1$ and $$p_n (x) = \mathrm{e}^{x} \frac{\mathrm{d}}{\mathrm{d}x} \left(x p_{n-1} (x) \mathrm{e}^{-x}\right) \, , \, n \in \mathbb{N} \, ,$$ for $x \in \mathbb{R}$ . The exponential generating function $$ \sum \limits_{n=0}^\infty \frac{p_n(x)}{n!} t^n = \mathrm{e}^{t+x(1-\mathrm{e}^t)}$$ can actually be computed from a PDE and it turns out that the polynomials are given by $$p_n(x) = \frac{B_{n+1}(-x)}{-x} \, , \, x \in \mathbb{R} \, , \, n \in \mathbb{N}_0 \, , $$ where $(B_k)_{k \in \mathbb{N}_0}$ are the Bell or Touchard polynomials.