The Museum of HP Calculators

Fresnel Integrals for the HP-41

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°) Ascending Series
2°) Asymptotic Series
3°) Ascending Series & Complex Continued Fraction ( new )

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 doesn't 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)
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:

FIX 9 4 XEQ "FSC" >>>> S(4) = 0.420515754 X<>Y C(4) = 0.498426033 ( in 11 seconds )
and 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.

01 LBL "FRI"
02 CLA
03 STO M
04 ABS
05 PI
06 SQRT
07 X>Y?
08 GTO 04     ( a three-byte GTO )
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        ( synthetic )
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        ( synthetic )
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        ( synthetic )
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: Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
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