結果

問題 No.463 魔法使いのすごろく🎲
ユーザー rpy3cpprpy3cpp
提出日時 2017-04-18 23:24:03
言語 Python3
(3.11.6 + numpy 1.26.0 + scipy 1.11.3)
結果
AC  
実行時間 55 ms / 2,000 ms
コード長 4,308 bytes
コンパイル時間 90 ms
コンパイル使用メモリ 11,152 KB
実行使用メモリ 9,800 KB
最終ジャッジ日時 2023-09-26 13:28:28
合計ジャッジ時間 2,872 ms
ジャッジサーバーID
(参考情報)
judge15 / judge14
このコードへのチャレンジ(β)

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 21 ms
8,708 KB
testcase_01 AC 20 ms
8,652 KB
testcase_02 AC 20 ms
8,808 KB
testcase_03 AC 19 ms
8,748 KB
testcase_04 AC 22 ms
8,816 KB
testcase_05 AC 28 ms
8,836 KB
testcase_06 AC 19 ms
8,576 KB
testcase_07 AC 20 ms
8,784 KB
testcase_08 AC 20 ms
8,740 KB
testcase_09 AC 48 ms
9,260 KB
testcase_10 AC 20 ms
8,636 KB
testcase_11 AC 26 ms
8,856 KB
testcase_12 AC 22 ms
8,852 KB
testcase_13 AC 18 ms
8,552 KB
testcase_14 AC 46 ms
9,340 KB
testcase_15 AC 37 ms
9,080 KB
testcase_16 AC 19 ms
8,556 KB
testcase_17 AC 41 ms
9,264 KB
testcase_18 AC 19 ms
8,592 KB
testcase_19 AC 19 ms
8,616 KB
testcase_20 AC 19 ms
8,728 KB
testcase_21 AC 55 ms
9,524 KB
testcase_22 AC 55 ms
9,800 KB
testcase_23 AC 54 ms
9,676 KB
testcase_24 AC 52 ms
9,212 KB
testcase_25 AC 53 ms
9,276 KB
testcase_26 AC 53 ms
9,304 KB
testcase_27 AC 53 ms
9,288 KB
testcase_28 AC 55 ms
9,668 KB
testcase_29 AC 53 ms
9,372 KB
testcase_30 AC 55 ms
9,360 KB
testcase_31 AC 54 ms
9,288 KB
testcase_32 AC 55 ms
9,500 KB
testcase_33 AC 53 ms
9,392 KB
testcase_34 AC 53 ms
9,300 KB
testcase_35 AC 19 ms
8,708 KB
testcase_36 AC 19 ms
8,656 KB
testcase_37 AC 20 ms
8,608 KB
testcase_38 AC 19 ms
8,544 KB
権限があれば一括ダウンロードができます

ソースコード

diff #

import copy

def read_data():
    n, m = map(int, input().split())
    if n == 2 and m == 1:
        Cs = [0, 0]
    else:
        Cs = [0] + list(map(int, input().split())) + [0]
    return n, m, Cs

def solve(n, m, Cs):
    '''
    dp1[a]: 位置 a から魔法を使わずにゴールするときのコストの期待値
    dp1[a] = sum(dp1[c] + Cs[c] for b in range(a + 1, a + m + 1))/m
    ただし、 c = b if b < n else 2 * n - 2 - b
    
    A = (1/m) *[[0, 1, 1, 1, ..., 1, 0, 0, 0, 0, 0],
         [0, 0, 1, 1, ..., 1, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, ..., 1, 1, 1, 0, 0, 0],
         ...
         [0, 0, 0, 0, ..., 1, 1, 1, 1, 1, 0],
         [0, 0, 0, 0, ..., 0, 1, 1, 1, 1, 1],
         [0, 0, 0, 0, ..., 0, 0, 1, 1, 2, 1],
         [0, 0, 0, 0, ..., 0, 0, 0, 2, 2, 1],
         [0, 0, 0, 0, ..., 0, 0, 1, 1, 2, 1],
         [0, 0, 0, 0, ..., 0, 1, 1, 1, 1, 1],
         [0, 0, 0, 0, ..., 0, 0, 0, 0, 0, 0]]
    とすると、
    dp1 = A dp1 + A Cs
    <=> (I - A)dp1 = A Cs = (A - I) Cs + Cs = - (I - A)Cs + Cs
    <=> dp1 = -(I - A)^-1 (I - A) Cs + (I - A)^-1 Cs
            = - Cs + (I - A)^-1 Cs
    <=> dp1 + Cs = (I - A)^-1 Cs
    
    dp2 = dp1 + Cs とする。
    
    dp0[a]: 位置 a からゴールするまでのコストの期待値の最小値(魔法を1度使ってもよい)
    dp0[a] = min(min(dp1[b] + Cs[b] for b in range(a + 1, min(a + m + 1, n))),
                 sum(dp0[c] + Cs[c] for b in range(a + 1, a + m + 1))/m)
    ただし、 c = b if b < n else 2 * n - 2 - b
    
    dp0[n - 1] = 0      # ゴールしている。
    dp0[n - 2] = 0      # 魔法で 1 を出せば、ゴールできる。
    dp0[n - 3] = 0      # 魔法で 2 を出せば、ゴールできる。
    ...
    dp0[n - m - 1] = 0  # 魔法で m を出せば、ゴールできる。
    dp0[n - m - 2] = min(dp2[b] for b in range(n - m - 2 + 1, n - 1)),
                         sum(dp0[b] + Cs[b] for b in range(n - m - 2 + 1, n - m - 2 + m + 1))/m)
    dp0[a] = min(dp2[b] for b in range(a + 1, a + m + 1)),
                 sum(dp0[b] + Cs[b] for b in range(a + 1, a + m + 1))/m)
    '''
    A = calc_matrix_A(n, m)
    I = identity_matrix(n)
    ImA = mat_minus(I, A)
    colCs = [[c] for c in Cs]
    det, dp2 = gauss(ImA, colCs)
    dp2 = [a[0] for a in dp2]
    dp0 = [0] * n
    for a in range(n - m - 2, -1, -1):
        min_dp2 = min(dp2[a + 1:a + m + 1])
        avg_dp0 = sum(dp0[b] + Cs[b] for b in range(a + 1, a + m + 1))/m
        dp0[a] = min(min_dp2, avg_dp0)
    return dp0[0]
    
def calc_matrix_A(n, m):
    A = [[0] * n for _ in range(n)]
    for r, row in enumerate(A):
        if r == n - 1:
            break
        for c in range(r + 1, r + 1 + m):
            row[b2c(c, n)] += 1/m
    return A

def identity_matrix(n):
    I = [[0] * n for _ in range(n)]
    for i in range(n):
        I[i][i] = 1
    return I

def mat_minus(A, B):
    C = [[a - b for a, b in zip(rowA, rowB)] for rowA, rowB in zip(A, B)]
    return C

def b2c(b, n):
    if b < n:
        return b
    else:
        return 2 * n - 2 - b

def gauss(a, b):
    '''The 'gauss' function takes two matrices, 'a' and 'b', with 'a' square, and it return the determinant of 'a' and a matrix 'x' such that a*x = b.
    If 'b' is the identity, then 'x' is the inverse of 'a'.
    Source: Rosetta Code
    http://rosettacode.org/wiki/Gaussian_elimination#Python
    '''
    a = copy.deepcopy(a)
    b = copy.deepcopy(b)
    n = len(a)
    p = len(b[0])
    det = 1
    for i in range(n - 1):
        k = i
        for j in range(i + 1, n):
            if abs(a[j][i]) > abs(a[k][i]):
                k = j
        if k != i:
            a[i], a[k] = a[k], a[i]
            b[i], b[k] = b[k], b[i]
            det = -det

        for j in range(i + 1, n):
            t = a[j][i]/a[i][i]
            for k in range(i + 1, n):
                a[j][k] -= t*a[i][k]
            for k in range(p):
                b[j][k] -= t*b[i][k]

    for i in range(n - 1, -1, -1):
        for j in range(i + 1, n):
            t = a[i][j]
            for k in range(p):
                b[i][k] -= t*b[j][k]
        t = 1/a[i][i]
        det *= a[i][i]
        for j in range(p):
            b[i][j] *= t
    return det, b


if __name__ == '__main__':
    n, m, Cs = read_data()
    print(solve(n, m, Cs))
0