結果

問題 No.463 魔法使いのすごろく🎲
コンテスト
ユーザー rpy3cpp
提出日時 2017-04-18 23:19:14
言語 Python3
(3.13.1 + numpy 2.2.1 + scipy 1.14.1)
結果
RE  
実行時間 -
コード長 4,248 bytes
コンパイル時間 108 ms
コンパイル使用メモリ 12,800 KB
実行使用メモリ 12,032 KB
最終ジャッジ日時 2024-07-19 07:49:30
合計ジャッジ時間 3,061 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 35 RE * 1
権限があれば一括ダウンロードができます

ソースコード

diff # 

import copy

def read_data():
    n, m = map(int, input().split())
    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 からゴールするまでのコストの期待値の最小値（魔法を１度使ってもよい）
    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