[VA] SRC#001 - Spiky Integral

[Thomas Klemm]

(07-18-2018 05:32 PM)ijabbott Wrote:  Is there a neat formula for just the constant term when converting the product to a sum?

From A058377:

Quote:FORMULA
a(n) is half the coefficient of q^0 in product('(q^(-k)+q^k)', 'k'=1..n) for n >= 1.

Thus I doubt there is a "neat formula".
However here's a table of n, a(n) for n = 1..3342

From this we can calculate the value of I(1000) in accordance with Valentin's result:

Code:

>>> pi * 233854293495526890065238464248235751245240100131956760799187425375936277246​90522222877561004003130584906467004224646067900127143577719083836763175003199331​52618246357256959031783393818502766544769519367557819335776078670054986664044830​87885857942123647648138674089597298681941927332215904466338708 / 2**998
0.00027425815360837894

Since I had trouble with the Python program with bigger numbers I implemented it in Clojure:

Code:

(defn zeros [k]
  (replicate k 0))

(defn spiky [n]
  (loop [k 2 a '(0 1)]
    (if (> k n) (first a) 
        (let [s (count a)
              [b c] (split-at k a)
              u (concat (zeros k) a)
              v (concat (conj (reverse b) 0) (zeros (dec s)))
              w (concat c (zeros (* 2 k)))] 
          (recur (inc k) (map +' u v w))))))

You may notice similarities to the implementation in RPL.

It's fast for n=100, takes a couple of seconds for n=200 and multiple minutes for n=1000.
I've tried it for n=2000 and after a while that felt like eternity the correct result was given.
From these measurements I would assume that it would take days or even weeks to calculate the value for n=20,000. Thus I refrained from trying.

Kind regards
Thomas