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