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))