(Код внизу) Кому интересно, вот сама статья: http://numbers.computation.free.fr/Constants/Algorithms/nthdecimaldigit.pdf

Imgur

Real compute(Int n) {
        Real n0 = 10;
        Int M = 2*ceil(n/pow(log(n), 3));
        Int N = ceil((n + n0 + 1)*(log(10)/log(2*exp(1)*M)));

        Real b = 0;
        for (Int k = 0; k < (M + 1)*N; ++k) {
            Real x = 4*modpow<Int>(10, n, 2*k + 1);
            b = fractionalpart(b + pow(-1, k)*x/(2*k + 1));
        }

        Real c = 0, x = 0, y = 0;
        for (Int j = 0; j < N; ++j) {
            Int modulus = 2*M*N + 2*j + 1;
            x += modbinominals(N, j, modulus);
            y = Int(pow(5, N - 2) * pow(10, n - N + 2) * x) % modulus;
            c = fractionalpart(c + pow(-1, j)*y/modulus);
        }

        return fractionalpart(b - c) ;
    }

Остальные функции:

using Int = long long;
using Real = long double;

inline Real fractionalpart(Real x) {
    return x - floor(x);
}

Int modbinominals(Int n, Int j, Int modulus) {
    Int result = 1;
    for (; j > 0; --j)
        result = (result * (n - j + 1)/j) % modulus;
    return result;
}

template <typename T>
T modpow(T base, T exp, T modulus) {
    base %= modulus;
    T result = 1;    
    while (exp > 0) {
        if (exp & 1) result = (result * base) % modulus;
        base = (base * base) % modulus;
        exp >>= 1;
    }
    return result;
}

Но...

cout << compute(4) << endl;

0.430662

(3,14159265358979323846)

Edit Report
Pub: 19 Feb 2020 19:22 UTC
Edit: 19 Feb 2020 19:38 UTC
Views: 280