Efficient summation in Python

Question

I am trying to efficiently compute a summation of a summation in Python:

WolframAlpha is able to compute it too a high n value: sum of sum.

I have two approaches: a for loop method and an np.sum method. I thought the np.sum approach would be faster. However, they are the same until a large n, after which the np.sum has overflow errors and gives the wrong result.

I am trying to find the fastest way to compute this sum.

import numpy as np import time  def summation(start,end,func):     sum=0     for i in range(start,end+1):         sum+=func(i)     return sum  def x(y):     return y  def x2(y):     return y**2  def mysum(y):     return x2(y)*summation(0, y, x)  n=100  # method #1 start=time.time() summation(0,n,mysum) print('Slow method:',time.time()-start)  # method #2 start=time.time() w=np.arange(0,n+1) (w**2*np.cumsum(w)).sum() print('Fast method:',time.time()-start) 

Answer 1

Here's a very fast way:

result = ((((12 * n + 45) * n + 50) * n + 15) * n - 2) * n // 120 

How I got there:

1. Rewrite the inner sum as the well-known x*(x+1)//2. So the whole thing becomes sum(x**2 * x*(x+1)//2 for x in range(n+1)).
2. Rewrite to sum(x**4 + x**3 for x in range(n+1)) // 2.
3. Look up formulas for sum(x**4) and sum(x**3).
4. Simplify the resulting mess to (12*n**5 + 45*n**4 + 50*n**3 + 15*n**2 - 2*n) // 120.
5. Horner it.

Another way to derive it if after steps 1. and 2. you know it's a polynomial of degree 5:

3. Compute six values with a naive implementation.
4. Compute the polynomial from the six equations with six unknowns (the polynomial coefficients). I did it similarly to this, but my matrix A is left-right mirrored compared to that, and I called my y-vector b.

Code:

from fractions import Fraction import math from functools import reduce  def naive(n):     return sum(x**2 * sum(range(x+1)) for x in range(n+1))  def lcm(ints):     return reduce(lambda r, i: r * i // math.gcd(r, i), ints)  def polynomial(xys):     xs, ys = zip(*xys)     n = len(xs)     A = [[Fraction(x**i) for i in range(n)] for x in xs]     b = list(ys)     for _ in range(2):         for i0 in range(n):             for i in range(i0 + 1, n):                 f = A[i][i0] / A[i0][i0]                 for j in range(i0, n):                     A[i][j] -= f * A[i0][j]                 b[i] -= f * b[i0]         A = [row[::-1] for row in A[::-1]]         b.reverse()     coeffs = [b[i] / A[i][i] for i in range(n)]     denominator = lcm(c.denominator for c in coeffs)     coeffs = [int(c * denominator) for c in coeffs]     horner = str(coeffs[-1])     for c in coeffs[-2::-1]:         horner += ' * n'         if c:             horner = f"({horner} {'+' if c > 0 else '-'} {abs(c)})"     return f'{horner} // {denominator}'  print(polynomial((x, naive(x)) for x in range(6))) 

Output (Try it online!):

((((12 * n + 45) * n + 50) * n + 15) * n - 2) * n // 120