TopCoder SRM 614: TorusSailing

解法は TopCoder SRM614 解説 を参考にしてください。

流れは、だいたい以下のとおり。

1. 全てのループを断ち切るように N + M - 1 箇所のマスを変数化
   • 斜めに変数化すれば max(N, M) 箇所でよい
2. すべての箇所の期待値を、DP を使って (N + M - 1) 変数多項式で表す
3. 上記の結果を元に (N + M - 1) 変数連立方程式に帰着できるため、変数の数が減ることによりガウスジョルダンや掃き出し法で解ける。

よく期待値問題で、F(n) = a * F(n - 1) + b * F(n)のような自己ループを含む漸化式になった場合、F(n) = a * F(n - 1) / (1 - b)のように移項して解くが、このとき F(n) をいったん変数化して一次方程式を解いていることになる。本問題は、このテクを一次方程式から連立方程式に拡張した問題として位置づけられる。

public class TorusSailing {
    static public void Swap<T>(ref T a, ref T b) { T t = a; a = b; b = t; }

    class Polynomial {
        public double[] Coefficient;
        public double constantTerm;

        protected Polynomial() { }

        public Polynomial(int variableNum) {
            this.Coefficient = new double[variableNum];
            this.constantTerm = 0.0;
        }

        public void Add(Polynomial p) {
            for (int i = 0; i < Coefficient.Length; i++) {
                Coefficient[i] += p.Coefficient[i];
            }
            constantTerm += p.constantTerm;
        }

        public void Divide(double p) {
            for (int i = 0; i < Coefficient.Length; i++) {
                Coefficient[i] /= p;
            }
            constantTerm /= p;
        }
    }

    Polynomial GetBasePolynomial(int posX, int posY, int N, int M) {
        if (posY == 0) {
            var p = new Polynomial(N + M - 1);
            p.Coefficient[posX] = 1.0;
            return p;
        } else if (posX == 0) {
            var p = new Polynomial(N + M - 1);
            p.Coefficient[posY - 1 + N] = 1.0;
            return p;
        }
        return null;
    }

    Polynomial[,] memo = new Polynomial[100, 100];

    Polynomial ExpectedPolynomial(int posX, int posY, int goalX, int goalY, int N, int M) {
        if (memo[posX, posY] == null) {
            var p = new Polynomial(N + M - 1);
            if (posX == goalX && posY == goalY) {
                ; // p = 0
            } else {
                int nX = (posX + 1) % N;
                int nY = (posY + 1) % M;

                var p1 = GetBasePolynomial(nX, posY, N, M) ?? ExpectedPolynomial(nX, posY, goalX, goalY, N, M);
                var p2 = GetBasePolynomial(posX, nY, N, M) ?? ExpectedPolynomial(posX, nY, goalX, goalY, N, M);

                p.Add(p1); p.constantTerm += 1.0;
                p.Add(p2); p.constantTerm += 1.0;

            }
            p.Divide(2.0);
            memo[posX, posY] = p;
        }
        return memo[posX, posY];
    }

    double[] GausianJordanElimination(double[,] mat) {
        int N = mat.GetLength(1) - 1;
        int M = mat.GetLength(0);
        var index = new int[M];
        const double Error = 1e-9;
        for (int i = 0; i < M; i++) {
            index[i] = i;
        }
        for (int i = 0; i < N; i++) {
            int row = -1;
            double pivot = Error;
            for (int j = i; j < M; j++) {
                if (Math.Abs(mat[index[j], i]) > Math.Abs(pivot)) {
                    row = j;
                    pivot = mat[index[j], i];
                }
            }
            if (row == -1) continue;
            Swap(ref index[i], ref index[row]);
            for (int j = 0; j <= N; j++) {
                mat[index[i], j] /= pivot;
            }
            for (int j = 0; j < M; j++) {
                if (i != j) {
                    double weight = mat[index[j], i];
                    for (int k = 0; k <= N; k++) {
                        mat[index[j], k] -= weight * mat[index[i], k];
                    }
                }
            }
        }
        double[] ans = new double[N];
        for (int i = 0; i < N; i++) {
            ans[i] = mat[index[i], N];
        }
        return ans;
    }

    public double expectedTime(int N, int M, int goalX, int goalY) {
        int size = N + M - 1;
        double[,] matrix = new double[size, size + 1];
        for (int i = 0; i < size; i++) {
            int x, y;
            if (i < N) {
                x = i;
                y = 0;
            } else {
                x = 0;
                y = i - N + 1;
            }
            var p = ExpectedPolynomial(x, y, goalX, goalY, N, M);
            for (int j = 0; j < size; j++) {
                matrix[i, j] = p.Coefficient[j];
                if (i == j) matrix[i, j] -= 1.0;
            }
            matrix[i, size] = -p.constantTerm;
        }

        var ans = GausianJordanElimination(matrix);
        return ans[0];
    }
}