Discussion on: Daily Challenge #269 - Decompose n!

from math import log, factorial


def sundaram3(max_n):
    """ Returns a list of all primes under max_n """
    numbers = list(range(3, max_n + 1, 2))
    half = (max_n) // 2
    initial = 4
    for step in range(3, max_n + 1, 2):
        for i in range(initial, half, step):
            numbers[i - 1] = 0
        initial += 2 * (step + 1)
        if initial > half:
            return [2] + list([_f for _f in numbers if _f])


def decomp(n):
    """ Performs factorial prime factorization of n.
        Returns a factor: exponent dictionary.
    """
    return {p: sum(n // p**k for k in range(1, int(log(n, p)) + 1))
            for p in sundaram3(n + 1)}


def main():
    testcases = [5, 12, 14, 17, 22, 25]
    for case in testcases:
        factors_str = ' * '.join(f'{b}^{p}' for b, p in decomp(case).items())
        print(f"{case}! = {factors_str}")


if __name__ == '__main__':
    main()