hp41programs

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 En(x)
 

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 )
 
 

      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!/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: /

-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          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