30分求条
查看原帖
30分求条
615573
HopeIsntHere楼主2025/1/12 20:04
// 7, 8, 18, 19, 20 MLE
// 9 ~ 17 TLE
#include <cstdio>
#include <vector>
#include <utility>
#include <cmath>

#define pb push_back

using std::vector;
using std::swap;

struct {
    operator int() {
        int x;
        scanf("%d", &x);
        return x;
    }
} it;

template <typename T>
struct mat {
    int n, m;
    T v;
    vector<T> w;

    mat(int n, int m, T v = T {}) : n(n), m(m), v(v), w(n * m) {}

    T& operator()(int i, int j) {
        if (i < 1 || i > n || j < 1 || j > m)
            return v;

        return w[(i - 1) * m + (j - 1)];
    }

    const T& operator()(int i, int j) const {
        if (i < 1 || i > n || j < 1 || j > m)
            return v;

        return w[(i - 1) * m + (j - 1)];
    }

    mat<T> t() const {
        mat c(m, n);
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < m; ++j)
                c.w[j * n + i] = w[i * m + j];
        return c;
    }
};

struct vec {
    vector<double> w;

    vec() = default;

    explicit vec(int n) : w(n) {}

    int size() const {
        return w.size();
    }

    double& operator[](int i) {
        return w[i];
    }

    double operator[](int i) const {
        return w[i];
    }

    double& at(int i) {
        return w[i < 0 ? i + w.size() : i];
    }

    const double& at(int i) const {
        return w[i < 0 ? i + w.size() : i];
    }

    vec& operator+=(const vec& r) {
        for (size_t i = 0; i != w.size(); ++i)
            w[i] += r.w[i];

        return *this;
    }

    vec& operator-=(const vec& r) {
        for (size_t i = 0; i != w.size(); ++i)
            w[i] -= r.w[i];

        return *this;
    }

    vec& operator*=(double k) {
        for (size_t i = 0; i != w.size(); ++i)
            w[i] *= k;

        return *this;
    }
};

vec operator+(const vec& a, const vec& b) {
    vec c = a;
    c += b;
    return c;
}

vec operator-(const vec& a, const vec& b) {
    vec c = a;
    c -= b;
    return c;
}

vec operator*(double k, const vec& a) {
    vec c = a;
    c *= k;
    return c;
}

double dot(const vec& a, const vec& b) {
    double s = 0;

    for (int i = 0; i < a.size(); ++i)
        s += a[i] * b[i];

    return s;
}

vec sol(const vector<vec>& a) {
    auto f = a;
    int n = f.size();

    for (int i = 0; i < n; ++i) {
        if (fabs(f[i][i] < 1e-4)) {
            int k = -1;

            for (int j = i + 1; j < n; ++j)
                if (fabs(f[j][i]) > 1e-4)
                    k = j;

            swap(f[i], f[k]);
        }

        f[i] *= 1 / f[i][i];

        for (int j = 0; j < n; ++j)
            if (i != j)
                f[j] -= f[j][i] * f[i];

    }

    int m = f[0].size();
    vec x(m);

    for (int i = 0; i < n; ++i)
        x[i] = -f[i].at(-1);

    x.at(-1) = 1;
    return x;
}

template <typename F1, typename F2>
mat<double> exp(int n, int m, F1 d, F2 b) {
    auto base = [n, m] (int i) {
        vec c(n * ((m - 1) / 2 + 1) + 2);
        c.at(i) = 1;
        return c;
    };

    mat<vec> f(n, m);
    int k = 0;

    for (int i = 1; i <= n; ++i)
        f(i, 1) = base(k++);

    vector<vec> r;
    for (int j = 1; j <= m; ++j) {
        for (int i = 1; i <= n; ++i) {
            /*
            f[i,j] = f[i,j-1]/d[i,j-1] + f[i-1,j]/d[i-1,j] + f[i+1,j]/d[i+1,j] + f[i,j+1]/d[i,j+1] + b[i,j]
            f[i,j] - f[i,j-1]/d[i,j-1] - f[i-1,j]/d[i-1,j] - f[i+1,j]/d[i+1,j] - b[i,j] = f[i,j+1]/d[i,j+1]
            */
            vec s = f(i, j);

            if (b(i, j))
                s.at(-1) -= b(i, j);

            if (d(i-1, j))
                s -= (1.0 / d(i-1, j)) * f(i-1, j);

            if (d(i+1, j))
                s -= (1.0 / d(i+1, j)) * f(i+1, j);

            if (d(i, j-1))
                s -= (1.0 / d(i, j-1)) * f(i, j-1);

            if (j == m)
                r.pb(s);

            else if (!d(i, j+1)) {
                f(i, j+1) = base(k++);
                r.pb(s);
            }

            else if ((j + 1) % 4 < 2) {
                f(i, j+1) = base(k++);
                s -= (1.0 / d(i, j+1)) * f(i, j+1);
                r.pb(s);
            }
            else
                f(i, j+1) = d(i, j+1) * s;
        }
    }
    vec x = sol(r);
    mat<double> e(n, m);
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= m; ++j) {
            e(i, j) = dot(f(i, j), x);
        }
    return e;
}

int main() {
    int n = it, m = it;
    mat<int> v(n - 1, m);
    mat<int> h(n, m - 1);
    for (int i = 1; i <= n - 1; ++i)
        for (int j = 1; j <= m; ++j) {
            v(i, j) = it;
        }
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= m - 1; ++j) {
            h(i, j) = it;
        }
    bool t = n > m;
    if (t) {
        swap(n, m);
        v = v.t();
        h = h.t();
    }
    for (int l = it; l >= 1; --l) {
        int x1 = it + 1, y1 = it + 1, x2 = it + 1, y2 = it + 1;
        if (t) {
            swap(x1, y1);
            swap(x2, y2);
        }
        auto d = [n, m, x2, y2] (int i, int j) -> int {
            if (i < 1 || i > n || j < 1 || j > m || (i == x2 && j == y2)) {
                return 0;
            }
            return 4 - (i == 1) - (i == n) - (j == 1) - (j == m);
        };
        auto b = [x1, y1] (int i, int j) -> int {
            return i == x1 && j == y1;
        };
        mat<double> e = exp(n, m, d, b);
        auto c = [&e, d] (int i, int j) -> double {
            return d(i, j) ? e(i, j) / d(i, j) : 0;
        };
        double s = 0;
        for (int i = 1; i <= n - 1; ++i)
            for (int j = 1; j <= m; ++j) {
                s += v(i, j) * (c(i, j) + c(i+1, j));
            }
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= m - 1; ++j) {
                s += h(i, j) * (c(i, j) + c(i, j+1));
            }
        printf("%.4f\n", s);
    }
}
2025/1/12 20:04
加载中...