walk.cpp

#include <vector>
#include <iostream>
#include <iomanip>

struct Row {
    int off;
    std::vector<double> buf{};
    double b = -4.;
};

void show(std::vector<Row> const& row, std::vector<int> const& maxy, int T) {
    for (auto y = 0; y != T; ++y) {
        auto const& r = row[y];
        for (auto x = 0; x != T; ++x) {
            if (x < row[y].off || x - row[y].off >= row[y].buf.size())
                std::cout << "  ";
            else if (r.buf[x - r.off] == 0)
                std::cout << ". ";
            else
                std::cout << r.buf[x - r.off] << ' ';
        }
        std::cout << "| " << r.b << '\n';
    }
    for (auto x = 0; x != T; ++x)
        std::cout << maxy[x] << ' ';
    std::cout << '\n';
};

double factor(int n) {
    auto const N = n;
    auto const T = N * (N + 1) / 2;
    auto row = std::vector<Row>(T);
    auto maxy = std::vector<int>(T);
    auto A = [&](int x, int y) -> double& { auto& r = row[y]; return r.buf[x - r.off]; };
    auto B = [&](int y) -> double& { return row[y].b; };
    auto p = [=](int i, int j) { return i + (j * (j + 1) / 2); };
    auto maxx = 0;
    for (auto j = 0; j != N; ++j)
        for (auto i = 0; i != j + 1; ++i) {
            auto const y = p(i, j);
            auto& r = row[y];
            for (auto l = 0; l != N; ++l)
                for (auto k = 0; k != l + 1; ++k) {
                    auto const x = p(k, l);
                    auto const c = (k == i && l == j) ? -4 :
                        (i == N - 1 && j == N - 1 && k == N - 2 && l == N - 1) ? 4 :
                        (((k == i - 1 || k == i + 1) && (l == j)) || (k == i && (l == j - 1 || l == j + 1))) +
                        ((j == N - 1 && l == N - 2 && k == i)) +
                        (i == j && ((k == i && l == j + 1) || (k == i - 1 && l == j)));
                    if (c != 0) {
                        if (r.buf.empty())
                            r.off = x;
                        maxx = std::max(maxx, x);
                        r.buf.resize(maxx - r.off + 1);
                        r.buf[x - r.off] = c;
                        maxy[x] = std::max(maxy[x], y);
                    }
                    else if (!r.buf.empty())
                        maxy[x] = std::max(maxy[x], y);
                }
        }

    for (auto x = 0; x != T; ++x)
        for (auto y = x + 1; y != maxy[x] + 1; ++y)
            if (row[y].off <= x && A(x, y) != 0) {
                auto const k = A(x, y) / A(x, x);
                auto const maxx = row[x].off + row[x].buf.size();
                for (auto i = x; i != maxx; ++i)
                    A(i, y) -= A(i, x) * k;
                B(y) -= B(x) * k;
            }
    return B(T - 1) / A(T - 1, T - 1) / (N*N);
}

int main()
{
    for (int i = 2; i != 10; ++i)
        std::cout << i << " " << std::setprecision(17) << factor(i) << std::endl;
    for (int i = 10; i != 100; i += 10)
        std::cout << i << " " << std::setprecision(17) << factor(i) << std::endl;
    for (int i = 100; i != 1000; i += 100)
        std::cout << i << " " << std::setprecision(17) << factor(i) << std::endl;
    for (int i = 1000; i != 10000; i += 1000)
        std::cout << i << " " << std::setprecision(17) << factor(i) << std::endl;
    for (int i = 10000; i != 100000; i += 10000)
        std::cout << i << " " << std::setprecision(17) << factor(i) << std::endl;
}