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)