polynomial-theorem.py

import sys
from functools import reduce

N = 4
VARS = 'abcd'

mult = lambda a, b: a * b

def showN(i):
	if i == 1:
		return None
	else:
		return str(i)

def showPow(a, b):
	if b == 0:
		return None
	elif b == 1:
		return a
	else:
		return f'{a}^{b}'

def fact(i):
	return reduce(mult, range(1, i + 1), 1)

def compute(n, i, v):
	if i < 0 or v == 0:
		return

	if v == 1:
		yield [i]
		return

	for j in reversed(range(n + 1)):
		for u in compute(n, i - j, v - 1):
			yield [j] + u

print('(' + ' + '.join(VARS) + ')^' + str(N), end='')
print(' = ', end='')
print(' + '.join([
	(
		'*'.join(filter(None, [
			showN(fact(N) // reduce(mult, map(fact, u))),
		] + [
			showPow(VARS[i], u[i]) for i in range(len(VARS))
		]))
	) for u in compute(N, N, len(VARS))
]))

# print([*compute(N, N, VARS)])