Fresnel Integrals for the HP-41

Overview
 

 1°) Ascending Series
 2°) Asymptotic Series
 3°) Ascending Series & Complex Continued Fraction
 
 

1°) Ascending Series
 

-Fresnel integrals are defined by:    S(x) = §₀^x  sin(pi.t²/2).dt      and      C(x) = §₀^x  cos(pi.t²/2).dt

Formulae:

     S(x) = SUM_{n=0,1,2,.....}  (-1)ⁿ (pi/2)²ⁿ⁺¹ x⁴ⁿ⁺³ / ((2n+1)!(4n+3))

     C(x) = SUM_{n=0,1,2,.....}  (-1)ⁿ (pi/2)²ⁿ x⁴ⁿ⁺¹ / ((2n)!(4n+1))
 

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

01  LBL "SX" 02  STO 00 03  3 04  Y^X 05  PI 06  6 07  / 08  * 09  1 10  STO 02 11  GTO 00

12  LBL "CX" 13  STO 00 14  0 15  STO 02 16  LBL 00 17  RCL 00 18  X^2 19  PI 20  * 21  2 22  /

23  X^2 24  CHS 25  STO 03     26  R^ 27  STO 01 28  LBL 01 29  RCL 01 30  RCL 02 31  2 32  + 33  STO 02

34  ST+ X 35  3 36  - 37  ST* Y      38  4 39  + 40  RCL 02 41  ST* Y 42  DSE X 43  * 44  /

45  RCL 03 46  * 47  STO 01     48  + 49  X#Y? 50  GTO 01 51  END

 
  ( 69 bytes / SIZE 004 )
 
 

STACK	INPUTS	OUTPUTS
X	x	S(x) or C(x)

Example:

    2  XEQ "SX"  >>>>  S(2) = 0.343415677    ( in 14 seconds )
    2  XEQ "CX" >>>>  C(2) = 0.488253412      --------------
 

Note:   This method does not give accurate results for large arguments:

    S(3) and C(3) are correct to 5 decimals
 S(3.5) and C(3.5) -----------  4 --------
    S(4) and C(4) are not exact even to 1D!
 

2°) Asymptotic Series ( large arguments )
 

Formulae:       S(x) = 1/2 - f(x).cos(pi.x²/2) - g(x).sin(pi.x²/2)     &     C(x) = 1/2 + f(x).sin(pi.x²/2) - g(x).cos(pi.x²/2)

   where               f(x) ~  ( 1 - 1*3/(pi.x²)² + 1*3*5*7/(pi.x²)⁴ - ....... ) / (pi.x)
     and                g(x) ~  ( 1 - 1*3*5/(pi.x²)² + 1*3*5*7*9/(pi.x²)⁴ - ....... ) / (pi².x³)
 

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

01  LBL "FSC" 02  STO 00 03  X^2 04  PI 05  * 06  STO 03 07  SIGN 08  STO 01 09  LASTX 10  1/X 11  STO 02 12  XEQ 01 13  STO 04

14  1 15  CHS 16  STO 01      17  CHS 18  STO 02 19  XEQ 01 20  PI 21  RCL 00 22  * 23  ST/ 04 24  / 25  RCL 03 26  2

27  / 28  STO 03 29  X<>Y 30  RAD 31  P-R 32  RCL 03      33  RCL 04 34  P-R 35  ST- T 36  RDN 37  + 38  .5 39  ST+ Z

40  X<>Y 41  - 42  DEG 43  RTN 44  LBL 01 45  RCL 02      46  RCL 01 47  2 48  + 49  ST* Y 50  LASTX 51  + 52  STO 01

53  * 54  RCL 03 55  X^2 56  / 57  CHS 58  STO 02      59  + 60  X#Y? 61  GTO 01 62  END

 
   ( 80 bytes / SIZE 005 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	C(x)
X	x	S(x)

Examples:

     10   XEQ FSC  >>>>    S(10) = 0.468169979    X<>Y    C(10) =  0.499898695    ( in 5 seconds )
     4.1  XEQ FSC  >>>>   S(4.1) = 0.475798258    X<>Y    C(4.1) = 0.573695632    ( in 14 seconds )
 

-The series are already divergent too soon for  x = 4   but if you add     RND  X<>Y   RND   X<>Y   after line 59, you'll get:

   in    FIX 9     4   XEQ "FSC"   >>>>   S(4) = 0.420515754   X<>Y   C(4) = 0.498426033   ( in 11 seconds )
   in    FIX 5     3   XEQ "FSC"   >>>>   S(3) = 0.4963140       X<>Y   S(3) = 0.6057203   ( exact values are  0.4963130  and  0.6057208 to 7D )
 

3°) Ascending Series & Continued Fraction
 

-The following program uses the same series expansions as "SX" & "CX"  if  | x |  <  SQRT(PI)
-Otherwise, a complex continued fraction is employed:

          C(x) + i.S(x) = ((1+i)/2) erf z   with   z = (1-i).(x/2).(pi)^1/2

 and   ( 1 - erf z ) exp z² =  (pi) ^-1/2 ( 1/(z+0.5/(z+1/(z+1.5/(z+2/(z+ .... ))))) )
 
 

Data Registers: /
Flags: /
Subroutines: /

-Synthetic registers M N O P may of course be replaced by registers R01 R02 R03 R00  ( for instance )
-In this case, delete line 122 , replace lines 110-111 by 3 , add  CLX  STO 02  after line 93 and delete line 02.
-Line 08 is a three-byte GTO 04
 
 

01  LBL "FRI"   02  CLA   03  STO M   04  ABS   05  PI   06  SQRT   07  X>Y?   08  GTO 04   09  DEG   10  *   11  24   12  STO N   13  CLX   14  2   15  /   16  STO O   17  ENTER^   18  LBL 01   19  X<>Y   20  CHS   21  STO Z   22  X^2   23  RCL Y   24  X^2   25  +

26  DSE N   27  X<0?   28  GTO 02   29  RCL N   30  X<>Y   31  ST+ X   32  /   33  ST* Z   34  *   35  RCL O   36  ST+ Z   37  +   38  GTO 01       39  LBL 02   40  PI   41  SQRT   42  *   43  ST/ Z   44  /   45  R-P   46  RCL M   47  X^2   48  90   49  *   50  R^

51  -   52  X<>Y   53  CHS   54  P-R   55  1   56  +   57  STO Z   58  X<>Y   59  ST+ Z   60  -   61  2   62  ST/ Z   63  /   64  GTO 05       65  LBL 03   66  X<> P   67  PI   68  2   69  ST+ N   70  /   71  RCL M   72  X^2   73  *   74  X^2   75  *

76  RCL N   77  ENTER^   78  DSE X   79  *   80  /   81  CHS   82  STO P   83  RCL N   84  ST+ X   85  ISG X   86  CLX   87  /   88  X<>Y   89  ST+ Y   90  X#Y?   91  GTO 03       92  RTN   93  LBL 04   94  X<>Y   95  STO P   96  ENTER^   97  XEQ 03   98  STO O   99  1 100  STO N

101  RCL M 102  ABS 103  3 104  Y^X 105  PI 106  * 107  2 108  / 109  STO P 110  PI 111  INT 112  / 113  ENTER^ 114  XEQ 03     115  RCL O 116  LBL 05 117  RCL M 118  SIGN 119  ST* Y 120  ST* Z 121  RDN 122  CLA 123  END

 
   ( 178 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	S(x)
X	x	C(x)
L	/	x

 

Examples:

    1.5  XEQ "FRI"  >>>>   C(1.5) = 0.445261176      X<>Y     S(1.5) = 0.697504959         ( in 22 seconds )
      4       R/S          >>>>     C(4)  = 0.498426033      X<>Y       S(4)  = 0.420515754         ( in 19 seconds )
 

-Thus, "FRI" produces accurate results for all arguments.
-However, "FSC" runs faster if x is very large.
 
 

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