The Museum of HP Calculators
The Museum of HP Calculators
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 /
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-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
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: 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
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall