factors = 20 upper = math.factorial(factors) divisors = range(2, factors+1) current = upper #repeatedly attempt to divide current number by prime factors ordered #from largest to smallest as long as the result has a remainder of 0 while True: found = False for p in reversed(divisors): c = current / p if c % p == 0: found = True current = c break if not found: break print 'divided by', p, 'got', current def squareDiff(x): s = range(1, x+1) sumSquares = sum([x*x for x in s]) squareSum = math.pow(sum(s),2) diff = squareSum - sumSquares return diff squareDiff(100) brute = lambda n: sum([x*x for x in xrange(1,n+1)]) poly = lambda n: round(1./6 * n + 1./2 * pow(n, 2) + 1./3 * pow(n, 3)) x = np.array(range(1,2200)) brute_y = np.array([brute(t) for t in x]) poly_y = np.array([poly(t) for t in x]) plt.plot(x, brute_y, label='Brute Force', color='blue') plt.plot(x, poly_y, label='Polynomial', color='red') plt.legend() print 'max error:', max(brute_y - poly_y) brute = lambda n: sum([1.*x*x for x in xrange(1,n+1)]) poly = lambda n: round(1./6 * n + 1./2 * pow(n, 2.) + 1./3 * pow(n, 3.)) x = np.array(range(1,2200)) brute_y = np.array([brute(t) for t in x]) poly_y = np.array([poly(t) for t in x]) plt.plot(x, brute_y, label='Brute Force', color='blue') plt.plot(x, poly_y, label='Polynomial', color='red') plt.legend() print 'max error:', max(brute_y - poly_y) def squareDiff2(x): sumSquares = round(1./6 * x + 1./2 * pow(x, 2.) + 1./3 * pow(x, 3.)) squareSum = pow(x*(x+1)/2., 2) return squareSum - sumSquares squareDiff2(100)