The Museum of HP Calculators

Series ( including the Euler acceleration method ) for the HP-41

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

-Six programs are listed hereafter:

      1°) a) "SIGMA1" computes the "infinite" sum S = u_k + u_k+1 + u_k+2 + u_k+3 + .........
            b) "SIGMAR" -------------- same series sum when u_n is defined by a recurrence relation
            c) "EULER" uses the Euler transformation for alternating series

      2°) "SIGMA2" calculates double series
      3°) "SIGMA3" --------- triple series
      4°) "SIGMA"   deals with the general case of multiple series.

Warning: -In these programs, the iterations stop when 2 consecutive partial sums are equal.
Therefore, wrong results may be obtained if one term is negligible and the next one is not!

1°) a) Simple Series

-The following program calculates S = u_k + u_k+1 + u_k+2 + u_k+3 + ......... = Sigma u_n for n >= k
defined by u_n = f(n) where f is a known function.

Data Registers:           • R00 = f name ( this register is to be initialized before executing "SIGMA1" )
                                          R01 = S
                                          R02 = n

Flags: /
Subroutine: a program which computes u_n = f(n) assuming n is in X-register upon entry

01 LBL "SIGMA1"
02 STO 02
03 CLX
04 STO 01
05 LBL 01
06 RCL 02
07 XEQ IND 00
08 ISG 02
09 LBL 02
10 RCL 01
11 +
12 STO 01
13 LASTX
14 X#Y?
15 GTO 01
16 END

( 30 bytes / SIZE 003 )

      STACK        INPUTS      OUTPUTS

           Y             /             S

           X             k             S

Example: Calculate S = 1 + 1/4 + 1/27 + ....... + 1/nⁿ + .........

-Here we have f(n) = 1/nⁿ and k = 1 and this function may be obtained by

LBL "T"
ENTER^
CHS
Y^X
RTN

-Then, "T" ASTO 00
1 XEQ "SIGMA1" yields X = 1.291285997 = R01 ( in 9 seconds )

1°) b) Simple Series: when u_n is defined by a recurrence relation

-The following program calculates S = u_k + u_k+1 + u_k+2 + u_k+3 + ......... = Sigma u_n for n >= k
when u_k is given and u_n+1 = f(u_n;n) where f is a known function.

Data Registers:           • R00 = f name ( this register is to be initialized before executing "SIGMAR" )
                                          R01 = S
                                          R02 = n
                                          R03 = u_n
Flags: /
Subroutine: a program which computes u_n+1 = f(u_n;n) assuming u_n is in X-register and n is in Y-register upon entry

01 LBL "SIGMAR"
02 STO 01
03 STO 03
04 X<>Y
05 STO 02
06 LBL 01
07 RCL 02
08 RCL 03
09 XEQ IND 00
10 STO 03
11 ISG 02
12 LBL 02
13 RCL 01
14 +
15 STO 01
16 LASTX
17 X#Y?
18 GTO 01
19 END

( 33 bytes / SIZE 004 )

      STACK        INPUTS      OUTPUTS

           Y             k             S

           X             u_k             S

Example: Calculate S = 1 + 1/2! + 1/3! + ....... + 1/n! + .........

-It may be obtained using the first program but we can also define u₁ = 1 and u_n+1 = u_n/(n+1)

LBL "T"
X<>Y
1
+
/
RTN

"T" ASTO 00
1 ENTER^ ( k = 1 = u₁ )
XEQ "SIGMAR" produces X = 1.718281830 = R01 in 7 seconds ( the exact result is e-1 )

1°) c) The Euler Transformation

-The Euler transformation is an acceleration method developed for alternating series:    S = u₀ - u₁ + u₂ - u₃ + ...... + (-1)ⁿ u_n + .......
-The sum is rewritten:   S = a₀/2 + ( C₁¹ a₀ - C₁⁰ a₁ ) / 2² + ( C₂² a₀ - C₂¹ a₁ + C₂⁰ a₂ ) / 2³ + ( C₃³ a₀ - C₃² a₁ + C₃¹ a₂ - C₃⁰ a₃ ) / 2⁴ + ................
   where   C_n^p = n! / ( p! ( n-p )! )   are the binomial coefficients.

-This may produce superb acceleration but it can also fail...

Data Registers:           • R00 = f name ( this register is to be initialized before executing "EULER" )
                                          R01 = S
                                          R02 thru R05: temp

Flags: /
Subroutine: a program which computes u_n = f(n) assuming n is in X-register upon entry

01 LBL "EULER"
02 CLX
03 XEQ IND 00
04 2
05 /
06 STO 01
07 1
08 STO 02
09 LBL 01
10 RCL 02
11 STO 03
12 SIGN
13 STO 04
14 CLX
15 STO 05
16 LBL 02
17 RCL 02
18 RCL 03
19 -
20 XEQ IND 00
21 RCL 04
22 *
23 ST+ 05
24 LASTX
25 CHS
26 RCL 03
27 *
28 RCL 02
29 LASTX
30 -
31 1
32 +
33 /
34 STO 04
35 DSE 03
36 GTO 02
37 RCL 02
38 XEQ IND 00
39 RCL 04
40 *
41 RCL 05
42 +
43 2
44 RCL 02
45 1
46 +
47 STO 02
48 Y^X
49 /
50 RCL 01
51 +
52 STO 01
53 LASTX
54 VIEW X
55 X#Y?
56 GTO 01
57 END

( 75 bytes / SIZE 006 )

Example: Compute S = 1 - 2^-1/2 + 3^-1/2 - 4^-1/2 + ....... + (-1)ⁿ (n+1)^-1/2 + ......

-Here, we have f(n) = (n+1)^-1/2

LBL "T"
SIGN
LASTX
+
SQRT
1/X
RTN

"T" ASTO 00 XEQ "EULER" the successive approximations are displayed and 6mn27s later, X = 0.6048986431 = R01

-The error is -3 10^-10 (!)
-Actually, only 28 terms are calculated
-Without an acceleration method, more than 10¹⁸ terms would be necessary to achieve the same accuracy...
( execution time would be much greater than the age of the Universe! )

Notes: a) This program computes the different f(n) several times and this may be time-consuming if f is a complicated function.
An alternative is to store the successive terms into data registers:

-Add   ST+ 06          after line 45
-Add   STO IND 06 after line 38
-Replace line 20 by   7 + RCL IND X
-Add   8 STO 06     after line 08
-Add   STO 07         after line 03

-The same result is now obtained in 5mn19s ( but the SIZE must be at least 035 in this example )

   b) It is often better to sum the first terms and apply the Euler transformation to the rest:
       -With the example above, if you sum u₀ - u₁ + u₂ - u₃ + ...... - u₉ = 1 - 2^-1/2 + 3^-1/2 - 4^-1/2 + ....... - 10^-1/2
         and execute "EULER" with

LBL "T"
11
+
SQRT
1/X
RTN the result is obtained in 1mn22s

2°) Double Series

-The following program calculates S = Sigma u_n;m for n >= n₀ ; m >= m₀

Data Registers:           • R00 = f name ( this register is to be initialized )
                                          R01 = n₀    R03 = S    R05 = m
                                          R02 = m₀   R04 = n

Flags: /
Subroutine: a program which computes u_n;m = f(n;m) assuming n is in X-register and m is in Y-register upon entry

01 LBL "SIGMA2"
02 STO 01
03 STO 04
04 X<>Y
05 STO 02
06 STO 05
07 CLX
08 STO 03
09 DSE 05
10 LBL 01
11 ISG 05
12 CLX
13 RCL 05
14 RCL 04
15 XEQ IND 00
16 RCL 03
17 +
18 STO 03
19 LASTX
20 X#Y?
21 GTO 01
22 ISG 04
23 CLX
24 RCL 02
25 STO 05
26 RCL 04
27 XEQ IND 00
28 RCL 03
29 +
30 STO 03
31 LASTX
32 X#Y?
33 GTO 01
34 END

( 52 bytes / SIZE 006 )

      STACK        INPUTS      OUTPUTS

           Y             m₀             S

           X             n₀             S

Example: Calculate S = Sigma 1 / ( nⁿ m! ) with n >= 1 and m >= 1

-Here, we have f(n;m) = 1 / ( nⁿ m! ) and this function may be obtained by

LBL "T"
ENTER^
Y^X
X<>Y
FACT
*
1/X
RTN

-Then "T" ASTO 00
1 ENTER^
XEQ "SIGMA2" and 1mn18s later X = 2.218793264 = R03

Note:     If the series converges too slowly, execution time becomes prohibitive.
           Lines 16 to 19 and 28 to 31 may then be replaced by    X<> 03   ST+ 03 RND RCL 03 RND
       and the accuracy will be determined by the display format.
-Similar modifications apply to the other programs listed in this page.

3°) Triple Series

-The following program calculates S = Sigma u_n;m;p for n >= n₀ ; m >= m₀ ; p >= p₀

Data Registers:           • R00 = f name ( this register is to be initialized )
                                          R01 = n₀    R04 = S    R05 = n
                                          R02 = m₀                    R06 = m
                                          R03 = p₀                     R07 = p

Flags: /
Subroutine: a program which computes u_n;m;p = f(n;m;p) assuming n is in X-register ; m is in Y-register and p is in Z-register upon entry.

01 LBL "SIGMA3"
02 STO 01
03 STO 05
04 RDN
05 STO 02
06 STO 06
07 X<>Y
08 STO 03
09 STO 07
10 CLX
11 STO 04
12 DSE 07
13 LBL 01
14 ISG 07
15 CLX
16 RCL 07
17 RCL 06
18 RCL 05
19 XEQ IND 00
20 RCL 04
21 +
22 STO 04
23 LASTX
24 X#Y?
25 GTO 01
26 ISG 06
27 CLX
28 RCL 03
29 STO 07
30 RCL 06
31 RCL 05
32 XEQ IND 00
33 RCL 04
34 +
35 STO 04
36 LASTX
37 X#Y?
38 GTO 01
39 ISG 05
40 CLX
41 RCL 03
42 STO 07
43 RCL 02
44 STO 06
45 RCL 05
46 XEQ IND 00
47 RCL 04
48 +
49 STO 04
50 LASTX
51 X#Y?
52 GTO 01
53 END

( 74 bytes / SIZE 008 )

      STACK        INPUTS      OUTPUTS

           Z             p₀             /

           Y             m₀             S

           X             n₀             S

Example: Calculate S = Sigma 1 / ( nⁿ m! (p!)² ) with n >= 1 ; m >= 1 ; p >= 1

-Here, we have f(n;m;p) = 1 / ( nⁿ m! (p!)² ) and this function may be obtained by

LBL "T"
ENTER^
Y^X
X<>Y
FACT
*
X<>Y
FACT
X^2
*
1/X
RTN

-Then "T" ASTO 00
1 ENTER^ ENTER^
XEQ "SIGMA3" and 6mn53s later, X = 2.839135243 = R04 ( error = -7 10^-9 )

4°) Multiple Series

-The following program calculates S = Sigma u(n₁;n₂; ..... ;n_k) with n₁>= 0 ; n₂>= 0 ; ........ ; n_k>= 0

Data Registers: • R00 = f name ( this register is to be initialized )
R01 = n₁ R02 = n₂ ......... Rkk = n_k

Flags: /
Subroutine: a program which computes u(n₁;n₂; ..... ;n_k) = f( n₁;n₂; ..... ;n_k) assuming n₁is in R01 ; n₂ is in R02 ; ..... ; n_k is in Rkk upon entry.

Note: Synthetic registers M N O may be replaced by any unused data registers.

01 LBL "SIGMA"
02 CLA
03 STO N
04 STO O
05 E3
06 /
07 ISG X
08 CLRGX if you don't have an HP-41CX , replace this line with 0 LBL 00 STO IND Y ISG Y GTO 00
09 DSE IND N
10 LBL 01
11 ISG IND N
12 CLX
13 XEQ IND 00
14 RCL M
15 +
16 STO M
17 LASTX
18 X#Y?
19 GTO 01
20 LBL 02
21 CLX
22 STO IND N
23 DSE N
24 X<0?
25 GTO 03
26 ISG IND N
27 CLX
28 XEQ IND 00
29 RCL M
30 +
31 STO M
32 LASTX
33 X=Y?
34 GTO 02
35 RCL O
36 STO N
37 GTO 01
38 LBL 03
39 X<>Y
40 CLA
41 END

( 73 bytes / SIZE kkk+1 )

STACK INPUTS OUTPUTS

X k S

k = 1 for simple series ; 2 for double series ; 3 for triple series ; 4 for quadruple series ... etc ...

Example: Let's take once again the triple series S = Sigma 1 / ( nⁿ m! (p!)² ) for n >= 1 ; m >= 1 ; p >= 1

-We make a change of arguments to reduce to the standard: n >= 0 ; m >= 0 ; p >= 0 by replacing n with (n+1) ; m with (m+1) ; p with (p+1)

LBL "T"
RCL 01
1
+
ENTER^
Y^X
RCL 02
1
+
FACT
*
RCL 03
1
+
FACT
X^2
*
1/X
RTN

-Then 3 XEQ "SIGMA" yields X = 2.839135243 ( in 8mn39s )

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall