Number of Rubik's Cube states with a solved side. Uses the inclusion–exclusion principle. Posted at https://www.speedsolving.com/forum/showthread.php?20384-Probability-Thread&p=1025447&viewfull=1#post1025447 and some more explanation at https://www.speedsolving.com/forum/showthread.php?20384-Probability-Thread&p=1025268&viewfull=1#post1025268

3x3_states_with_solved_side.py

from math import factorial
from re import sub

def calc(expression):
    expression = expression.replace('^', '**')
    expression = expression.replace('/', '//')
    expression = sub(r'(\w+)!', r'factorial(\1)', expression)
    return eval(expression)

# Numbers for certain cases (always sideSubsets * variationInsideSolvedsides * variationInRest)
a1 = calc(' 6 * (4! ^ 2)     * (4! * 3^3 * 8! * 2^7 / 2)') # any one (side)
a2 = calc('12 * (2 * 3!)^2   * (2! * 3^1 * 5! * 2^4 / 2)') # adjacent two (sides)
o2 = calc(' 3 * (4! ^ 4)     * (4! * 2^3) / 2')            # opposite two (sides)
a3 = calc(' 8 * (2^3)        * (3! * 2^2) / 2')            # adjacent three (sides)
s3 = calc('12 * (2^3 * 3!^2) * (2)')                       # straight three (sides)
b4 = calc('12 * (1)          * (2)')                       # block four (sides)
r4 = calc(' 3 * (2^3)        * (1)')                       # ring four (sides)
a5 = calc(' 6 * (1)          * (1)')                       # any five (sides)
a6 = calc(' 1 * (1)          * (1)')                       # all six (sides)

# Multipliers
ma1 = 1
ma2 = 1 - (2 * ma1)
mo2 = 1 - (2 * ma1)
ma3 = 1 - (3 * ma1 + 3 * ma2)
ms3 = 1 - (3 * ma1 + 1 * mo2)
mb4 = 1 - (4 * ma1 + 5 * ma2 + 1 * mo2 + 2 * ma3 + 2 * ms3)
mr4 = 1 - (4 * ma1 + 4 * ma2 + 2 * mo2 + 4 * ms3)
ma5 = 1 - (5 * ma1 + 8 * ma2 + 2 * mo2 + 4 * ma3 + 6 * ms3 + 4 * mb4 + 1 * mr4)
ma6 = 1 - (6 * ma1 + 12 * ma2 + 3 * mo2 + 8 * ma3 + 12 * ms3 + 12 * mb4 + 3 * mr4 + 6 * ma5)
print('Multipliers:', ma1, ma2, mo2, ma3, ms3, mb4, mr4, ma5, ma6)

total = ma1 * a1 + \
        ma2 * a2 + \
        ma2 * o2 + \
        ma3 * a3 + \
        ma3 * s3 + \
        mb4 * b4 + \
        mr4 * r4 + \
        ma5 * a5 + \
        ma6 * a6
print('There are {:,} states with a solved side.'.format(total))
print('Probability:  About 1 / {:,}'.format(round(calc('12! * 8! / 2 * 2^11 * 3^7') / total)))