ExpInt

Exponential , Sine , Cosine Integrals for the HP-41

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 )
4°) Exponential Integrals E_n(x)

1°) Series Expansions

-These functions are defined by:

Ei(x) = §_-infinity^x e^t/t dt ; Si(x) = §₀^x sin(t)/t dt ; Ci(x) = C + Ln(x) + §₀^x (cos(t)-1)/t dt

Shi(x) = §₀^x sinh(t)/t dt ; Chi(x) = C + Ln(x) + §₀^x (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) + Sum_n=1,2,..... xⁿ/(n.n!)
   Si(x) = Sum_{n=0,1,2,.....} (-1)ⁿ x²ⁿ⁺¹/((2n+1).(2n+1)!)
Ci(x) = C + ln(x) + Sum_n=1,2,..... (-1)ⁿ x²ⁿ/(2n.(2n)!)                 where C = 0.5772156649 is Euler's constant.
Shi(x) = Sum_{n=0,1,2,.....} x²ⁿ⁺¹/((2n+1).(2n+1)!)
Chi(x)= C + ln(x) + Sum_n=1,2,..... x²ⁿ/(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 )

STACK	INPUTS	OUTPUTS
X	x	f(x)

Examples:

   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!/x² + 4!/x⁴ - 6!/x⁶ + ..... )/x
Ci(x) = f(x) sin x - g(x) cos x and g(x) ~ ( 1 - 3!/x² + 5!/x⁴ - 7!/x⁶ + ..... )/x²

-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 /
19 RAD
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-1²/(3+i.x-2²/(5+i.x-3²/(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
06 RAD
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: E_n(x)

-We have E_n(x) = §₁^+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

E_n(x) = (-x)^n-1 ( -Ln x - gamma + Sum_k=1,...,n-1 1/k ) / (n-1)! - Sum_k#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 )

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

-We also have: E₀(x) = (1/x).exp(-x) and E_n(0) = 1/(n-1) if n > 1

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

-Lines 04-05 are only useful to produce a DATA ERROR if x = n = 0

01 LBL "ENX"
02 X#0?
03 GTO 00
04 X=Y?
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 ""
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
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	E_n(x)

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

Examples:

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

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

100 ENTER^
1.4 XEQ "ENX" >>>> E₁₀₀(1.4) = 0.0024558006 ( in 9 seconds )

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

100 ENTER^
XEQ "ENX" >>>> E₁₀₀(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

hp41programs

Exponential , Sine , Cosine Integrals for the HP-41