[VA] SRC#001 - Spiky Integral

[Gerson W. Barbosa]

(07-14-2018 08:38 PM)Thomas Klemm Wrote:  These are the examples for \(N=39\) and \(N=71\):

Code:

>>> b = spiky(100)
>>> print b[39][0], '/', 2**37
1512776590 / 137438953472
>>> print b[71][0], '/', 2**69
2681644149792639400 / 590295810358705651712

However the ratios haven't been reduced.

In exact mode, the HP 50g does it automatically.

(07-14-2018 08:38 PM)Thomas Klemm Wrote:  

(07-12-2018 08:10 PM)Gerson W. Barbosa Wrote:  Expand Product{ k=1..n, x^N + 1/x^N } and take the coefficient of the power of x corresponding to the Nth triangular number in the numerator (if there is no correspondence, then the result will be zero). That's your numerator. Your denominator is 2^(N - 1). Multiply the resulting fraction by \(\pi\).

Not sure if I understood that correctly but it might explain the expression:
\[x^N+\frac{1}{x^N}\]

I haven't trodden any of the hard (and beautiful) paths you all have done. Instead, I took an unallowed (according to Valentin's rules) and dull shortcut. It was not difficult to recognize pi/2 and pi/4 in his his three-digit results for n = 3 and n = 4. Then I evaluated the integrals for n up to 12 on my CASIO fx-991 LA X, which doesn't take too long to return exact results in terms of pi for small values of n (1/8*pi after 3m 30s, for n = 7; 0.3436116965 after 4m 31s, for n = 8). The pattern soon became apparent: a fraction of pi, the denominator being a power of 2. For denominators = 2^(n - 1), the first numerators, starting with n = 1, are 0, 0, 2, 2, 0, 0, 8, 14, 0, 0, 70, 124, 0..., that is, OEIS sequence A063865.


Quoting from the formula section:


"a(n) = constant term in expansion of Product_{ k = 1..n } (x^k + 1/x^k). - N. J. A. Sloane, Jul 07 2008"

Thus, since

(X + 1/X)(X^2 + 1/X^2)(X^3 + 1/X^3)(X^4 + 1/X^4)(X^5 + 1/X^5)(X^6 + 1/X^6)(X^7 + 1/X^7)(X^8 + 1/X^8)

=

X^36 + 1/X^36 + X^34 + 1/X^34 + X^32 + 1/X^32 + 2 X^30 + 2/X^30 + 2 X^28 + 2/X^28 + 3 X^26 + 3/X^26 + 4 X^24 + 4/X^24 + 5 X^22 + 5/X^22 + 6 X^20 + 6/X^20 + 7 X^18 + 7/X^18 + 8 X^16 + 8/X^16 + 9 X^14 + 9/X^14 + 10 X^12 + 10/X^12 + 11 X^10 + 11/X^10 + 12 X^8 + 12/X^8 + 13 X^6 + 13/X^6 + 13 X^4 + 13/X^4 + 13 X^2 + 13/X^2 + 14

a(8) = 14

So, the numerator can be found by means of a polynomial expansion. However, when expanding this polynomial on the HP 50g, using the EXPAND command, I got

'(X^72+X^70+X^68+2*X^66+2*X^64+3*X^62+4*X^60+5*X^58+6*X^56+7*X^54+8*X^52+9*X^50+​10*X^48+11*X^46+12*X^44+13*X^42+13*X^40+13*X^38+14*X^36+13*X^34+13*X^32+13*X^30+12*X^28+11*X^26+10*X^24+9*X^22+8*X^20+7*X^18+6*X^1​6+5*X^14+4*X^12+3*X^10+2*X^8+2*X^6+X^4+X^2+1)/X^36',

which is equivalent to the previous polynomial, except that the constant term in the numerator is 1. But notice 14 is the coefficient in the numerator that corresponds to the power of the denominator (...+14*X^36+.../X^36). Also, 36 = 8*(8 + 1)/2, that is, the 8th triangular number. This has worked for some othe values of n I tried, so I assumed it is a valid property.

The degree of the numerator polynomial is n*(n + 1), which means its size grows quadratically as n increases, which is both memory and time-consuming. Hopefully your method of generating the coefficients is faster (since I don't know python, I can't properly decode your algorithm).

I am disappointed there isn't a formula to directly compute the terms of the sequence. There is an asymptotic formula at OEIS, but it is not good enough for practical purposes. So I made some adjustments, which are by no means exact, but might give 5 or 6 correct significant digits for relatively low n, when computing the integral:

a(n) ~ sqrt(6/pi)*2^n*(1-6/(5*n)+21/(20*n^2)-1/(8*n^3)+3/n^4)/(n*sqrt(n))

or

a(n) ~ (sqrt(3/p) 2^(n - 5/2) (n (2 n (4 n (5 n - 6) + 21) - 5) + 120))/(5 n^(11/2))

The following HP-42S program computes the integral for n >= 1 and n up to about 20000 (on Free42) and returns exact results for n < 15.

Code:


00 { 90-Byte Prgm }

01▸LBL "INCS"
02 ENTER
03 ENTER
04 ENTER
05 5
06 ×
07 6
08 -
09 ×
10 4
11 ×
12 21
13 +
14 ×
15 2
16 ×
17 5
18 -
19 3
20 PI
21 ÷
22 SQRT
23 ×
24 X<>Y
25 4.5
26 Y↑X
27 5
28 ×
29 ÷
30 X<>Y
31 X↑2
32 RCL+ ST L
33 2
34 ÷
35 1
36 -
37 2
38 MOD
39 ×
40 X<>Y
41 2.5
42 -
43 2
44 X<>Y
45 Y↑X
46 ×
47 0.5
48 +
49 IP
50 2
51 R↑
52 1
53 -
54 Y↑X
55 ÷
56 PI
57 ×
58 END

Gerson.