[VA] SRC#001 - Spiky Integral

[Gerson W. Barbosa]

(07-12-2018 10:32 AM)Pjwum Wrote:  With HP Prime

we can go beyond 10

0
0
1/2 pi
1/4 pi
0
0
1/8 pi
7/64 pi
0
0
35/512 pi
31/512 pi
0
0
361/8192 pi
657/16384 pi
0
0
2055/65536 pi
1909/65536 pi

Sure we can:

%%HP: T(3)A(R)F(.);
\<< DUPDUP 'X' 1 ROT OVER SWAP
  FOR i OVER i ^ DUP INV + *
  NEXT NIP EXPAND FXND DROP \->STR "*X^" ROT DUPDUP * + 2 / DUP 2 + \->STR "+" + UNROT \->STR + "+" + "X^" ROT + PICK3 SWAP POS PICK3 ROT POS 1 - SUB DUP SIZE 1 -
  IF NOT
  THEN DROP2 0
  ELSE DUP SIZE OVER "+" POS 1 + SWAP SUB OBJ\-> 2 ROT 1 - ^ / \pi *
  END
\>>

'VASRC1' STO

40

%%HP: T(3)A(R)F(.);
\<< { } 1 ROT
  FOR i i VASRC1 +
  NEXT
\>>

EVAL

-->

{ 0 0 '1/2*π' '1/4*π' 0 0 '1/8*π' '7/64*π' 0 0 '35/512*π' '31/512*π' 0 0 '361/8192*π' '657/16384*π' 0 0 '2055/65536*π' '1909/65536*π' 0 0 '24955/1048576*π' '46923/2097152*π' 0 0 '316301/16777216*π' '299973/16777216*π' 0 0 '4136805/268435456*π' '15796439/1073741824*π' 0 0 '13853361/1073741824*π' '26585247/2147483648*π' 0 0 '756388295/68719476736*π' '182188585/17179869184*π' }

Not the best method, I fear. 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\). The cases where the results are zero should be handled more cleverly, as you have suggested, but this is only a test. The RPL program might not be fast enough on the real HP 50g as N get larger. For N = 12, it takes 16.75 seconds; for N = 20 it takes 118.41 seconds. The evaluation of the integrals would take much, much longer, I guess.

Gerson.