結果

問題 No.463 魔法使いのすごろく🎲
ユーザー izuru_matsuuraizuru_matsuura
提出日時 2016-12-15 02:25:16
言語 D
(dmd 2.109.1)
結果
AC  
実行時間 9 ms / 2,000 ms
コード長 4,604 bytes
コンパイル時間 3,563 ms
コンパイル使用メモリ 171,904 KB
実行使用メモリ 5,376 KB
最終ジャッジ日時 2024-06-12 05:33:50
合計ジャッジ時間 4,976 ms
ジャッジサーバーID
（参考情報） 		judge2 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 36
権限があれば一括ダウンロードができます

ソースコード

diff # 

import std.algorithm;
import std.array;
import std.ascii;
import std.container;
import std.conv;
import std.math;
import std.numeric;
import std.range;
import std.stdio;
import std.string;
import std.typecons;
void log(A...)(A arg) { stderr.writeln(arg); }
int size(T)(in T s) { return cast(int)s.length; }

real dot(T)(in T[] a, in T[] b) {
    assert(a.length == b.length);
    return reduce!((sum, i) => sum + a[i] * b[i])(cast(T)0, iota(0, a.length, 1));
}
class Matrix(T) {
    int H, W;
    T[] data;
    this(int H, int W) {
        this.H = H;
        this.W = W;
        this.data = new T[H * W];
    }
    this(int H, int W, T init) {
        this(H, W);
        foreach (i; 0 .. H) foreach (j; 0 .. W) at(i, j) = init;
    }
    this(int H, int W, T delegate(int, int) init) {
        this(H, W);
        foreach (i; 0 .. H) foreach (j; 0 .. W) at(i, j) = init(i, j);
    }
    //pragma(inline, true) {
        ref T at(int i, int j) { return data[i * W + j]; }
        T get(int i, int j) const { return data[i * W + j]; }
        const(T[]) get(int i) const { return data[i * W .. (i + 1) * W]; }
    //}
    Matrix dup() const { return new Matrix(H, W, (i, j) => get(i, j)); }
    Matrix swapRows(int i, int j) { foreach (k; 0 .. W) swap(at(i, k), at(j, k)); return this; }
    Matrix swapCols(int i, int j) { foreach (k; 0 .. H) swap(at(k, i), at(k, j)); return this; }
    Matrix transpose() const { return new Matrix(H, W, (i, j) => get(j, i)); }
    Matrix opUnary(string op)() const if (op == "-") {
        return new Matrix(H, W, (i, j) => -get(i, j));
    }
    Matrix opBinary(string op)(in Matrix B) const if (op == "+" || op == "-") {
        assert(H == B.H && W == B.W);
        return new Matrix(H, W, (i, j) => mixin("get(i, j) " ~ op ~ " B.get(i, j)"));
    }
    Matrix opBinary(string op)(in Matrix B) const if (op == "*") {
        assert(W == B.H);
        auto X = B.transpose;
        auto A = new Matrix(H, B.W, 0);
        foreach (i; 0 .. H) foreach (j; 0 .. X.H) {
            foreach (k; 0 .. W) {
                A.at(i, j) += get(i, k) * X.get(j, k);
            }
        }
        return A;
    }
    T[] opBinary(string op)(in T[] b) const if (op == "*") {
        auto y = new T[b.length];
        foreach (i; 0 .. H) y[i] = dot(get(i), b);
        return y;
    }
    override string toString() const { return to!string(iota(0, H, 1).map!((i) => data[i * W .. (i + 1) * W]).array); }
}
immutable EPS = 1e-4;
bool EQ(real a, real b) { return abs(a - b) < EPS; }
real[] gaussianElimination(in Matrix!real X, in real[] b) {
    // 方程式Xy = bの解yを求める. 
    assert(X.H == b.size);
    assert(X.H == X.W);
    auto N = X.H;
    auto A = X.dup;
    auto c = b.dup;
    foreach (int i; 0 .. N) {
        int p = i;
        foreach (int j; i .. N) {
            if (abs(A.get(i, j)) > abs(A.get(p, j))) {
                p = j;
            }
        }
        A.swapRows(p, i);
        real x = A.get(i, i);
        if (abs(x) < EPS) continue;
        foreach (j; i .. N) A.at(i, j) /= x;
        c[i] /= x;
        foreach (k; i + 1 .. N) {
            real y = A.get(k, i);
            c[k] -= y * c[i];
            foreach (l; i .. N) {
                A.at(k, l) -= y * A.get(i, l);
            }
        }
    }
    auto x = new real[N];
    foreach_reverse (j; 0 .. N) {
        auto k = A.get(j, j);
        x[j] = c[j] / k;
        foreach (i; 0 .. j) {
            c[i] -= A.get(i, j) * x[j];
        }
    }
    return x;
}

void main() {
    int N, M; readf("%s %s\n", &N, &M);
    auto C = [cast(real)0] ~ readln.chomp.split(" ").map!(to!real).array;
    auto A = new Matrix!real(N - 1, N - 1, 0);
    auto forward = (int j, int c) => (j + c < N ? j + c : (N - 1) - (j + c - (N - 1)));
    foreach (i; 0 .. N - 1) {
        for (int c = 1, j = forward(i, c); c <= M; c++, j = forward(i, c)) {
            if (j == N - 1) continue;
            A.at(i, j) += 1.0 / M;
        }
    }
    auto X = A.dup;
    foreach (i; 0 .. N - 1) X.at(i, i) -= 1.0;
    //log("A: ", A);
    //log("X: ", X);
    //log("-A*C: ", -A * C);
    auto E = gaussianElimination(X, -A * C);
    //log(E);

    E ~= cast(real)0;
    C ~= cast(real)0;
    auto F = E.dup;
    foreach_reverse (i; 0 .. N - 1) {
        for (int k = i + 1; k <= min(i + M, N - 1); k++) {
            F[i] = min(F[i], E[k] + C[k]);
            real s = 0;
            for (int c = 1, j = forward(i, c); c <= M; c++, j = forward(i, c)) {
                s += F[j] + C[j];
            }
            s /= M;
            F[i] = min(F[i], s);
        }
    }
    writefln("%.12f", F[0]);
}
0