def MDSsubd(t):
# starting with a degee 2 vertex blue; u1 and u2 are the two centers
d=t+1 # degree of HD vertices u1 and u2
# states are (0),(1) and (2,k1,k2) for k1,k2 = 0,...,t; (0) means neither u1 nor u2 blue
# (1) means one of u2,u1 blue and other white; for book-keeping u2 blue
# states (2,k1,k2) means u1,u2 blue and k1,k2 are the # of blue leaves
m=zero_matrix(RR,(t+1)^2+2) # m is the matrix; use QQ for exact, RR for decimal
# matrix row 0 is neither of u1,u2 blue = state (0)
# matrix row 1 is 1 of u2,u1 blue (say u2) = state (1)
# for both u1 and u2 blue, use lex order for states (2,k1,k2), i.e.,
# matrix row 2+k1(t+1)+k2 is state (2,k1,k2)
#
# forcing with both u1 and u2 white
# row 0 = state (0): the degree 2 blue vertex can force 0, 1, or 2 of u1 and u2
m[0,0]=1/4
m[0,1]=1/2
m[0,2]=1/4
# forcing with u2 blue and u1 white
# u2 blue implies u1 will become blue in next round, forced by the degree 2 vertex w
pfu=2/d # probability of u2 forcing any one specific leaf neighbor
pnfu=1-pfu # probability of u not forcing any one specific leaf neighbor
for r in [0..t]: # r is the number of new blue leaves
m[1,2+r]=binomial(t,r)*pfu^r*pnfu^(t-r)
# starting with u1 and u2 both blue
# state (2,k1,k2) is matrix index 2+k1(t+1)+k2 for k1,k2=0,...,t
prm=zero_matrix(QQ,t+1) # prm is the leaf probability matrix; use QQ for exact, RR for decimal
# prm[i,j] is the probability of going from i blue leaves to j blue leaf neighbprs of u2 (or u1)
for i in [0..t-2]: # construct prm
pfu=(i+2)/d # probability center forcing one specific leaf when i leaf neighbors are blue
pnfu=1-pfu # probability center not forcing any one specific leaf neighbor when i are blue
for j in [i..t]:
x=t-i # x is # white vertices
r=j-i # r is the number of new blue leaves
prm[i,j]=binomial(x,r)*pfu^r*pnfu^(x-r) # prob r additional blue u-leaves
prm[t-1,t]=1
prm[t,t]=1
# construct transition matrix m
for k1 in [0..t]:
for k2 in [0..t]:
for r1 in [0..t-k1]:
for r2 in [0..t-k2]:
m[2+k1*(t+1)+k2,2+(k1+r1)*(t+1)+k2+r2]=prm[k1,k1+r1]*prm[k2,k2+r2]
return m