Actually much more than "four or six digits" is reachable. I'm showing this for $(1)$ below.
Let $\Xi=\displaystyle\sum_{n=1}^{\infty}\mu(n)|G_n|$. Using $|G_n|=\displaystyle\int_{0}^{\infty}\frac{dx}{(1+x)^n(\pi^2+\ln^2x)}$, we get
$$\Xi = \int_{0}^{\infty}F\Big(\frac{1}{1+x}\Big)\frac{dx}{\pi^2+\ln^2x},\quad\color{blue}{F(z)=\sum_{n=1}^{\infty}\mu(n)z^n};$$
to apply numeric integration (say, with double-exponential method), we must have a fast enough computation of $F(z)$, especially for $z$ close to $1$. For any $x\in\mathbb{R}_{>0}$ and $c\in\mathbb{R}_{>1}$ we have
$$F(e^{-x})=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}\frac{\Gamma(s)}{x^s\zeta(s)}\,ds$$
which, after somewhat boring computations of residues, arrives at
$$\begin{gather}\color{blue}{F(e^{-x})}=-2+\sum_{\omega\in\Omega}\operatorname*{Res}_{s=\omega}\frac{\Gamma(s)}{x^s\zeta(s)}+\sum_{n=1}^{\infty}\frac{2nx^{2n-1}}{(2n-1)!B_{2n}}\\ {}+2\sum_{n=1}^{\infty}\left[\frac{(-1)^n(2\pi x)^{2n}}{(2n)!^2\zeta(2n+1)}\left(\frac{\zeta'(2n+1)}{\zeta(2n+1)}-\ln 2\pi x + 2(H_{2n}-\gamma)\right)\right]\end{gather}$$
with $\Omega=\{\omega\in\mathbb{C}\setminus\mathbb{R}:\zeta(\omega)=0\}$ the set of "nontrivial zeros of Riemann zeta", and other known species ($B_{2n}$ are Bernoulli numbers, $H_{2n}$ are harmonic numbers, and $\gamma$ is Euler's constant). With $x$ small, this converges much faster than the original series (despite looking that complicated — and I've checked it numerically before going further).
This allows to compute $\Xi$ using, e.g., PARI/GP. I've started with
$$\Xi=\frac{1}{2}-\int_{0}^{\infty}\frac{1-e^x F(e^{-x})}{\pi^2+\ln^2(e^x-1)}dx,$$
split $\int_{0}^{\infty}=\int_{0}^{1}+\int_{1}^{\infty}$, computed the second (via PARI's $\texttt{intnum}$) using the definition of $F$, and the first using the cumbersome formula above (each integral in the sum over $\omega\in\Omega$ needs to be computed separately, due to oscillating behaviour of the integrand, and moreover, yet another substitution $x=e^{-t}$ is needed for this to keep the accuracy).
mgsDenom(z) = Pi^2 + log(z)^2;
mgsDoubleExpo(x) =
{
my (z = exp(-x));
return (exp(z - x / 2) / mgsDenom(exp(z) - 1))
};
mgsNewtonRoot(f, z) =
{
my (e, r = f(z));
until (e <= norm(r), e = norm(r); z -= r / f'(z); r = f(z));
return (z)
};
mgsRegularPartInit(rbp) =
{
my (N = 0, c = 2.0 ^ (rbp + 10)); while (1, N += 1;
c *= (Pi / N / (N + N - 1)) ^ 2; if (c < 1, break));
my (ctx = matrix(N, 3)); c = 2 * Euler + log(2 * Pi);
for (n = 1, N,
my (m = n + n, h = sum(k = 1, m, 2.0 / k));
my (zv = zeta(m + 1), zp = zeta'(m + 1));
ctx[n, 1] = n / factorial(m - 1) / bernreal(m);
ctx[n, 2] = (-1)^n * (2 * Pi)^m / factorial(m)^2 / zv;
ctx[n, 3] = h - c + zp / zv);
return (ctx)
};
mgsRegularPart(x, ctx) =
{
my (ex = exp(x), lx = log(x), rs = 0);
forstep (n = matsize(ctx)[1], 1, -1, rs = x * (ctx[n, 1]
+ x * (rs + ctx[n, 2] * (ctx[n, 3] - lx))));
return ((1 + 2 * ex * (1 - rs)) / mgsDenom(ex - 1))
};
MoebiusGregorySum() =
{
my (rbp = default(realbitprecision));
my (zzz = exp(1) - 1, eps = 0.5 ^ rbp);
my (result = 0.5 + sum(n = 2, rbp, moebius(n) * intnum(
z = zzz, [+oo, -n], (1 + z)^(-n) / mgsDenom(z))));
my (ctx = mgsRegularPartInit(rbp));
result -= intnum(x = 0, 1, mgsRegularPart(x, ctx));
my (a = 0.5 + 14.0 * I, h = 0.1 * I);
my (Pv, Cv = +oo, Nv = norm(zeta(a)));
while (1, Pv = Cv; Cv = Nv; Nv = norm(zeta(a + h));
if (Cv < Pv && Cv < Nv,
my (z = mgsNewtonRoot(zeta, a), t = imag(z), c = gamma(z) / zeta'(z));
my (rv = real(c) * intnum(x = 0, [+oo, +t * I], mgsDoubleExpo(x) * cos(t * x)));
my (iv = imag(c) * intnum(x = 0, [+oo, -t * I], mgsDoubleExpo(x) * sin(t * x)));
my (d = 2 * (rv + iv)); result += d; if (d < eps, break));
a += h);
return (result)
};
MoebiusGregorySum()
This way I get
$$\color{blue}{\sum_{n=1}^{\infty}\mu(n)|G_n|}=0.3600138625016611865745170005656289245070028602995555633\ldots$$
