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} 1/(p^s log(p)). */
SumEulerlog(s, N = max(2, 30/abs(s))) =
{ my(B = getlocalbitprec(), S = 0, LN, T, lim);

  localbitprec(32); LN = log(N);
  lim = ceil(B*log(2)/LN);
  localbitprec(B + 32);
  forprime(p = 2, N, S += 1 / (p^s * log(p)));
  LN = bitprecision(LN, B + 32);
  T = intnuminit(0,[oo, 1]);
  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, LogzetaN(k * (t + s) , N), Tk));
  return (S);
}