The Museum of HP Calculators


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

This program is Copyright © 2003 by Jean-Marc Baillard and is used here by permission.

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 = uk + uk+1 +  uk+2 + uk+3 + .........
            b)  "SIGMAR" -------------- same series sum when  un 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 = uk + uk+1 +  uk+2 + uk+3 + .........    =  Sigma  un    for  n >= k
  defined by    un = 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  un = 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/nn + .........

-Here we have  f(n) = 1/nn  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  un  is defined by a recurrence relation
 

-The following program calculates  S = uk + uk+1 +  uk+2 + uk+3 + .........    =  Sigma  un    for  n >= k
  when  uk is given  and   un+1 =  f(un;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 = un
Flags: /
Subroutine:  a program which computes  un+1 = f(un;n)  assuming un 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             uk             S

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

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

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

"T"  ASTO 00
 1    ENTER^   ( k = 1 = u1 )
       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 = u0 - u1 +  u2 - u3 + ...... + (-1)n  un + .......
-The sum is rewritten:   S = a0/2 + ( C11 a0 - C10 a1 ) / 22 + ( C22 a0 - C21 a1 + C20 a2  ) / 23 +  ( C33 a0 - C32 a1 + C31 a2 - C30 a3  ) / 24 + ................
   where   Cnp =  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  un = 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 (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 1018 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  u0 - u1 +  u2 - u3 + ...... -  u9 =  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  un;m    for  n >= n0 ; m >= m0
 

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

Flags: /
Subroutine:  a program which computes  un;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             m0             S
           X             n0             S

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

-Here, we have   f(n;m) = 1 / ( nn 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  un;m;p    for  n >= n0 ; m >= m0 ; p >= p0
 

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

Flags: /
Subroutine:  a program which computes  un;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             p0             /
           Y             m0             S
           X             n0             S

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

-Here, we have   f(n;m;p) = 1 / ( nn m! (p!)2 )  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(n1;n2; ..... ;nk)   with  n1>= 0 ;  n2>= 0 ; ........ ;  nk>= 0
 

Data Registers:           •  R00 =  f name  ( this register is to be initialized )
                                          R01 = n1   R02 = n2  .........  Rkk = nk
 

Flags: /
Subroutine:  a program which computes  u(n1;n2; ..... ;nk)  = f( n1;n2; ..... ;nk)  assuming n1is in R01 ; n2 is in R02 ; ..... ; nk 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 / ( nn m! (p!)2 )      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