The idea is simple. Notice that given integers n and t, where t >= 2, the product of x1~xt (sum x1~xt = n, xi > 0, xi is an integer) reaches its maximum when xi are (almost) uniform. For example, when n = 10 and t = 3, the maximum product is given by [3,3,4]. Thus, iterating over all possible t gives the final result.
return max([reduce(lambda x,y: x*y, [int(n / t) + (1 if i < n % t else 0) for i in range(t)]) for t in range(2, n+1)])