Threaded Mode | Linear Mode

Valentin Albillo · 07-14-2018, 10:54 PM

.
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 (HP-71B and HP-41C) plus relevant comments:

As already mentioned, an educated conjecture might be that I(N) is non-zero 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 HP-71B is the following 220-byte 6-liner:

      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 multi-line user-defined 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 J-F 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 HP-41C/Advantage ROM is a verbatim port of the one I gave for the HP-15C, 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 HP-15C'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 HP-71B 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 user-defined 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:

the integral of a sum of functions is the sum of the integrals of each function separately
the integral of a simple exponential function is another simple exponential function so no need for numerical integration

In this way, the integral of the product of cosines becomes a sum of exponential functions and no numerical approximations are required. There's of course the drudgery of determining the components of the sum but that's another story ! XD

Again, thanks for your interest and regards.
V.
.