import math
from math import comb as nCk
from math import factorial as F

#1s, 2s, 3s - four of each card.
#ones, twos, threes - 14,13,13 times each one is spoken outloud.  Slots.

#Coefficients.
coef = [1,4,6,4,1]



#Place (a) 1s in twos and (4-a) 1s in threes.
a_all = 0
for a in range(5):
    
    #There are 13 twos to place the 1s in.
    n2 = 13
    #There are 13 threes to place the 1s in.
    n3 = 13
    
    #There are coef[a] ways to order four 1s such that (a) are in the twos and (4-a) are in the threes.
    #a_ is the total number of ways for the product.
    a_ = coef[a] * nCk(n2, a) * nCk(n3, 4-a)
    #Twos remaining.
    n2 = n2 - a
    #Threes remaining.
    n3 = n3 - (4-a)
    
    #Place (b) 2s in ones and (4-b) 2s in threes.
    b_all = 0
    for b in range(5):
        
        #There are 14 ones to place the 2s in.
        n1 = 14
        
        #There are coef[b] ways to order four 2s such that (b) are in the ones and (4-b) are in the threes.
        #b_ is the total number of ways for the product.
        b_ = coef[b] * nCk(n1, b) * nCk(n3, 4-b)
        #Ones remaining.
        n1 = n1 - b
        
        #Place (c) 3s in ones and (4-c) 3s in twos.
        c_all = 0
        for c in range(5):
            
            #There are coef[c] ways to order four 3s such that (c) are in the ones and (4-c) are in the twos.
            #c_ is the total number of ways for the product.
            c_ = coef[c] * nCk(n1, c) * nCk(n2, 4-c)
            
            #All five (c) distributions for a particular (a) and (b). 
            c_all += c_
        b_ = b_ * c_all

        #All five (b) distributions for a particular (a). 
        b_all += b_
    a_ = a_ * b_all

    #All five (a) distributions yields the total.
    a_all += a_
print(a_all)

149967879125350

149967879125350* F(28) / F(40)

5.603946800634387e-05

Rohan Lewis

2023.02.20