[VA] SRC#001 - Spiky Integral

[Albert Chan]

Using Python + numpy, I managed to calculate I(1000)

Code:

import numpy

def spike(n):
    terms = 1 + n*(n+1) // 2    # terms count of Product[z^i + 1/z^i, {i, n}]
    if terms % 2 == 0: return 0 # no constant term
    b = numpy.array([0] * (terms//2 + 1), dtype=object)
    b[-2:] = 1                  # 1 + z
    for i in range(2, n+1):     # multiply by 1 + z^i
        b[:-i] += b[i:]         # build half terms
    return b[0]                 # constant term

>>> spike.spike(1000)
46770858699105378013047692849647150249048020026391352159837485075187255449381044​
44575512200800626116981293400844929213580025428715543816767352635000639866305236​
49271451391806356678763700553308953903873511563867155215734010997332808966175771​
715884247295296277348179194597363883854664431808932677416L

Above calculation take 49 seconds on my laptop

I(1000) = b[0] / 2^1000 * (2 Pi) =

0.000274258153608378926807734432669808007979394749673091726358234027755841714506​
72423455696454538012082538178315765975675889323840397403322977190964502744011004​
81739611552026042903356881417709016110968635764441594831973619603002437175558485​
42910006760212326308258399120935101281203608176357009606892564575924775067719 ...

Rounded to 45 digits, it match Valentin value exactly

Is it possible to explain how the integration code work ?
How does it handle so many spikes ?