Initially I thought this was an impossible task. If placing cubelets one by one there are 27!/(8!*12!*6!) = 783,055,973,700 distinct orders in which types of cubelets could be arranged. Even dividing by the 6*4 = 24 ways that a single arrangement of cubes shows up in the counting, leaves more than 32,627,332,237.5 ways; the fraction resulting from some of the arrangements being symmetrical so as not to have 24 ways of showing up. That's a lot of calculation.
You can't just calculate the corner pieces, side pieces and face pieces (classified by their new positions) separately and multiply them together, as the probabilities are not independent, with the probability of the edge pieces being all-white-showing depends on the number of pieces of various types (from original position) remaining in the mix.
Then I realized the way around this problem is to keep a matrix of how many of each type remains at each stage and multiply the specific instances to calculate the next generation of the matrix (a 9x13x7x2 array, with the highest subscripts being the 8, 12, 6 and 1 vertex, edge, face and center pieces respectively and all having the possibility of a zero subscript for zero of that type remaining).
The binary distribution is used to find the likelihood of a given combination of V, E, F and C pieces at a given time. It needs be multiplied by the likelihood of having a given set available at any time (with the requisite feature of not having a black square showing up previously). And of course the probability that a given type in a given position will result in a success.
That latter probability is kept in a matrix of its own:
new type (original position)
placement vertex edge face center
V 1/8 3/12 1/2 1
E 3/12 5/12 2/3 1
F 1/2 2/3 5/6 1
C 1 1 1 1
The program below reports:
phase 1
phase 2
1.20158738370485E-06
phase 3
2.15921476422768E-15
final tally
4.92046901590714E-18 = 1/2.03232658668747E+17
The 1.20158738370485E-06 indicates that after phase 1, the assignment of vertex pieces, the overall probability of success is about 1.2016 x 10^-6.
Likewise after that the probability that success continues through phase 2 is about 2.1592 x 10^-15.
The overall probability of going without any black squares showing up on the outside is about 4.9205 x 10^-18, or 1 in 2.0323 x 10^17, or 1 in 203,232,658,668,747,000 (remember, that's about; the three zeros at the end are not part of the significant figures).
There are a lot of places for possible error here: first in the theory, second in the calculation of the above probabilities, but thirdly in the programming. I put in debugging statements (wherever you see xx=xx was a point for pausing the program to look for something wrong), and went over the listing to look for places where, say, I didn't change a 19 to a 7 in the section I copied code from the edge portion to the face portion. So other sets of eyes would be helpful.
DefDbl A-Z
Dim crlf$, prob(3, 3)
Private Sub Form_Load()
Form1.Visible = True
Text1.Text = ""
crlf = Chr$(13) + Chr$(10)
' prob(a,b) is that of all-white showing at position a
' when taking a piece from the original at position b
' where positions are:
' 0 = corner (vertex)
' 1 = edge
' 2 = face center
' 3 = cube center
prob(0, 0) = 1 / 8: prob(0, 1) = 3 / 12: prob(0, 2) = 1 / 2: prob(0, 3) = 1
prob(1, 0) = 3 / 12: prob(1, 1) = 5 / 12: prob(1, 2) = 2 / 3: prob(1, 3) = 1
prob(2, 0) = 1 / 2: prob(2, 1) = 2 / 3: prob(2, 2) = 5 / 6: prob(2, 3) = 1
prob(3, 0) = 1: prob(3, 1) = 1: prob(3, 2) = 1: prob(3, 3) = 1
' phase 1: transition to states after selection of four corners
Text1.Text = Text1.Text & "phase 1" & crlf
ReDim newState(8, 12, 6, 1)
For vPieces = 0 To 8
remain = 8 - vPieces
For ePieces = 0 To remain
remain = remain - ePieces
For fPieces = 0 To remain
If fPieces <= 6 Then
remain = remain - fPieces
cPieces = remain
If cPieces <= 1 And cPieces >= 0 Then
p = prob(0, 0) ^ vPieces * prob(0, 1) ^ ePieces * prob(0, 2) ^ fPieces
newV = 8 - vPieces
newE = 12 - ePieces
newF = 6 - fPieces
newC = 1 - cPieces
p = p * binom(vPieces, 8, 27) * binom(ePieces, 12, 27 - vPieces) * binom(fPieces, 6, 27 - ePieces - vPieces)
If newV + newE + newF + newC <> 19 Then
xx = xx
End If
newState(newV, newE, newF, newC) = newState(newV, newE, newF, newC) + p
End If
remain = remain + fPieces
End If
Next fPieces
remain = remain + ePieces
Next ePieces
remain = remain + vPieces
Next vPieces
' phase 2 choose edge pieces for each possible availability
' of numbers of pieces
Text1.Text = Text1.Text & "phase 2" & crlf
ReDim oldState(8, 12, 6, 1)
For a = 0 To 8
For b = 0 To 12
For c = 0 To 6
For d = 0 To 1
oldState(a, b, c, d) = newState(a, b, c, d)
If newState(a, b, c, d) > 0 And a + b + c + d <> 19 Then
xx = xx
End If
tot1 = tot1 + newState(a, b, c, d)
Next
Next
Next
Next
Text1.Text = Text1.Text & tot1 & crlf
ReDim newState(8, 12, 6, 1)
For a = 0 To 8
For b = 0 To 12
For c = 0 To 6
For d = 0 To 1
If a + b + c + d = 19 Then
For vPieces = 0 To 12
If vPieces <= a Then
remain = 12 - vPieces
For ePieces = 0 To remain
If ePieces <= b Then
remain = remain - ePieces
For fPieces = 0 To remain
If fPieces <= c Then
remain = remain - fPieces
cPieces = remain
If cPieces <= d And cPieces >= 0 Then
p = prob(1, 0) ^ vPieces * prob(1, 1) ^ ePieces * prob(1, 2) ^ fPieces
newV = a - vPieces
newE = b - ePieces
newF = c - fPieces
newC = d - cPieces
p = p * binom(vPieces, a, 19) * binom(ePieces, b, 19 - vPieces) * binom(fPieces, c, 19 - ePieces - vPieces)
p = p * oldState(a, b, c, d)
If newV + newE + newF + newC <> 7 Then
xx = xx
End If
newState(newV, newE, newF, newC) = newState(newV, newE, newF, newC) + p
End If
remain = remain + fPieces
End If
Next fPieces
remain = remain + ePieces
End If
Next ePieces
remain = remain + vPieces
End If
Next vPieces
End If
Next
Next
Next
Next
' phase 3 choose face pieces for each possible availability
' of numbers of pieces
Text1.Text = Text1.Text & "phase 3" & crlf
ReDim oldState(8, 12, 6, 1)
For a = 0 To 8
For b = 0 To 12
For c = 0 To 6
For d = 0 To 1
oldState(a, b, c, d) = newState(a, b, c, d)
If newState(a, b, c, d) > 0 And a + b + c + d <> 7 Then
xx = xx
End If
tot2 = tot2 + newState(a, b, c, d)
Next
Next
Next
Next
Text1.Text = Text1.Text & tot2 & crlf
ReDim newState(8, 12, 6, 1)
For a = 0 To 8
For b = 0 To 12
For c = 0 To 6
For d = 0 To 1
If a + b + c + d = 7 Then
For vPieces = 0 To 6
If vPieces <= a Then
remain = 6 - vPieces
For ePieces = 0 To remain
If ePieces <= b Then
remain = remain - ePieces
For fPieces = 0 To remain
If fPieces <= c Then
remain = remain - fPieces
cPieces = remain
If cPieces <= d And cPieces >= 0 Then
p = prob(2, 0) ^ vPieces * prob(2, 1) ^ ePieces * prob(2, 2) ^ fPieces
newV = a - vPieces
newE = b - ePieces
newF = c - fPieces
newC = d - cPieces
p = p * binom(vPieces, a, 7) * binom(ePieces, b, 7 - vPieces) * binom(fPieces, c, 7 - ePieces - vPieces)
p = p * oldState(a, b, c, d)
If newV + newE + newF + newC <> 1 Then
xx = xx
End If
newState(newV, newE, newF, newC) = newState(newV, newE, newF, newC) + p
End If
remain = remain + fPieces
End If
Next fPieces
remain = remain + ePieces
End If
Next ePieces
remain = remain + vPieces
End If
Next vPieces
End If
Next
Next
Next
Next
'final tally
Text1.Text = Text1.Text & "final tally " & crlf
For a = 0 To 8
For b = 0 To 12
For c = 0 To 6
For d = 0 To 1
If newState(a, b, c, d) > 0 And a + b + c + d <> 1 Then
xx = xx
End If
tot3 = tot3 + newState(a, b, c, d)
Next
Next
Next
Next
Text1.Text = Text1.Text & tot3 & " = 1/" & 1 / tot3 & crlf
Text1.Text = Text1.Text & crlf & " done"
End Sub
Function binom(a, b, c)
p = 1
For i = 0 To a - 1
p = p * (b - i) / (c - i)
Next
binom = p
End Function
|
Posted by Charlie
on 2017-07-14 17:05:24 |