LogzetaN(s, N) =
{ my(B = getlocalbitprec());

B += ceil(abs(real(s) - 1) * exponent(2 * N));
localbitprec(B);
if (bitprecision(s) < B, s = bitprecision(s, B));
log(zeta(s) * prodeuler(p = 2, N, 1 - p^(-s)));
}

/* Compute sum_{p prime} log(log(p))/(p^s log(p)). */
SumEulerloglog(s, N = max(2, 30 / abs(s))) =
{ my(B = getlocalbitprec(), S = 0, LN, T, lim, E);

localbitprec(32); LN = log(N);
lim = ceil(B*log(2) / LN);
localbitprec(B + 32);
forprime (p = 2, N, S += log(log(p)) / (p^s * log(p)));
LN = bitprecision(LN, B + 32);
T = intnuminit(0,[oo, 1]);
E = Euler();
forsquarefree (K = 1, lim,
my([k] = K, m = moebius(K) / k, a = 1 / (k * LN));
/* Tk = intnuminit(0, [oo, 1/a]) */
my(Tk = vector(#T, i, if (i == 1, T[1], T[i] * a)));
S -= m * intnum(t = 0, oo, (log(t) + E)
* LogzetaN(k * (t+s) , N), Tk));
return (S);
}

SumEulerloglog(1)