First of all the binom function was misnamed, though it did the correct multiplications (choices without replacement).  But after those multiplications, the permutations should have been counted.

For example, if we were choosing the first 8 cubelets and designating them for the 8 vertices and specifically for the case of say choosing 3 vertex pieces, 2 edge pieces, 2 face pieces and the center piece, it was calculating the probability that those would be the resulting numbers as

(8/27 * 7/26 * 6/25) * (12/24 * 11/23) * (6/22 * 5/21)

But that should be multiplied by 8!/(3!*2!*2!), as the order of the choice does not matter, so all orders must count.

The new result:

phase 1

phase 2

4.91668896501502E-03

phase 3

3.69404873880623E-08

final tally 

2.23709604269576E-10  = 1/4470080769.50946

The added lines of code were

           p = p * fact(8) / (fact(vPieces) * fact(ePieces) * fact(vPieces))

and

           p = p * fact(12) / (fact(vPieces) * fact(ePieces) * fact(vPieces))

and

           p = p * fact(6) / (fact(vPieces) * fact(ePieces) * fact(vPieces))

in the sections on choosing new vertices, new edges and new faces respectively, and of course incorporating the fact() function.

The new listing is:

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)

           p = p * fact(8) / (fact(vPieces) * fact(ePieces) * fact(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 * fact(12) / (fact(vPieces) * fact(ePieces) * fact(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 * fact(6) / (fact(vPieces) * fact(ePieces) * fact(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

Function fact(n)

   f = 1

   For i = 2 To n

     f = f * i

   Next

   fact = f

End Function