The Museum of HP Calculators

# Exponential , Sine , Cosine Integrals for the HP-41

This program is Copyright © 2006-2007 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

1°) Series Expansions for Ei(x) , Si(x) , Ci(x) , Shi(x) , Chi(x)
2°) Asymptotic Series for Si(x) & Ci(x)
3°) A Continued Fraction for Ci(x) & Si(x)  ( x > 3 )   ( new )
4°) Exponential Integrals En(x) ( improved )

1°) Series Expansions

-These functions are defined by:

Ei(x) = §-infinityx  et/t dt      ;      Si(x) = §0x sin(t)/t dt      ;      Ci(x) = C + Ln(x) + §0x (cos(t)-1)/t dt

Shi(x) = §0x sinh(t)/t dt      ;      Chi(x) = C + Ln(x) + §0x (cosh(t)-1)/t dt

-The program computes  Ei(x) , Si(x) , Ci(x) , Shi(x) and Chi(x) with the following formulae:

Ei(x)  = C + ln(x) + Sumn=1,2,.....  xn/(n.n!)
Si(x)  = Sumn=0,1,2,..... (-1)n x2n+1/((2n+1).(2n+1)!)
Ci(x)  = C + ln(x) + Sumn=1,2,..... (-1)n x2n/(2n.(2n)!)                 where  C = 0.5772156649 is Euler's constant.
Shi(x) = Sumn=0,1,2,.....  x2n+1/((2n+1).(2n+1)!)
Chi(x)= C + ln(x) + Sumn=1,2,.....  x2n/(2n.(2n)!)

Data Registers:  R00 = x , R01 to R04: temp
Flags: /
Subroutines: /

01  LBL "EI"
02  STO 00
03  STO 04
04  1
05  STO 01
06  STO 03
07  X<>Y
08  GTO 03
09  LBL "SI"
10  STO 00
11  STO 02
12  X^2
13  CHS
14  GTO 01
15  LBL "CI"
16  STO 00
17  X^2
18  CHS
19  GTO 02
20  LBL "SHI"
21  STO 00
22  STO 02
23  X^2
24  LBL 01
25  STO 04
26  1
27  STO 01
28  2
29  STO 03
30  LASTX
31  GTO 04
32  LBL "CHI"
33  STO 00
34  X^2
35  LBL 02
36  STO 04
37  2
38  STO 01
39  STO 03
40  /
41  LBL 03
42  STO 02
43  RCL 01
44  /
45  XEQ 04
46  RCL 00
47  LN
48  .5772156649
49  +
50  +
51  RTN
52  LBL 04
53  RCL 02
54  RCL 04
55  *
56  RCL 01
57  RCL 03
58  ST+ 01
59  SIGN
60  +
61  RCL 01
62  X=Y?
63  SIGN
64  *
65  /
66  STO 02
67  RCL 01
68  /
69  +
70  X#Y?
71  GTO 04
72  END

( 119 bytes / SIZE 005 )

Example:    1.4   XEQ "EI"     >>>>  Ei(1.4) =  3.007207463    ( in 8 seconds )
1.4   XEQ "SI"     >>>>  Si(1.4) =  1.256226733    ( in 4 seconds )
1.4   XEQ "CI"    >>>>  Ci(1.4) =  0.462006585    ( in 5 seconds )
1.4   XEQ "SHI"  >>>> Shi(1.4) =  1.561713387   ( in 4 seconds )
1.4   XEQ "CHI" >>>> Chi(1.4) =  1.445494076   ( in 5 seconds )

Note:   Logarithmic Integral = Li(x) = Ei(Ln(x))    ( x > 1 )    e.g.   Li(100) = 30.12614160
( simply add  LBL "LI"   LN   before the first line if needed )

2°) Asymptotic Series  for  Si(x) & Ci(x)

-For large arguments, Si(x) and Ci(x) are not calculated with a good accuracy by the above program.
-The following one uses the formulae:

Si(x) = pi/2 - f(x) cos x - g(x) sin x                  where    f(x) ~  ( 1 - 2!/x2 + 4!/x4 - 6!/x6 + ..... )/x
Ci(x) = f(x) sin x - g(x) cos x                             and     g(x) ~  ( 1 - 3!/x2 + 5!/x4 - 7!/x6 + ..... )/x2

-Actually, these series are divergent but the error is small if we stop at a proper instant.

Data Registers:  R00 = x , R01 to R03: temp
Flags: /
Subroutines: /

01  LBL "SCI"
02  STO 00
03  CLX
04  STO 01
05  SIGN
06  STO 02
07  XEQ 01
08  RCL 00
09  /
10  STO 03
11  1
12  STO 01
13  STO 02
14  XEQ 01
15  RCL 00
16  STO Z
17  X^2
18  /
20  P-R
21  RCL 00
22  RCL 03
23  P-R
24  ST+ T
25  X<> Z
26  -
27  1
28  ASIN
29  R^
30  -
31  DEG
32  RTN
33  LBL 01
34  RCL 02
35  RCL 01
36  2
37  +
38  STO 01
39  ENTER^
40  DSE X
41  *
42  *
43  RCL 00
44  X^2
45  CHS
46  /
47  STO 02
48  +
49  X#Y?
50  GTO 01
51  END

( 68 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Y / Ci(x) X x Si(x)

Example:    30   XEQ "SCI"   >>>>   Si(30) =   1.566756540     ( in 14 seconds )
X<>Y               Ci(30) =  -0.033032417

Notes:

-This program will not work if x is too small:
-Divergence already appears too soon for x = 29 but if you add  RND  X<>Y  RND  X<>Y  after line 48
and run the program in FIX 6 ( for instance ) you will have a good approximation:

-With this modification ( and FIX 6 )    20  XEQ "SCI"  yields   Si(20) = 1.54824168  and  Ci(20) = 0.044419864  ( errors are smaller than 10-7 )

3°) A Continued Fraction for  Ci(x) & Si(x)   ( x > 3 )

Formula:      -Ci(x) + i.Si(x) - i.pi/2 = exp(-i.x) ( 1/(1+i.x-12/(3+i.x-22/(5+i.x-32/(7+i.x- ..... )))) )

-The program below uses this complex continued fraction to compute  Ci(x) & Si(x)    for x greater or equal to 3
( use the program at the top for smaller x-values )
-The number of loops is fixed to 22 ( line 03 ) to produce accurate results for x >= 3

Data Registers:   R00 & R01: temp
Flags: /
Subroutines: /

01  LBL "CISI"
02  STO 01
03  22
04  STO 00
05  CLST
07  LBL 01
08  RCL 00
09  ST+ X
10  1
11  +
12  +
13  X<>Y
14  RCL 01
15  +
16  X<>Y
17  R-P
18  X<>Y
19  CHS
20  RCL 00
21  X^2
22  CHS
23  RCL Z
24  /
25  P-R
26  DSE 00
27  GTO 01
28  1
29  +
30  X<>Y
31  RCL 01
32  +
33  X<>Y
34  R-P
35  1/X
36  X<>Y
37  RCL 01
38  +
39  CHS
40  X<>Y
41  P-R
42  CHS
43  1
44  ASIN
45  RCL Z
46  +
47  X<>Y
48  DEG
49  END

( 64 bytes / SIZE 002 )

 STACK INPUTS OUTPUTS Z / Si(x)-pi/2 Y / Si(x) X x >= 3 Ci(x)

Execution time = 34 seconds

Examples:

3  XEQ "CISI"  >>>>   Ci(3)  = 0.119629786
RDN         Si(3)  = 1.848652528
RDN  Si(3)-pi/2  = 0.277856201

100  R/S    >>>>     Ci(100)  = -0.005148824724
RDN         Si(100)  =   1.562225467
RDN  Si(100)-pi/2  = -0.008570860160

Note:

-In fact, the series expansions still produce 9 exact decimals for x = 6
-So, if you want to use the continued fraction for x > 6 only, the number of required loops may be decreased:
-Simply replace line 03 by 14 and the execution time is reduced to 23 seconds.

4°) Exponential Integrals:  En(x)

-We have  En(x) = §1+infinity  t -n e -x.t dt        ( x > 0 ;  n = a positive integer )
-This function is computed by a series expansion if x <= 1.5

En(x) = (-x)n-1 ( -Ln x - gamma + Sumk=1,...,n-1 1/k ) / (n-1)!  - Sumk#n-1 (-x)k / (k-n+1) / k!          where  gamma = Euler's Constant = 0.5772156649...

and by the following continued fraction if  x > 1.5  ( line 26 )

En(x) = exp (-x) ( 1/(x+n-1.n/(x+n+2-2(n+1)/(x+n+4- .... ))) )

-We also have:   E0(x) = (1/x).exp(-x)  and   En(0) = 1/(n-1)  if  n > 1

Data Registers:     R00 thru R04: temp
Flags: /
Subroutines: /

01  LBL "ENX"
02  X#0?
03  GTO 00
04  X=Y?                         lines 04-05 are only useful to produce a DATA ERROR if  x = n = 0
05  LN
06  CLX
07  SIGN
08  -
09  1/X
10  RTN
11  LBL 00
12  STO 01
13  X<>Y
14  X#0?
15  GTO 00
16  X<>Y
17  CHS
18  E^X
19  LASTX
20  CHS
21  /
22  RTN
23  LBL 00
24  STO 02
25  CLX
26  1.5
27  X<Y?
28  GTO 05
29  RCL 02
30  STO 04
31  1
32  STO 00
33  STO 03
34  -
35  ENTER^
36  X=0?
37  GTO 02
38  1/X
39  GTO 04
40  LBL 01
41  *
42  LBL 02
43  DSE 04
44  X<0?
45  GTO 03
46  RCL 04
47  1/X
48  +
49  GTO 02
50  LBL 03
51  RCL 01
52  LN
53  .5772156649
54  +
55  -
56  RCL 01
57  CHS
58  RCL 02
59  DSE X
60  ""                         ( TEXT 0 or another NOP instruction like STO X ... )
61  Y^X
62  LASTX
63  FACT
64  /
65  *
66  +
67  LBL 04
68  RCL 00
69  RCL 01
70  CHS
71  *
72  RCL 03
73  /
74  STO 00
75  ISG 03
76  CLX
77  RCL 03
78  RCL 02
79  -
80  X=0?
81  GTO 01
82  /
83  -
84  X#Y?
85  GTO 04
86  RTN
87  LBL 05
88  26                               the number of loops for the continued fraction
89  STO 00
90  CLX
91  LBL 06
92  RCL 00
93  ST+ X
94  RCL 01
95  +
96  RCL 02
97  +
98  X<>Y
99  -
100  RCL 02
101  1
102  -
103  RCL 00
104  ST+ Y
105  *
106  X<>Y
107  /
108  DSE 00
109  GTO 06
110  RCL 01
111  RCL 02
112  +
113  X<>Y
114  -
115  1/X
116  RCL 01
117  CHS
118  E^X
119  *
120  END

( 157 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS Y n / X x En(x)

where  x >= 0  &  n is a non-negative integer

Examples:

0    ENTER^
1.4  XEQ "ENX"  >>>>  E0(1.4) = 0.176140689        ( in 0.6 second )

2    ENTER^
1.4  XEQ "ENX"  >>>>  E2(1.4) = 0.0838899263      ( in 11 seconds )

100  ENTER^
1.4  XEQ "ENX"  >>>>  E100(1.4) = 0.0024558006    ( in 9 seconds )

3    ENTER^
2    XEQ "ENX"  >>>>  E3(2) = 0.03013337978       ( in 16 seconds )

100  ENTER^
XEQ "ENX"  >>>>  E100(100) = 1.864676429 E-46    ( in 16 seconds )

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4