math.prod is slow at multiplying arbitrary-precision numbers. E.g., compare the run time of factorial(50000) to prod(range(2, 50001)).

factorial has some special-case optimizations, but the bulk of the difference is due to prod evaluating an expression tree of depth n. If you re-parenthesize the product so that the tree has depth log n, as factorial does, it's much faster. The evaluation order of prod isn't documented, so I think the change would be safe.

factorial uses recursion to build the tree, but it can be done iteratively with no advance knowledge of the total number of nodes.

This trick is widely useful for turning a way of combining two things into a way of combining many things, so I wouldn't mind seeing a generic version of it in the standard library, e.g. reduce(..., order='mid'). For many specific cases there are more efficient alternatives (''.join, itertools.chain, set.unions, heapq.merge), but it's nice to have a recipe that saves you the trouble of writing special-case algorithms at the cost of a log factor that's often ignorable.