The expected value is the weighted average of all possible outcomes of an experiment, weighted with the probabilities of each particular outcome. For a random variable ${X}$, the expected value is written as $E[X]$.  

Intuitively, the expected value is the long run average value of repetitions of the experiment.  

The variance is the expected value of the outcome's squared deviation from its expected value. For a random variable ${X}$, the variance is written as ${Var}[X]$ and is defined as the expected value of $(X-E[X])^2$.  

Intuitively, the variance is a measure of how far the outcomes of an experiment are spread out. The higher the variance, the more spread out the outcomes.  

Let's say we perform the following experiment involving throwing a die:  

Throw the die, and record the outcome as d[1].  

For i from 2 to N:
    Repeatedly throw the die until the outcome is different from d[i-1].
    Record the outcome as d[i].  

Output d[1] + d[2] + ... + d[N].

The die used in this experiment is a standard 6-sided die with outcomes $1,2,\ldots,6$. However, it is biased. In each throw, the probability of ${i}$ appearing is $p_i$ for $i=1,2,\ldots6$.  

Find the expected value and variance of the outcome of this experiment.  

Note: Certain formulas for variance are not fit for computation because of loss of significance/numerical instability. This link contains a discussion about how to avoid/mitigate this problem.

Input Format

The first six lines contain the probabilities of the die's outcomes. Specifically, the ${i}$th line contains $p_i$, for $i=1,2,\ldots6$. 

The seventh (and final) line contains $N$, the number of times the die is thrown.  

Constraints 

$0.1\leq p_i\leq0.2$ 

$p_1+p_2+\ldots+p_6=1$  

For test cases worth $25\%$ of the total points: $1\leq N\leq8$ 

For test cases worth $25\%$ of the total points: $1\leq N\leq3000$ 

For test cases worth $50\%$ of the total points: $1\leq N\leq100000$  

Output Format

The first line of output contains the expected value. 

The second line contains the variance. 

The answer will be accepted if it is within an absolute error of ${10^{-5}}$ of the true answer.  

Sample Input
0.16666666667
0.16666666666
0.16666666667
0.16666666667
0.16666666666
0.16666666667
2

Sample Output
7.0
4.66666666666

Explanation

One can verify these results by writing code that performs the experiment, running it multiple times, and computing the expected value and variance from the outcomes. The more times the experiment is run, the more accurate the answer will be.

Difficulty: MEDIUM
Tags: Probability

Starter code:

Solutions:
def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print(mean + n)
print(var)


def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print('%.7f' % (mean + n))
print('%.7f' % (var,))


def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print(mean + n)
print(var)


import math

def stepP(z, p):
	z = tuple((zi / (1 - pi) for (zi, pi) in zip(z, p)))
	c = sum(z)
	return tuple((c * pi - zi * pi for (zi, pi) in zip(z, p)))

def stepP_n(z, p, n=1):
	for _ in range(n):
		z = stepP(z, p)
	return z

def calc_pi(p):
	PI = tuple((pi * (1 - pi) for pi in p))
	c = sum(PI)
	PI = tuple((pii / c for pii in PI))
	return PI

def calc_z(p0, p, PI):
	z = tuple((p0i * (1 - pi) for (p0i, pi) in zip(p0, p)))
	c = sum(z)
	z = tuple((zi - c * pii for (zi, pii) in zip(z, PI)))
	return z

def dot(a, b):
	return sum((ai * bi for (ai, bi) in zip(a, b)))

def mul(a, b):
	return tuple((ai * bi for (ai, bi) in zip(a, b)))

def add(a, b):
	return tuple((ai + bi for (ai, bi) in zip(a, b)))

def sub(a, b):
	return tuple((ai - bi for (ai, bi) in zip(a, b)))

def E_sum_n(p0, p, n, x, eps=1e-05):
	if n <= 0:
		return 0
	if n == 1:
		return dot(p0, x)
	PI = calc_pi(p)
	sum_p0 = sum(p0)
	p0 = tuple((p0i - sum_p0 * pii for (p0i, pii) in zip(p0, PI)))
	stationary = n * sum_p0 * dot(PI, x)
	z = calc_z(p0, p, PI)
	transient = dot(z, x)
	norm_prod = max((abs(xi) for xi in x)) * sum((abs(zi) for zi in z))
	n0 = (math.log(norm_prod) - math.log(eps)) / math.log(3) if norm_prod > 0 else 0
	finite_n = 0 if n > n0 else dot(stepP_n(z, p, n), x)
	return stationary + transient - finite_n

def expectation(p, n):
	return E_sum_n(p, p, n, tuple(range(1, 7)), eps=1e-06)

def variance(p, n):
	x = tuple(range(1, 7))
	PI = calc_pi(p)
	e = dot(x, PI)
	x = tuple((i - e for i in x))
	var_sum = E_sum_n(p, p, n, tuple((xi ** 2 for xi in x)), eps=1e-06)
	var_sum -= E_sum_n(p, p, n, x, eps=1e-06) ** 2
	z = mul(x, PI)
	x_p = tuple(((1 - pi) * xi for (xi, pi) in zip(x, p)))
	z0 = stepP(z, p)
	stat_cov_sum = n * dot(z0, x_p) - E_sum_n(z0, p, n, x_p, eps=1e-06)
	z = calc_z(sub(p, PI), p, PI)
	transient_cov_sum = E_sum_n(p0=stepP(mul(z, x), p), p=p, n=n - 1, x=x, eps=1e-06)
	err = n * max((xi ** 2 for xi in x)) * sum((abs(zi) for zi in z)) * (1 / 3) ** n
	finite_n = 0
	if err > 1e-06:
		for k in range(1, n):
			z0 = stepP_n(z, p, n - k)
			z0 = mul(z0, x)
			finite_n += dot(stepP_n(z0, p, k), x)
	return var_sum + 2 * stat_cov_sum + 2 * transient_cov_sum - 2 * finite_n
p = tuple((float(input()) for _ in range(6)))
n = int(input())
print(expectation(p, n))
print(variance(p, n))


def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print(mean + n)
print(var)


def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print(mean + n)
print(var)


def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print(mean + n)
print(var)


def reduce_probs(p):
	probs = [[0 for x in range(6)] for y in range(6)]
	for x in range(6):
		for y in range(6):
			if x == y:
				continue
			probs[x][y] = p[y] / (1.0 - p[x])
	return probs

def calc_stats(n, probs):
	mean = [[0 for x in range(6)] for y in range(6)]
	var = [[0 for x in range(6)] for y in range(6)]
	if n == 1:
		for x in range(6):
			for y in range(6):
				mean[x][y] = y
		return (mean, var, probs)
	p = [[0 for x in range(6)] for y in range(6)]
	if n % 2 == 1:
		(Xm, Xv, Xp) = calc_stats(n - 1, probs)
		for x in range(6):
			for y in range(6):
				for m in range(6):
					p[x][y] += Xp[x][m] * probs[m][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					mean[x][y] += (Xm[x][m] + y) * Xp[x][m] * probs[m][y] / p[x][y]
					var[x][y] += Xv[x][m] * Xp[x][m] * probs[m][y] / p[x][y]
		for x in range(6):
			for y in range(6):
				for m in range(6):
					var[x][y] += (Xm[x][m] + y - mean[x][y]) ** 2 * Xp[x][m] * probs[m][y] / p[x][y]
		return (mean, var, p)
	(Xm, Xv, Xp) = calc_stats(n // 2, probs)
	for x in range(6):
		for y in range(6):
			for m in range(6):
				p[x][y] += Xp[x][m] * Xp[m][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				mean[x][y] += (Xm[x][m] + Xm[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
				var[x][y] += (Xv[x][m] + Xv[m][y]) * Xp[x][m] * Xp[m][y] / p[x][y]
	for x in range(6):
		for y in range(6):
			for m in range(6):
				var[x][y] += (Xm[x][m] + Xm[m][y] - mean[x][y]) ** 2 * Xp[x][m] * Xp[m][y] / p[x][y]
	return (mean, var, p)
p = [float(input()) for x in range(6)]
n = int(input())
if n == 1:
	mean = 0
	var = 0
	for x in range(6):
		mean += p[x] * x
	for x in range(6):
		var += (x - mean) ** 2 * p[x]
else:
	probs = reduce_probs(p)
	(Xm, Xv, Xp) = calc_stats(n - 1, probs)
	mean = 0
	var = 0
	for x in range(6):
		for y in range(6):
			mean += (x + Xm[x][y]) * p[x] * Xp[x][y]
			var += Xv[x][y] * p[x] * Xp[x][y]
	for x in range(6):
		for y in range(6):
			var += (Xm[x][y] + x - mean) ** 2 * p[x] * Xp[x][y]
print(mean + n)
print(var)


def mul(A, B):
	M = [[0] * 6 for x in range(6)]
	for i in range(6):
		for j in range(6):
			s = 0
			for k in range(6):
				s += A[i][k] * B[k][j]
			M[i][j] = s
	return M

def make_trans(p):
	P = [[0] * 6 for x in range(6)]
	for i in range(6):
		div = 1 - p[i]
		for j in range(6):
			if i != j:
				P[i][j] = p[j] / div
	return P

def countMU(M, p):
	mu = [0] * 6
	for i in range(6):
		for j in range(6):
			mu[i] += p[j] * M[j][i]
	return mu

def countEX(mu):
	s = 0
	for i in range(6):
		s += (i + 1) * mu[i]
	return s

def EY(N, k, EXL):
	if N <= k:
		return sum([EXL[i] for i in range(N)])
	return sum(EXL) + (N - k) * EXL[k - 1]

def EXX(n, m, muL, PL):
	if n == m:
		s = 0
		for i in range(6):
			s += muL[n][i] * xx[i]
		return s
	if n > m:
		(n, m) = (m, n)
	s = 0
	for i in range(6):
		si = 0
		for j in range(6):
			si += PL[m - n - 1][i][j] * x[j]
		s += x[i] * muL[n][i] * si
	return s

def EYY(N, k, muL, PL):
	XX = [[0] * k for i in range(k)]
	for i in range(k):
		for j in range(i, k):
			XX[i][j] = EXX(i, j, muL, PL) - EXL[i] * EXL[j]
			XX[j][i] = XX[i][j]
	if N <= k:
		S = 0
		for i in range(N):
			for j in range(N):
				S += XX[i][j]
		return S
	A = 0
	for i in range(k):
		for j in range(k):
			A += XX[i][j]
	B = 0
	for i in range(k - 1):
		B += XX[k - 1][i]
	return A + (N - k) * (2 * B + XX[k - 1][k - 1])
p = []
for i in range(6):
	p.append(float(input()))
N = int(input())
P = make_trans(p)
x = [i for i in range(1, 7)]
xx = [i * i for i in x]
M = mul(P, P)
PL = [P, M]
mu = countMU(P, p)
muL = [p, mu]
EXL = [countEX(p), countEX(mu)]
eps = 10 ** (-15)
k = 2
for i in range(N):
	M = mul(M, P)
	PL.append(M)
	mu = countMU(P, mu)
	muL.append(mu)
	EXL.append(countEX(mu))
	k += 1
	S = sum([abs(PL[k - 2][i][j] - PL[k - 1][i][j]) for i in range(6) for j in range(6)])
	if S < eps:
		break
print(EY(N, k, EXL))
print(EYY(N, k, muL, PL))