[VA] SRC#001  Spiky Integral

07102018, 10:10 PM
Post: #1




[VA] SRC#001  Spiky Integral
.
Hi all, welcome to my SRC#001  Spiky Integral: Here I'll deal with a realworld math problem not of my own making but which did appear at a certain math competition. The problem introduces this "spiky" integral: I(N) = \(\int_{0}^{2\pi} cos(x)\, cos(2x)\, cos(3x)\, ...\, cos(Nx)\, dx \) and asks for which values of N in the interval [1,10] does I(N) have a nonzero value. It then specifies that the result must include proof of correctness but in my extensive experience with problems of this sort I've found that it's best to first obtain the correct result using whatever means (praying included, bribing not excluded), as this will usually provide a most helpful "hint" to afterwards try and get proof of its validity. With this strategy in mind (and if praying and bribing have already proved ineffective), we can use our beloved HP calcs (yes, HP calcs, not Excel, not Mathematica, not Wolfram Alpha, not Maple, not Matlab, not laptops, you get the drift) to get some numerical evidence first, then use it to make an educated conjecture on what the correct result might probably be. As for numerical evidence, this is how I'd obtain it using an HP15C and a very simple, straightforward little program I wrote which goes as follows: The program segment which computes Cos(x)*Cos(2x)*Cos(3x)*...*Cos(Nx), where N is assumed to be in R0 and x in stack register X (the display), is placed under LBL B and the main program which computes the integral I(N) for N from 1 to 10 is placed after it under LBL A: 001 LBL B 015 LBL A 002 STO 1 016 1.01 003 RCL 0 020 STO 0 004 INT 021 LBL 0 005 STO I 022 0 006 1 023 PI 007 LBL 1 024 2 008 RCL I 025 x 009 RCLx 1 026 INTEG B 010 COS 027 RCL 0 011 x 028 R/S 012 DSE I 029 X<>Y 013 GTO 1 030 R/S 014 RTN 031 ISG 0 032 GTO 0 To run it, simply press: [USER] [RAD] [FIX 3] [A] > 1.010 [R/S] 1.083 04 [R/S] > 2.010 [R/S] 5.974 05 [R/S] > 3.010 [R/S] 1.571 [R/S] > 4.010 [R/S] 0.785 [R/S] > 5.010 [R/S] 3.778 08 [R/S] > 6.010 [R/S] 2.733 09 [R/S] > 7.010 [R/S] 0.393 [R/S] > 8.010 [R/S] 0.344 [R/S] > 9.010 [R/S] 3.600 04 [R/S] > 10.010 [R/S] 1.664 09 Notes:  It's best to run this on much faster, modern HP15Cbased hardware or a fast software emulator as otherwise it will take several hours for a physical vintage HP15C to get the results above, mostly because the function being integrated gets more an more "spiky" as N grows and so it takes more and more time to integrate it accurately.  FIX 3 is recommended for speed as FIX 4 or higher takes significantly longer and the extra accuracy isn't needed here. On the other hand, FIX 2 can't be used in the HP15C or the HP41C/Advantage versions because for N=4 we get: I(4) = 0.86 + 0.03 which is wrong, the uncertainty is lying. This is similar to what Mr. Kahan said about one of his sample integrals in some issue of the HP Journal: the error is too small to notice but too big to ignore. The solution is of course to use FIX 3 as above or SCI 2, which also gives the correct result, namely: I(4) = 7.85E01 + 1.48E03 With this evidence in mind and regarding anything smaller than 1E3 as 0 (as we are computing the integrals in FIX 3) we can see that I(N) is nonzero for N=3, 4, 7 and 8. Using this numerical evidence I have my own conjecture on what the result will be for general N, which I'll post in a few days together with my verbatim conversion of this HP15C program for the HP41C/Advantage ROM (which runs very fast in modern hardware/emulators) as well as a completely different 6line enhanced version for the HP71B/Math ROM which profitably uses a technique discussed in my latest S&SMC#23 and further includes a couple of very nice extras not present in the other versions. First though, I'd like to see what your own conjectures are (computing for N>10 might help to check them out) and of course your very own programs for your preferred hardware (RPL comes to mind, Prime's PPL, even SHARP or Casio models ...) or, if you're not up to the task (you know who you are), at least your very own personal comments or ideas ! (Googling for the solution is totally laaaaaame, boo !). Regards. V. . 

07112018, 11:10 AM
Post: #2




RE: [VA] SRC#001  Spiky Integral
Thanks for the post. What does SRC means in this case? Source problem?
Wikis are great, Contribute :) 

07122018, 04:23 AM
Post: #3




RE: [VA] SRC#001  Spiky Integral
(07102018 10:10 PM)Valentin Albillo Wrote: Using this numerical evidence I have my own conjecture on what the result will be for general N, which I'll post in a few days... Hello, Valentin, Here is my conjecture for the exact result when N = 39: \(\frac{756388295}{68719476736}\pi\) No programs. I’ve evaluated the integrals for N up to 12 on my CASIO fx991 LA X, which does it fast enough and gives exact results for N = 1, 2, 3, 4, 7 and 8. No googling, except for a OEIS sequence. Best regards, Gerson. 

07122018, 06:10 AM
Post: #4




RE: [VA] SRC#001  Spiky Integral
(07122018 04:23 AM)Gerson W. Barbosa Wrote: No googling, except for a OEIS sequence. I bet is one of these! Greetings, Massimo +×÷ ↔ left is right and right is wrong 

07122018, 10:32 AM
Post: #5




RE: [VA] SRC#001  Spiky Integral
With HP Prime
Code:
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 suggesting that I(N) will be nonzero if and only if N = 4k+3 and N = 4k+4, for k = 0, 1, 2.. 

07122018, 08:10 PM
Post: #6




RE: [VA] SRC#001  Spiky Integral
(07122018 10:32 AM)Pjwum Wrote: With HP Prime 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. 

07122018, 08:29 PM
(This post was last modified: 07122018 10:19 PM by Gerson W. Barbosa.)
Post: #7




RE: [VA] SRC#001  Spiky Integral
(07122018 06:10 AM)Massimo Gnerucci Wrote:(07122018 04:23 AM)Gerson W. Barbosa Wrote: No googling, except for a OEIS sequence. I thought this one was not at OEIS yet, but it is: 2, 10, 12, 17, 18, 19, 200, 201, 202... BTW, at first I searched for 140, 248, but these didn't match anything in the table. 70, 124 was successful. Edited to fix a typo Gerson. 

07122018, 10:43 PM
Post: #8




RE: [VA] SRC#001  Spiky Integral
.
Hi, all: Thanks for your interest in my SRC#001, some excellent contributions so far, much appreciated. Next Sunday night (GMT+1) I'll post my HP41C and HP71B solutions plus extras, but meanwhile I'll comment on some of your recent posts, read on: pier4r Wrote:Thanks for the post. What does SRC means in this case? Source problem? Nope, SRC = SemiRegular Column. Gerson W. Barbosa Wrote:Here is my conjecture for the exact result when N = 39: 756388295 / 68719476736*Pi Fully correct. You might want to check the exact result for N=71, namely: I(71) = 335205518724079925 / 73786976294838206464 * Pi which I got by actually computing the integral, then identifying the resulting constant, using dedicated programs I wrote myself (no OEIS involved). It's the hard way but hey, it's fun ! Pjwum Wrote:we can go beyond 10 [...] suggesting that I(N) will be nonzero if and only if N = 4k+3 and N = 4k+4, for k = 0, 1, 2.. Fully correct as well, we have a winner ! Well done, congratulations. 8) Gerson (again) Wrote:Sure we can: [...] It would take much, much longer or not, depending on how you go about computing them. Your results are correct, congratulations, but without explaining why would you compute the numerators and denominators that way it all seems a "magic trick" unrelated to the problem at hand as the relation to the integral seems shrouded in mystery ... 8D Regards. V. . 

07132018, 03:54 AM
Post: #9




RE: [VA] SRC#001  Spiky Integral  
07132018, 06:52 PM
(This post was last modified: 07142018 01:47 AM by ijabbott.)
Post: #10




RE: [VA] SRC#001  Spiky Integral
You don't really need to do any integration to solve this problem. You just need to use the trig identities to convert the product of cosines into a sum of cosines. When there is a nonzero constant term in the sum, the integral from \( 0 \) to \( 2\pi \) will be nonzero. The integrals of all the \( \cos(kx) \) terms from \( 0 \) to \( 2\pi \) will be zero.
The general form for converting the product of cosines to a sum of cosines is: \[ \cos(s)\cos(t) = \frac{1}{2}\cos(st) + \frac{1}{2}\cos(s+t) \] So: \[ \begin{eqnarray} \cos(1x) &=& \cos(1x) \end{eqnarray} \] \[ \begin{eqnarray} \cos(1x)\cos(2x) &=& \frac{1}{2}\cos(1x) + \frac{1}{2}\cos(3x) \end{eqnarray} \] \[ \begin{eqnarray} \cos(1x)\cos(2x)\cos(3x) &=& \frac{1}{2}\cos(1x)\cos(3x) + \frac{1}{2}\cos(3x)\cos(3x) \\ &=& \frac{1}{2}\left(\frac{1}{2}\cos(2x) + \frac{1}{2}\cos(4x)\right) + \frac{1}{2}\left(\frac{1}{2}\cos(0x) + \frac{1}{2}\cos(6x)\right) \\ &=& \frac{1}{4}\cos(2x) + \frac{1}{4}\cos(4x) + \frac{1}{4}\cos(0x) + \frac{1}{4}\cos(6x) \\ &=& \frac{1}{4}\cos(0x) + \frac{1}{4}\cos(2x) + \frac{1}{4}\cos(4x) + \frac{1}{4}\cos(6x) \end{eqnarray} \] \[ \begin{eqnarray} \cos(1x)\cos(2x)\cos(3x)\cos(4x) &=& \frac{1}{4}\cos(0x)\cos(4x) + \frac{1}{4}\cos(2x)\cos(4x) + \frac{1}{4}\cos(4x)\cos(4x) + \frac{1}{4}\cos(6x)\cos(4x) \\ &=& \frac{1}{4}\left(\frac{1}{2}\cos(4x) + \frac{1}{2}\cos(4x)\right) + \frac{1}{4}\left(\frac{1}{2}\cos(2x) + \frac{1}{2}\cos(6x)\right) + \frac{1}{4}\left(\frac{1}{2}\cos(0x) + \frac{1}{2}\cos(8x)\right) + \frac{1}{4}\left(\frac{1}{2}\cos(2x) + \frac{1}{2}\cos(10x)\right) \\ &=& \frac{1}{8}\cos(4x) + \frac{1}{8}\cos(4x) + \frac{1}{8}\cos(2x) + \frac{1}{8}\cos(6x) + \frac{1}{8}\cos(0x) + \frac{1}{8}\cos(8x) + \frac{1}{8}\cos(2x) + \frac{1}{8}\cos(10x) \\ &=& \frac{1}{8}\cos(0x) + \frac{1}{4}\cos(2x) + \frac{1}{4}\cos(4x) + \frac{1}{8}\cos(6x) + \frac{1}{8}\cos(8x) + \frac{1}{8}\cos(10x) \end{eqnarray} \] \[ \begin{eqnarray} \cos(1x)\cos(2x)\cos(3x)\cos(4x)\cos(5x) &=& \frac{1}{8}\cos(0x)\cos(5x) + \frac{1}{4}\cos(2x)\cos(5x) + \frac{1}{4}\cos(4x)\cos(5x) + \frac{1}{8}\cos(6x)\cos(5x) + \frac{1}{8}\cos(8x)\cos(5x) + \frac{1}{8}\cos(10x)\cos(5x) \\ &=& \frac{1}{8}\left(\frac{1}{2}\cos(5x) + \frac{1}{2}\cos(5x)\right) + \frac{1}{4}\left(\frac{1}{2}\cos(3x) + \frac{1}{2}\cos(7x)\right) + \frac{1}{4}\left(\frac{1}{2}\cos(1x) + \frac{1}{2}\cos(9x)\right) + \frac{1}{8}\left(\frac{1}{2}\cos(1x) + \frac{1}{2}\cos(11x)\right) + \frac{1}{8}\left(\frac{1}{2}\cos(3x) + \frac{1}{2}\cos(13x)\right) + \frac{1}{8}\left(\frac{1}{2}\cos(5x) + \frac{1}{2}\cos(15x)\right) \\ &=& \frac{1}{16}\cos(5x) + \frac{1}{16}\cos(5x) + \frac{1}{8}\cos(3x) + \frac{1}{8}\cos(7x) + \frac{1}{8}\cos(1x) + \frac{1}{8}\cos(9x) + \frac{1}{16}\cos(1x) + \frac{1}{16}\cos(11x) + \frac{1}{16}\cos(3x) + \frac{1}{16}\cos(13x) + \frac{1}{16}\cos(5x) + \frac{1}{16}\cos(15x) \\ &=& \frac{3}{16}\cos(1x) + \frac{3}{16}\cos(3x) + \frac{3}{16}\cos(5x) + \frac{1}{8}\cos(7x) + \frac{1}{8}\cos(9x) + \frac{1}{16}\cos(11x) + \frac{1}{16}\cos(13x) + \frac{1}{16}\cos(15x) \end{eqnarray} \] It's a bit tedious to continue, but you can see a constant term \( \frac{1}{4}\cos(0x) = \frac{1}{4} \) in the expansion of \( \prod_{i=1}^{3} \cos(ix) \), a constant term of \( \frac{1}{8}\cos(0x) = \frac{1}{8} \) in the expansion of \( \prod_{i=1}^{4} \cos(ix) \), and the constant term goes away in the expansion of \( \prod_{i=1}^{5} \cos(ix) \). The constant terms are the coefficients of \( \cos(0x) \) in those expansions whose terms consist of cosines of even multiples of \( x \). Your challenge, should you choose to accept it, is to come up with a general formula for each term (or a summation formula for all the terms). Of course, once you've found the coefficient for the constant \(\cos(0x)\) term, calculating the integral from \(0\) to \(2\pi\) is a piece of cake! (Just multiply the coefficient by \(2\pi\).) — Ian Abbott 

07142018, 01:30 AM
Post: #11




RE: [VA] SRC#001  Spiky Integral
I'm truly impressed at your diligence writing all those equations before essentially saying "... and so on". I'd probably be at least as impressed at the math, if I followed it, but I just wanted to give you kudos for all the work to explain your analysis.
Bob Prosperi 

07142018, 01:39 AM
Post: #12




RE: [VA] SRC#001  Spiky Integral
(07142018 01:30 AM)rprosperi Wrote: I'm truly impressed at your diligence writing all those equations before essentially saying "... and so on". I'd probably be at least as impressed at the math, if I followed it, but I just wanted to give you kudos for all the work to explain your analysis. Thanks! I'm slowly getting the hang of this MathJax lark! — Ian Abbott 

07142018, 02:08 PM
(This post was last modified: 07142018 02:09 PM by ijabbott.)
Post: #13




RE: [VA] SRC#001  Spiky Integral
Converting the product of cosines \(\prod_{i=1}^{n}\cos(ix) \) into a sum of cosines results in the angle multipliers being in the sum of cosines being all odd (when \( (n \bmod 4) \in \{1,2\} \)) or all even (when \( (n \bmod 4) \in \{0,3\} \)).
Identity: \[ \cos(s)\cos(t) = \frac{1}{2}\cos(st) + \frac{1}{2}\cos(s+t) \] Let \( s = \textrm{a}x, t = \textrm{b}x \). Then: \[ \cos(\textrm{a}x)\cos(\textrm{b}x) = \frac{1}{2}\cos((\textrm{a}\textrm{b})x) + \frac{1}{2}\cos((\textrm{a}+\textrm{b})x) \] When \(\textrm{a}\) and \(\textrm{b}\) are both odd or both even, then \((\textrm{a}\textrm{b})\) and \((\textrm{a}+\textrm{b})\) are both even, otherwise \((\textrm{a}\textrm{b})\) and \((\textrm{a}+\textrm{b})\) are both odd. This results in the factors \(\textrm{k}\) of \(x\) in the \(\cos(\textrm{k}x)\) terms of the summation switching between all odd and all even after every two \(\cos(ix)\) factors are appended to the product of cosines. As discussed in my earlier post, the summations with all odd \(\textrm{k}\), \(\cos(\textrm{k}x)\) terms all integrate to 0 over the interval \( [0,2\pi] \), but the summations with all even \(\textrm{k}\), \(\cos(\textrm{k}x)\) terms all include a constant term \(\textrm{q}\cos(0x)\) for some positive rational factor \(\textrm{q}\) which integrates to \(2\textrm{q}\pi\) over the interval \( [0,2\pi] \). — Ian Abbott 

07142018, 08:38 PM
Post: #14




RE: [VA] SRC#001  Spiky Integral
(07132018 06:52 PM)ijabbott Wrote: It's a bit tedious to continue We can use: \[cos(nx)=\frac{e^{inx}+e^{inx}}{2}\] Thus: \[ \begin{eqnarray} cos(mx)cos(nx) &=&(\frac{e^{imx}+e^{imx}}{2})(\frac{e^{inx}+e^{inx}}{2}) \\ &=& \frac{e^{i(m+n)x}+e^{i(m+n)x}}{4}+\frac{e^{i(mn)x}+e^{i(mn)x}}{4} \\ &=&\frac{cos(m+n)}{2}+\frac{cos(mn)}{2} \end{eqnarray} \] Let's forget the factor \(\frac{1}{2}\) for a moment and define: \[a_{k}=a^{k}+a^{k}\] For the same reason as above we have: \[ \begin{eqnarray} a_{m}a_{n}&=&(a^{m}+a^{m})(a^{n}+a^{n}) \\ &=&a^{m+n}+a^{(m+n)}+a^{mn}+a^{(mn)} \\ &=&a_{m+n}+a_{mn} \end{eqnarray} \] But of course \(a_{k}=a_{k}\). This allows us to rewrite it as: \[a_{m}a_{n}=a_{m+n}+a_{mn}\] We want to calculate the product: \[p_N=\Pi_{k=1}^{N}a_{k}\] Let's assume we already have \(p_{N1}\) represented as a sum of \(a_k\) with coefficients \(b_k\): \[p_{N1}=\Sigma_{k=1}^{M}b_{k}a_{k}\] Then \(p_{N}=a_{N}p_{N1}\) and thus: \[ \begin{eqnarray} p_{N} &=& a_{N}\Sigma_{k=1}^{M}b_{k}a_{k} \\ &=& \Sigma_{k=1}^{M}b_{k}a_{N}a_{k} \\ &=& \Sigma_{k=1}^{M}b_{k}(a_{N+k}+a_{Nk}) \end{eqnarray} \] This Python program allows us to calculate the coefficients \(b_k\): Code: def spiky(N): Here's the result for the first couple of values: Code: >>> for b in spiky(11): But we're only interested in the coefficient with the index 0: Code: >>> for b in spiky(40): Now it's time to deal with the factor \(\frac{1}{2}\) that we neglected. But that's trivial. We just have to add it to each factor \(a_{k}\). This leads to: \(\frac{1}{2^{N1}}\) And since we integrate the constant over \(2\pi\) we loose another 2. These are the examples for \(N=39\) and \(N=71\): Code: >>> b = spiky(100) (07122018 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}\] Thanks both for the challenge and the contributions. Cheers Thomas 

07142018, 10:54 PM
Post: #15




RE: [VA] SRC#001  Spiky Integral
.
Hi all, First of all, thanks for your interest and outstanding contributions, I've really enjoyed them all, much appreciated. As promised, these are my versions for other HP models (HP71B and HP41C) plus relevant comments: As already mentioned, an educated conjecture might be that I(N) is nonzero if the reminder of dividing N by 4 (that is, N mod 4) is either 0 or 3 (i.e.: N=3, 4, 7, 8, 11, 12, 15, 16,...) which actually is the correct answer. 1) My version for the HP71B is the following 220byte 6liner: 1 DESTROY ALL @ DIM S$[200] @ T=1/10^8 @ S$="COS(IVAR)" @ FOR N=1 TO 16 2 IF N>1 THEN S$=S$&"*COS("&STR$(N)&"*IVAR)" 3 S=VAL("INTEGRAL(0,2*PI,"&STR$(T)&","&S$&")") @ IF ABS(S)<T THEN 6 4 CALL IDENTIFY(S/PI,T$) @ T$=T$&"*Pi" 5 DISP USING "3D,':',2DZ.10D,' = ',20A";N,S,T$ 6 NEXT N which makes good use of a new technique I discovered long ago and then explained in my recent S&SMC#23, which in this case consist of incrementally creating an arbitrarily large string representing the function to be integrated, and then prepend to it the INTEGRAL keyword, the limits of integration and the tolerance before passing it to VAL to perform the actual computation. For instance, for N=16 this is the generated string passed to VAL: INTEGRAL(0,2*PI,.00000001,COS(IVAR)*COS(2*IVAR)*COS(3*IVAR)*COS(4*IVAR)*COS(5*IVAR)*COS(6*IVAR)*COS(7*IVAR)*COS(8*IVAR)*COS(9*IVAR) *COS(10*IVAR)*COS(11*IVAR)*COS(12*IVAR)*COS(13*IVAR)*COS(14*IVAR)*COS(15*IVAR)*COS(16*IVAR)) This means that, for each N, the string is parsed just the one time when passed to VAL but the function itself doesn't need any additional parsing no matter how many times the INTEGRAL keyword evaluates it, nor are there any loops whatsoever while evaluating it. This all results in extra simplicity (no loops) and much faster execution (neither loops nor calls to a multiline userdefined function or extra parsing). As a nice extra, line 4 (which is optional, can be omitted) does CALL my IDENTIFY subprogram (see my "Boldly Going  Identifying Constants" article) to identify every numeric result after integration and display the symbolic value. Finally, taking advantage of JF Garnier's Emu71's much greater speed, the results are provided not just up to N=10 but up to N=16, and results below the tolerance are considered to be 0 and aren't output. Let's see: >RUN 3: 1.5707963268 = 1/2*Pi 4: 0.7853981634 = 1/4*Pi 7: 0.3926990817 = 1/8*Pi 8: 0.3436116965 = 7/64*Pi 11: 0.2147573103 = 35/512*Pi 12: 0.1902136177 = 31/512*Pi 15: 0.1384417661 = 361/8192*Pi 16: 0.1259781722 = 657/16384*Pi Without the identification calls this runs in 19 sec. in my POPS system. 2) My version for the HP41C/Advantage ROM is a verbatim port of the one I gave for the HP15C, namely: 01 LBL "FX" 01 LBL "CINT" 02 STO 01 02 1.01 03 RCL 00 03 STO 00 04 INT 04 "FX" 05 STO 02 05 LBL 00 06 1 06 0 07 LBL 01 07 PI 08 RCL 02 08 ST+ X 09 RCL 01 09 INTEG 10 x 10 RCL 00 11 COS 11 STOP 12 x 12 X<>Y 13 DSE 02 13 STOP 14 GTO 01 14 ISG 00 15 END 15 GTO 00 16 END To run it (the results are exactly the same as the HP15C's but obtained much faster): RAD, FIX 3 XEQ "CINT" > 1.010 [R/S] > 1.083 04 [R/S] > 2.010 [R/S] > 5.974 05 [R/S] > 3.010 [R/S] > 1.571 [R/S] > 4.010 [R/S] > 0.785 [R/S] > 5.010 [R/S] > 3.778 08 [R/S] > 6.010 [R/S] > 2.733 09 [R/S] > 7.010 [R/S] > 0.393 [R/S] > 8.010 [R/S] > 0.344 [R/S] > 9.010 [R/S] > 3.600 04 [R/S] > 10.010 [R/S] > 1.664 09 With FIX 2, for N=4 we'd get I(4) = 0.86 + 0.03, which is 10% wrong. Details: N=4, f(x) = cos(x)*cos(2*x)*cos(3*x)*cos(4*x) FIX 3: 0.7853981871 + 0.003141592705, correct FIX 2: 0.8645215735 + 0.03141592659 , 10% wrong SCI 2: 0.7853981871 + 0.001484303379 . correct The exact value is Pi/4 = 0.785398163397 so we get almost 8 correct digits with FIX 3 or SCI 2 but hardly one with FIX 2. 3) Extra comments: A little theory now. After some algebraic massaging, any product of cosines of real arguments can be converted into a sum of exponentials of complex arguments using the following identity: Cos(x) = (e^(i*x) + e^(i*x))/2 This can be checked with this little HP71B routine: 10 DESTROY ALL @ COMPLEX A,B @ INPUT X 15 ! 20 DEF FNC1(X)=COS(X) 25 DEF FNC2(X)=COS(X)*COS(2*X) 30 DEF FNC3(X)=COS(X)*COS(2*X)*COS(3*X) 35 DEF FNC4(X)=COS(X)*COS(2*X)*COS(3*X)*COS(4*X) 40 ! 45 DEF FNE(X)=REPT(EXP((0,X))+EXP((0,X))) 50 ! 55 DEF FNE1(X)=1/2*FNE(X) 60 DEF FNE2(X)=1/4*(FNE(X)+FNE(3*X)) 65 DEF FNE3(X)=1/8*(FNE(0)+FNE(2*X)+FNE(4*X)+FNE(6*X)) 70 DEF FNE4(X)=1/16*(FNE(0)+2*FNE(2*X)+2*FNE(4*X)+FNE(6*X)+FNE(8*X)+FNE(10*X)) 75 ! 80 DISP 1;FNC1(X),FNE1(X) @ DISP 2;FNC2(X),FNE2(X) @ DISP 3;FNC3(X),FNE3(X) @ DISP 4;FNC4(X),FNE4(X) where the FNC1, FNC2, FNC3, FNC4 are userdefined functions (UDFs) implementing the products of 1, 2, 3 and 4 cosines, while the FNE1, FNE2, FNE3, FNE4 are the equivalent sums of 1, 2, 4 and 8 complex exponentials (not all necessarily distinct). Let's run it: >RUN ? .2018 1 .979707385531 .97970738553 2 .900990956269 .90099095627 3 .740861912405 .740861912405 4 .512323594236 .512323594236 As you may see, every FNE1..4 gives the same result as the corresponding FNC1..4 for an arbitrary argument (.2018). What is to be gained by converting a product of cosines into a sum of exponentials ? Two main things:
Again, thanks for your interest and regards. V. . 

07152018, 02:53 PM
(This post was last modified: 07152018 02:55 PM by Gerson W. Barbosa.)
Post: #16




RE: [VA] SRC#001  Spiky Integral
(07142018 08:38 PM)Thomas Klemm Wrote: These are the examples for \(N=39\) and \(N=71\): In exact mode, the HP 50g does it automatically. (07142018 08:38 PM)Thomas Klemm Wrote:(07122018 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\). 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 threedigit results for n = 3 and n = 4. Then I evaluated the integrals for n up to 12 on my CASIO fx991 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^16+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 timeconsuming. 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*(16/(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 HP42S program computes the integral for n >= 1 and n up to about 20000 (on Free42) and returns exact results for n < 15. Code:
Gerson. 

07152018, 06:26 PM
Post: #17




RE: [VA] SRC#001  Spiky Integral
This RPL program calculates the coefficients:
Code: « Example: The value for 4 is: { 1 0 2 0 2 0 1 0 1 0 1 } To get the value for 5 we create the following lists: { 0 0 0 0 0 1 0 2 0 2 0 1 0 1 0 1 } { 0 2 0 2 0 1 0 0 0 0 0 0 0 0 0 0 } { 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 } And then we just ADD them up: { 0 3 0 3 0 3 0 2 0 2 0 1 0 1 0 1 } The 2nd and the 3rd list is just the 1st list reversed and then again mirrored at the left border. That's a consequence of \(a_k=a^k+a^{k}\) being symmetric, that is \(a_k=a_{k}\). We don't want negative indices. Cheers Thomas PS: Is there a better way to create a list of m zeros? 

07152018, 07:10 PM
(This post was last modified: 07152018 07:16 PM by Gerson W. Barbosa.)
Post: #18




RE: [VA] SRC#001  Spiky Integral
The RPL program has been rewritten so as not to needlessly waste time expanding a polynomial when results are supposed to be zero. %%HP: T(3)A(R)F(.); \<< DUPDUP DUPDUP * + 2 / 2 MOD IF NOT THEN DUP '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 OVER "+" POS 1 + SWAP SUB OBJ\> 2 ROT 1  ^ / \pi * ELSE DROP2 0 END \>> 71 > '335205518724079925/73786976294838206464*π' > NUM > 1.42718843886E2 267.5 bytes, but it takes too long for this example (28m 34s... on the emulator!) As a comparison, the RPN program on Free42 returns 1.4271(9054304)E2 instantly (maybe a couple of seconds on a real 42S). 

07152018, 07:35 PM
(This post was last modified: 07152018 07:36 PM by Gerson W. Barbosa.)
Post: #19




RE: [VA] SRC#001  Spiky Integral  
07152018, 07:53 PM
Post: #20




RE: [VA] SRC#001  Spiky Integral
(07152018 07:35 PM)Gerson W. Barbosa Wrote:(07152018 06:26 PM)Thomas Klemm Wrote: PS: Is there a better way to create a list of m zeros? Gerson's method above is also how I would do it if the ListExt library isn't available. Using the ListExt library: « 0 m LMRPT » Creating a list of 1000 0s with the CON method: 0.77s Creating the same list with the ListExt command: 0.07s Either should give reasonable performance in this scenario. 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)