#include #include #include #include #include #include typedef __float128 real; template struct Adder { std::vector 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 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; }