Pasted

#include <quadmath.h>
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <vector>
  typedef __float128 real;
  template <typename T> struct Adder {
	std::vector<T> v;
 	void reset() {
		v.clear();
		v.push_back(0);
	}
 	Adder() { reset(); }
 	T sum() const {
		T res = v[0];
		size_t l = v.size();
		for (size_t i = 1; i < l; i++)
			res += v[i];
		return res;
	}
 	void operator+= (T x) {
		size_t i = 0;
		for (i = 0; v[i]; i++) {
			x += v[i];
			v[i] = 0;
		}
		v[i] = x;
		if (i + 1 == v.size())
			v.push_back(0);
	}
 };
  int main()
{
	uint64_t n = 100000000;
	Adder<real> a;
	for (uint64_t i = 2; i <= n; i++) {
		real x = i;
		real y = expm1q(logq(x) / x) / x;
		a += y;
	}
	real s = a.sum();
 	{
		real x = n + 0.5q;
		real lx = logq(x);
		real y1 = (lx + 1) / x;
		real y2 = (2 * lx * (lx + 1) + 1) / (2 * 4 * x * x);
		real y3 = (((9 * lx + 9) * lx + 6) * lx + 2) / (6 * 27 * x * x * x);
		real y4 = ((((32 * lx + 32) * lx + 24) * lx + 12) * lx + 3) / (24 * 128 * (x * x) * (x * x));
		real y5 = (((((625 * lx + 625) * lx + 500) * lx + 300) * lx + 120) * lx + 24) / (120 * 3125 * x * (x * x) * (x * x));
		real xrx = expm1q(lx / x);
		real yd = (x * xrx + (xrx + 1) * (lx - 1)) / (24 * x * x * x);
		s += y1;
		s += y2;
		s += y3;
		s += y4;
		s += y5;
		s -= yd;
	}
 	printf("n=%lu, %.36Qf\n", n, s);
	return 0;
}