Ytukyg

We denote for integers $ngeq 1$ the $n$th Gregory coefficient as $G_n$, and the Möbius function as $mu(n)$. You've here the Wikipedia's article dedicated to the Gregory coefficients.

Using an argument of absolute convergence, and the information of previous Wikipedia for the first related series to the Gregory coefficients and the result due to Candelperger, Coppo and Young, it is obvious to prove that $$sum_{n=1}^inftymu(n)|G_n|tag{1}$$ and $$sum_{n=1}^inftyfrac{|G_n|cdot m(n)}{n}tag{2}$$ are convergent series, where $m(x)$ denotes the function $$m(x)=sum_{1leq kleq x}frac{mu(k)}{k}.tag{3}$$

Question. Have you an idea/hint to get a good approximation (the first four or six digits) of $(1)$ and $(2)$? Many thanks.

I know that there are upper bounds for the absolute value of $(3)$ for large values of $x$.

Actually much more than "four or six digits" is reachable. I'm showing this for $(1)$ below.

Let $Xi=displaystylesum_{n=1}^{infty}mu(n)|G_n|$. Using $|G_n|=displaystyleint_{0}^{infty}frac{dx}{(1+x)^n(pi^2+ln^2x)}$, we get
$$Xi = int_{0}^{infty}FBig(frac{1}{1+x}Big)frac{dx}{pi^2+ln^2x},quadcolor{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 $xinmathbb{R}_{>0}$ and $cinmathbb{R}_{>1}$ we have
$$F(e^{-x})=frac{1}{2pi i}int_{c-iinfty}^{c+iinfty}frac{Gamma(s)}{x^szeta(s)},ds$$
which, after somewhat boring computations of residues, arrives at
$$begin{gather}color{blue}{F(e^{-x})}=-2+sum_{omegainOmega}operatorname*{Res}_{s=omega}frac{Gamma(s)}{x^szeta(s)}+sum_{n=1}^{infty}frac{2nx^{2n-1}}{(2n-1)!B_{2n}}\ {}+2sum_{n=1}^{infty}left[frac{(-1)^n(2pi x)^{2n}}{(2n)!^2zeta(2n+1)}left(frac{zeta'(2n+1)}{zeta(2n+1)}-ln 2pi x + 2(H_{2n}-gamma)right)right]end{gather}$$
with $Omega={omegainmathbb{C}setminusmathbb{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 $omegainOmega$ 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.3600138625016611865745170005656289245070028602995555633ldots$$