(los métodos más rápidos, 3 y 4, están al final)

En un método rápido de NumPy, debe especificar dtype=np.object para que NumPy no convierta Python int en sus propios dtypes ( np.int64 u otros). Ahora le dará resultados correctos (lo verificó hasta N = 100000).

 # method #2 start=time.time() w=np.arange(0, n+1, dtype=np.object) result2 = (w**2*np.cumsum(w)).sum() print('Fast method:', time.time()-start)

Su solución rápida es significativamente más rápida que la lenta. Sí, para N grandes, pero ya en N=100 es como 8 veces más rápido:

 start=time.time() for i in range(100): result1 = summation(0, n, mysum) print('Slow method:', time.time()-start) # method #2 start=time.time() for i in range(100): w=np.arange(0, n+1, dtype=np.object) result2 = (w**2*np.cumsum(w)).sum() print('Fast method:', time.time()-start)

 Slow method: 0.06906533241271973 Fast method: 0.008007287979125977

EDITAR: un método aún más rápido (por KellyBundy, la calabaza) es usar Python puro. Resulta que NumPy no tiene ventaja aquí, porque no tiene código vectorizado para np.objects .

 # method #3 import itertools start=time.time() for i in range(100): result3 = sum(x*x * ysum for x, ysum in enumerate(itertools.accumulate(range(n+1)))) print('Faster, pure python:', (time.time()-start))

 Faster, pure python: 0.0009944438934326172

EDIT2: Forss notó que el método rápido numpy se puede optimizar usando x*x en lugar de x**2 . Para N > 200 es más rápido que el método Python puro. Para N < 200 es más lento que el método puro de Python (el valor exacto del límite puede depender de la máquina, en el mío fue 200, es mejor que lo verifique usted mismo):

 # method #4 start=time.time() for i in range(100): w = np.arange(0, n+1, dtype=np.object) result2 = (w*w*np.cumsum(w)).sum() print('Fast method x*x:', time.time()-start)

Aquí hay una manera muy rápida:

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

Cómo llegué allí:

1. Reescriba la suma interna como el bien conocido x*(x+1)//2 . Entonces todo se convierte en sum(x**2 * x*(x+1)//2 for x in range(n+1)) .
2. Reescriba a sum(x**4 + x**3 for x in range(n+1)) // 2 .
3. Busque fórmulas para sum(x**4) y sum(x**3) .
4. Simplifica el desorden resultante a (12*n**5 + 45*n**4 + 50*n**3 + 15*n**2 - 2*n) // 120 .
5. Hornéalo .

Otra forma de derivarlo si después de los pasos 1 y 2 sabes que es un polinomio de grado 5:

3. Calcule seis valores con una implementación ingenua.
4. Calcule el polinomio a partir de las seis ecuaciones con seis incógnitas (los coeficientes del polinomio). Lo hice de manera similar a this , pero mi matriz A está reflejada de izquierda a derecha en comparación con eso, y llamé a mi vector y b .

Código:

 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)))

Salida ( ¡Pruébelo en línea! ):

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

Comparar Python con WolframAlpha de esa manera es injusto, ya que Wolfram simplificará la ecuación antes de calcular.

Afortunadamente, el ecosistema de Python no conoce límites, por lo que puedes usar SymPy :

 from sympy import summation from sympy import symbols n, x, y = symbols("n,x,y") eq = summation(x ** 2 * summation(y, (y, 0, x)), (x, 0, n)) eq.evalf(subs={"n": 1000})

Calculará el resultado esperado casi instantáneamente: 100375416791650 . Esto se debe a que SymPy simplifica la ecuación para usted, tal como lo hace Wolfram. Ver el valor de eq :

La respuesta de @Kelly Bundy es increíble, pero si eres como yo y usas una calculadora para calcular 2 + 2 , entonces te encantará SymPy ❤. Como puede ver, obtiene los mismos resultados con solo 3 líneas de código y es una solución que también funcionaría para otros casos más complejos.

Todas las respuestas usan matemáticas para simplificar o implementar el ciclo en Python tratando de ser óptimos para la CPU, pero no son óptimos para la memoria.

Aquí una implementación ingenua sin usar ninguna simplificación matemática que sea eficiente en memoria

 def function5(): inner_sum = float() result = float() for x in range(0, n + 1): inner_sum += x result += x ** 2 * inner_sum return result

Es bastante lento con respecto a las otras soluciones de dankal444:

 method 2 | 31 µs ± 2.06 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) method 3 | 116 µs ± 538 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) method 4 | 91 µs ± 356 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) function 5 | 217 µs ± 1.14 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

por cierto, si deja la función con numba (puede haber mejores opciones):

 from numba import jit function5 = jit(nopython=True)(function5)

usted obtiene

 59.8 ns ± 0.209 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

En un comentario, mencionas que en realidad son f(x) y g(y) en lugar de x ² e y. Si solo necesita una aproximación a esa suma, puede pretender que las sumas son sumas de Riemann de punto medio, de modo que su suma se aproxime mediante la integral doble ∫ _-.5 ^n+.5 f(x) ∫ _-.5 ^x+.5 g(y) dy dx.

Con su original f(x)=x ² y g(y)=y, esto se simplifica a n ⁵ /10+3n ⁴ /8+n ³ /2+5n ² /16+3n/32+1/160, que difiere del resultado correcto en n ³ /12+3n ² /16+53n/480+1/160.

Basado en esto, sospecho que (real-integral)/real sería max(f'',g'')*O(n ^-2 ), pero no pude probarlo.