Template

Carlson Elliptic Integrals for the HP-41

Overview

1°) RF+RFZ+RJ+RJZ
2°) RF+RFZ+RJ+RJZ ( shorter )
3°) RF+RFZ+RJ+RJZ ( slightly shorter )
4°) Slightly Different Formulae

a) Program#1
b) Program#2

5°) Shorter Program RF+RJ ( Exp-Sinh Quadrature )
6°) M-Code Routines

a) RF+RFZ
b) RJ+RJZ ( 2 versions )

Latest Update: paragraph 6°)b)

-Several programs are already listed in "Elliptic Integrals for the HP41"
-The routines in paragraphs 1°) to 4°) apply formulas in reference [2] & [3] when 2 arguments are conjugate complex.

-The different RF RFZ RJ RJZ are put together to save bytes.

-Of course, the M-Code routines in "Elliptic Integrals for the HP41" (§2°)e)) are much faster, but take also a look at paragraph 6°) ...

1°) RF+RFZ+RJ+RJZ

-We have: R_F(x;y;z) = (1/2) §₀^infinity ( ( t + x ).( t + y ).( t + z ) )^-1/2 dt with x , y , z non-negative and at most one is zero

and R_J(x;y;z;p) = (3/2) §₀^infinity ( t + p )^-1 ( ( t + x ).( t + y ).( t + z ) )^-1/2 dt with x , y , z non-negative and at most one is zero p > 0

-Cf "Elliptic Integrals for the HP41" if you want to calculate Cauchy principal value ( if p < 0 ).

-To compute I = R_F( x , y+i.z , y-i.z ) this program employs:

I = R_F( u , v , w ) with:

2 u = y + ( y² + z² )^1/2
2 v = x + ( y² + z² )^1/2 + [ (y-x)² + z² ]^1/2
2 w = x + ( y² + z² )^1/2 - [ (y-x)² + z² ]^1/2

( cf reference [2] )

Data Registers: R00 to R09: temp
Flag: F07
Subroutines: /

01 LBL "RFZ"
02 STO 01
03 X<>Y
04 STO 02
05 -
06 X<>Y
07 STO 03
08 R-P
09 ENTER
10 CHS
11 RCL 01
12 ST+ Z
13 +
14 RCL 02
15 RCL 03
16 R-P
17 X<>Y
18 RDN
19 ST+ Y
20 ST+ Z
21 RCL 02
22 +
23 2
24 ST/ Z
25 ST/ T
26 /
27 LBL "RF"
28 STO 01
29 RDN
30 STO 02
31 X<>Y
32 STO 03
33 CLX
34 LBL 01

35 RCL 02
36 SQRT
37 RCL X
38 RCL 03
39 SQRT
40 ST* Z
41 +
42 RCL 01
43 SQRT
44 *
45 +
46 ST+ 01
47 ST+ 02
48 ST+ 03
49 4
50 ST/ 01
51 ST/ 02
52 ST/ 03
53 R^
54 3
55 RCL 01
56 RCL 02
57 +
58 RCL 03
59 +
60 /
61 SQRT
62 X>Y?
63 GTO 01
64 RTN
65 LBL "RJ"
66 XEQ 09
67 RAD
68 LBL 02

69 RCL 01
70 SQRT
71 ENTER
72 STO 05
73 RCL 02
74 SQRT
75 ST* Z
76 STO 06
77 RCL 03
78 SQRT
79 ST* T
80 ST* 06
81 +
82 ST* 05
83 +
84 RCL 04
85 *
86 +
87 X^2
88 RCL 04
89 RCL 05
90 RCL 06
91 +
92 ST+ 01
93 ST+ 02
94 ST+ 03
95 RCL 04
96 +
97 STO 04
98 X^2
99 *
100 RCL 03
101 XEQ 08
102 X>Y?

103 GTO 02
104 DEG
105 RTN
106 LBL "RJZ"
107 XEQ 09
108 SF 07
109 RAD
110 LBL 03
111 RCL 03
112 RCL 02
113 R-P
114 ENTER
115 STO Z
116 LASTX
117 +
118 ST+ X
119 SQRT
120 STO 05
121 RCL 01
122 SQRT
123 ST* T
124 ST+ 05
125 *
126 +
127 ST+ 01
128 ST+ 02
129 X<> 04
130 ST+ 04
131 ENTER
132 X<> 05
133 *
134 +
135 X^2
136 RCL 05

137 RCL 04
138 ENTER
139 *
140 *
141 RCL 02
142 XEQ 08
143 ST- Y
144 /
145 ABS
146 2 E-9
147 X<Y?
148 GTO 03
149 FS?C 07
150 GTO 03
151 RCL 00
152 DEG
153 RTN
154 LBL 08
155 STO 06
156 RDN
157 X=Y?
158 GTO 04
159 -
160 STO 05
161 X<>Y
162 /
163 X<0?
164 GTO 05
165 SQRT
166 ENTER
167 ST+ Y
168 SIGN
169 LASTX

170 -
171 /
172 LN1+X
173 2
174 /
175 GTO 06
176 LBL 04
177 SQRT
178 1/X
179 GTO 07
180 LBL 05
181 CHS
182 SQRT
183 ATAN
184 LBL 06
185 RCL 05
186 ABS
187 SQRT
188 /
189 LBL 07
190 RCL 07
191 *
192 ST+ 08
193 4
194 ST/ 01
195 ST/ 02
196 ST/ 03
197 ST/ 04
198 ST/ 06
199 ST/ 07
200 RCL 00
201 RCL 07
202 RCL 01

203 RCL 02
204 +
205 RCL 04
206 ST+ X
207 +
208 RCL 06
209 +
210 5
211 /
212 ENTER
213 SQRT
214 *
215 /
216 RCL 08
217 3
218 *
219 +
220 STO 00
221 RTN
222 LBL 09
223 STO 01
224 RDN
225 STO 02
226 RDN
227 STO 03
228 X<>Y
229 STO 04
230 CLX
231 STO 00
232 STO 08
233 SIGN
234 STO 07
235 END

( 319 bytes / SIZE 009 )

STACK	INPUTS	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	R_F or R_J

p is only useful to compute R_J

Examples:

4 ENTER^
2 ENTER^
1 XEQ "RF" >>>> R_F( 1 , 2 , 4 ) = 0.685085816 ---Exectution time = 11s---

4 ENTER^
3 ENTER^
2 XEQ "RFZ" >>>>   R_F( 2 , 3+4.i , 3-4.i ) = 0.543421944    ---Exectution time = 12s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJ" >>>>   R_J( 1 , 2 , 3 , 4 ) = 0.239848100 ---Exectution time = 25s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJZ" >>>>   R_J( 1 , 2+3.i , 2-3.i , 4 ) = 0.205644141 ---Exectution time = 27s---

Note:

-If these programs are called as subroutines, registers R00-R01 .... may be replaced for instance by R40-R41 ..... to avoid registers conflicts

Data Registers: R40 to R49: temp
Flag: F07
Subroutines: /

01 LBL "RFZ"
02 STO 41
03 X<>Y
04 STO 42
05 -
06 X<>Y
07 STO 43
08 R-P
09 ENTER
10 CHS
11 RCL 41
12 ST+ Z
13 +
14 RCL 42
15 RCL 43
16 R-P
17 X<>Y
18 RDN
19 ST+ Y
20 ST+ Z
21 RCL 42
22 +
23 2
24 ST/ Z
25 ST/ T
26 /
27 LBL "RF"
28 STO 41
29 RDN
30 STO 42
31 X<>Y
32 STO 43
33 CLX
34 LBL 01

35 RCL 42
36 SQRT
37 RCL X
38 RCL 43
39 SQRT
40 ST* Z
41 +
42 RCL 41
43 SQRT
44 *
45 +
46 ST+ 41
47 ST+ 42
48 ST+ 43
49 4
50 ST/ 41
51 ST/ 42
52 ST/ 43
53 R^
54 3
55 RCL 41
56 RCL 42
57 +
58 RCL 43
59 +
60 /
61 SQRT
62 X>Y?
63 GTO 01
64 RTN
65 LBL "RJ"
66 XEQ 09
67 RAD
68 LBL 02

69 RCL 41
70 SQRT
71 ENTER
72 STO 45
73 RCL 42
74 SQRT
75 ST* Z
76 STO 46
77 RCL 43
78 SQRT
79 ST* T
80 ST* 46
81 +
82 ST* 45
83 +
84 RCL 44
85 *
86 +
87 X^2
88 RCL 44
89 RCL 45
90 RCL 46
91 +
92 ST+ 41
93 ST+ 42
94 ST+ 43
95 RCL 44
96 +
97 STO 44
98 X^2
99 *
100 RCL 43
101 XEQ 08
102 X>Y?

103 GTO 02
104 DEG
105 RTN
106 LBL "RJZ"
107 XEQ 09
108 SF 07
109 RAD
110 LBL 03
111 RCL 43
112 RCL 42
113 R-P
114 ENTER
115 STO Z
116 LASTX
117 +
118 ST+ X
119 SQRT
120 STO 45
121 RCL 41
122 SQRT
123 ST* T
124 ST+ 45
125 *
126 +
127 ST+ 41
128 ST+ 42
129 X<> 44
130 ST+ 44
131 ENTER
132 X<> 45
133 *
134 +
135 X^2
136 RCL 45

137 RCL 44
138 ENTER
139 *
140 *
141 RCL 42
142 XEQ 08
143 ST- Y
144 /
145 ABS
146 2 E-9
147 X<Y?
148 GTO 03
149 FS?C 07
150 GTO 03
151 RCL 40
152 DEG
153 RTN
154 LBL 08
155 STO 46
156 RDN
157 X=Y?
158 GTO 04
159 -
160 STO 45
161 X<>Y
162 /
163 X<0?
164 GTO 05
165 SQRT
166 ENTER
167 ST+ Y
168 SIGN
169 LASTX

170 -
171 /
172 LN1+X
173 2
174 /
175 GTO 06
176 LBL 04
177 SQRT
178 1/X
179 GTO 07
180 LBL 05
181 CHS
182 SQRT
183 ATAN
184 LBL 06
185 RCL 45
186 ABS
187 SQRT
188 /
189 LBL 07
190 RCL 47
191 *
192 ST+ 48
193 4
194 ST/ 41
195 ST/ 42
196 ST/ 43
197 ST/ 44
198 ST/ 46
199 ST/ 47
200 RCL 40
201 RCL 47
202 RCL 41

203 RCL 42
204 +
205 RCL 44
206 ST+ X
207 +
208 RCL 46
209 +
210 5
211 /
212 ENTER
213 SQRT
214 *
215 /
216 RCL 48
217 3
218 *
219 +
220 STO 40
221 RTN
222 LBL 09
223 STO 41
224 RDN
225 STO 42
226 RDN
227 STO 43
228 X<>Y
229 STO 44
230 CLX
231 STO 40
232 STO 48
233 SIGN
234 STO 47
235 END

2°) RF+RFZ+RJ+RJZ ( shorter program )

-We can save several bytes if we compute R_F( a , b , b ) with "RF" instead of the formulae:

     R_F(a,b,b) = ( a - b )^-1/2 Arctanh ( 1-b/a )^1/2     if   a > b
     R_F(a,b,b) =   a ^-1/2                                             if   a = b
     R_F(a,b,b) = ( b - a )^-1/2 Arctan ( b/a - 1 )^1/2      if   a < b

-But of course, "RJ" will be slower.

Data Registers: R00 to R10: temp
Flag: F07
Subroutines: /

29 STO 01
30 RDN
31 STO 02
32 X<>Y
33 STO 03
34 CLX
35 LBL 01
36 RCL 02
37 SQRT
38 RCL X
39 RCL 03
40 SQRT
41 ST* Z
42 +
43 RCL 01
44 SQRT
45 *
46 +
47 ST+ 01
48 ST+ 02
49 ST+ 03
50 4
51 ST/ 01
52 ST/ 02
53 ST/ 03
54 R^
55 3
56 RCL 01

57 RCL 02
58 +
59 RCL 03
60 +
61 /
62 SQRT
63 X>Y?
64 GTO 01
65 RTN
66 LBL "RJ"
67 XEQ 14
68 LBL 02
69 RCL 04
70 SQRT
71 ENTER
72 STO 01
73 RCL 05
74 SQRT
75 ST* Z
76 STO 02
77 RCL 06
78 SQRT
79 ST* T
80 ST* 02
81 +
82 ST* 01
83 +
84 RCL 07

85 *
86 +
87 X^2
88 RCL 07
89 RCL 01
90 RCL 02
91 +
92 ST+ 04
93 ST+ 05
94 ST+ 06
95 RCL 07
96 +
97 STO 07
98 RCL 06
99 XEQ 12
100 X>Y?
101 GTO 02
102 RTN
103 LBL "RJZ"
104 XEQ 14
105 SF 07
106 LBL 03
107 RCL 06
108 RCL 05
109 R-P
110 ENTER
111 STO Z
112 LASTX

113 +
114 ST+ X
115 SQRT
116 STO 01
117 RCL 04
118 SQRT
119 ST* T
120 ST+ 01
121 *
122 +
123 ST+ 04
124 ST+ 05
125 X<> 07
126 ST+ 07
127 ENTER
128 X<> 01
129 *
130 +
131 X^2
132 RCL 01
133 RCL 07
134 RCL 05
135 XEQ 12
136 ST- Y
137 /
138 ABS
139 2 E-9
140 X<Y?

141 GTO 03
142 FS?C 07
143 GTO 03
144 RCL 00
145 RTN
146 LBL 12
147 STO 10
148 RDN
149 X^2
150 *
151 ENTER
152 XEQ 13
153 RCL 09
154 *
155 ST+ 08
156 4
157 ST/ 04
158 ST/ 05
159 ST/ 06
160 ST/ 07
161 ST/ 09
162 ST/ 10
163 RCL 00
164 RCL 09
165 RCL 04
166 RCL 05
167 +
168 RCL 07
169 ST+ X

170 +
171 RCL 10
172 +
173 5
174 /
175 ENTER
176 SQRT
177 *
178 /
179 RCL 08
180 3
181 *
182 +
183 STO 00
184 RTN
185 LBL 14
186 STO 04
187 RDN
188 STO 05
189 RDN
190 STO 06
191 X<>Y
192 STO 07
193 CLX
194 STO 00
195 STO 08
196 SIGN
197 STO 09
198 END

( 279 bytes / SIZE 011 )

STACK	INPUTS	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	R_F or R_J

p is only useful to compute R_J

Examples:

    4 ENTER^
2 ENTER^
1 XEQ "RF" >>>> R_F( 1 , 2 , 4 ) = 0.685085816 ---Exectution time = 11s---

4 ENTER^
3 ENTER^
2 XEQ "RFZ" >>>>   R_F( 2 , 3+4.i , 3-4.i ) = 0.543421944    ---Exectution time = 12s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJ" >>>>   R_J( 1 , 2 , 3 , 4 ) = 0.239848100 ---Exectution time = 39s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJZ" >>>>   R_J( 1 , 2+3.i , 2-3.i , 4 ) = 0.205644140 ---Exectution time = 52s---

3°) RF+RFZ+RJ+RJZ ( slightly shorter program )

-We can again save a few bytes if we use the formulae given in reference [3] to compute I = R_J( x , y+i.z , y-i.z , p )

I = [ s / ( s² - u p ) ] [ 2 ( s - u ) R_J( u , v , w , s ) - 3 R_F( x , y+i.z , y-i.z ) + 3 R_F( x , p , p ) }] with

2 u = y + ( y² + z² )^1/2
2 v = x + ( y² + z² )^1/2 + [ (y-x)² + z² ]^1/2
2 w = x + ( y² + z² )^1/2 - [ (y-x)² + z² ]^1/2

2 s = p + ( y² + z² )^1/2 + [ (y-p)² + z² ]^1/2

Data Registers: R00 to R13: temp
Flags: /
Subroutines: /

01 LBL "RFZ"
02 LBL 12
03 STO 01
04 X<>Y
05 STO 02
06 -
07 X<>Y
08 STO 03
09 R-P
10 ENTER
11 CHS
12 RCL 01
13 ST+ Z
14 +
15 RCL 02
16 RCL 03
17 R-P
18 STO 00
19 X<>Y
20 RDN
21 ST+ Y
22 ST+ Z
23 RCL 02
24 +
25 2
26 ST/ Z
27 ST/ T
28 /
29 STO 08
30 RDN

31 STO 09
32 X<>Y
33 STO 10
34 X<>Y
35 R^
36 LBL "RF"
37 LBL 13
38 STO 01
39 RDN
40 STO 02
41 X<>Y
42 STO 03
43 CLX
44 LBL 01
45 RCL 02
46 SQRT
47 RCL X
48 RCL 03
49 SQRT
50 ST* Z
51 +
52 RCL 01
53 SQRT
54 *
55 +
56 ST+ 01
57 ST+ 02
58 ST+ 03
59 4
60 ST/ 01

61 ST/ 02
62 ST/ 03
63 R^
64 3
65 RCL 01
66 RCL 02
67 +
68 RCL 03
69 +
70 /
71 SQRT
72 X>Y?
73 GTO 01
74 RTN
75 LBL "RJZ"
76 XEQ 14
77 RCL 07
78 RCL 07
79 RCL 04
80 XEQ 13
81 STO 11
82 RCL 06
83 RCL 05
84 RCL 04
85 XEQ 12
86 RCL 11
87 -
88 3
89 *
90 STO 11

91 RCL 05
92 RCL 07
93 ST+ 00
94 -
95 RCL 06
96 R-P
97 RCL 00
98 +
99 2
100 /
101 STO 12
102 ENTER
103 X^2
104 RCL 08
105 RCL 07
106 *
107 -
108 /
109 STO 13
110 RCL 12
111 RCL 10
112 RCL 09
113 RCL 08
114 ST- 12
115 XEQ 10
116 RCL 12
117 ST+ X
118 *
119 RCL 11
120 -

121 RCL 13
122 *
123 RTN
124 LBL "RJ"
125 LBL 10
126 XEQ 14
127 LBL 02
128 RCL 04
129 SQRT
130 ENTER
131 STO 01
132 RCL 05
133 SQRT
134 ST* Z
135 STO 02
136 RCL 06
137 SQRT
138 ST* T
139 ST* 02
140 +
141 ST* 01
142 +
143 RCL 07
144 *
145 +
146 X^2
147 RCL 07
148 RCL 01
149 RCL 02

150 +
151 ST+ 04
152 ST+ 05
153 ST+ 06
154 RCL 07
155 +
156 STO 07
157 X^2
158 *
159 ENTER
160 XEQ 13
161 RCL 09
162 *
163 ST+ 08
164 4
165 ST/ 04
166 ST/ 05
167 ST/ 06
168 ST/ 07
169 ST/ 09
170 RCL 00
171 RCL 09
172 RCL 04
173 RCL 05
174 +
175 RCL 06
176 +
177 RCL 07
178 ST+ X

179 +
180 5
181 /
182 ENTER
183 SQRT
184 *
185 /
186 RCL 08
187 3
188 *
189 +
190 STO 00
191 X>Y?
192 GTO 02
193 RTN
194 LBL 14
195 STO 04
196 RDN
197 STO 05
198 RDN
199 STO 06
200 X<>Y
201 STO 07
202 CLX
203 STO 00
204 STO 08
205 SIGN
206 STO 09
207 END

( 275 bytes / SIZE 014 )

STACK	INPUTS	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	R_F or R_J

p is only useful to compute R_J

Examples:

    4 ENTER^
2 ENTER^
1 XEQ "RF" >>>> R_F( 1 , 2 , 4 ) = 0.685085816 ---Exectution time = 11s---

4 ENTER^
3 ENTER^
2 XEQ "RFZ" >>>>   R_F( 2 , 3+4.i , 3-4.i ) = 0.543421944    ---Exectution time = 12s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJ" >>>>   R_J( 1 , 2 , 3 , 4 ) = 0.239848100 ---Exectution time = 38s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJZ" >>>>   R_J( 1 , 2+3.i , 2-3.i , 4 ) = 0.205644140 ---Exectution time = 69s---

Note:

"RJZ" is much slower but we have saved 4 bytes (!)

4°) Slightly Different Formulae

a) Program#1

-In paragraphs 1°) 2°) 3°) we use the formulae:

R_F(x;y;z) = R_F((x+L)/4;(y+L)/4;(z+L)/4) where L = x^1/2y^1/2 + x^1/2z^1/2 + y^1/2z^1/2

and R_J(x;y;z;p) = 3 R_F(a,b,b) + (1/4) R_J((x+L)/4;(y+L)/4;(z+L)/4,(p+L)/4) where a = ( p( x^1/2 + y^1/2 + z^1/2 ) + ( xyz )^1/2 )² ; b = p ( p + L )²

-In paragraph 4°) , the following formulas are employed:

R_F(x;y;z) = 2 R_F(x+L;y+L;z+L) and

R_J(x;y;z;p) = 3 R_F(a,b,b) + 2 R_J(x+L;y+L;z+L,p+L) 12 bytes and a few seconds are saved !

Data Registers: R00 to R08: temp
Flag: F07
Subroutines: /

34 STO 00
35 CLX
36 LBL 01
37 RCL 02
38 SQRT
39 RCL X
40 RCL 03
41 SQRT
42 ST* Z
43 +
44 RCL 01
45 SQRT
46 *
47 +
48 ST+ 01
49 ST+ 02
50 ST+ 03
51 R^
52 3
53 RCL 01
54 RCL 02
55 +
56 RCL 03
57 +
58 /
59 SQRT
60 RCL 00
61 ST+ 00
62 *
63 X>Y?
64 GTO 01
65 RTN
66 LBL "RJ"

67 XEQ 09
68 RAD
69 LBL 02
70 RCL 01
71 SQRT
72 ENTER
73 STO 05
74 RCL 02
75 SQRT
76 ST* Z
77 STO 06
78 RCL 03
79 SQRT
80 ST* T
81 ST* 06
82 +
83 ST* 05
84 +
85 RCL 04
86 *
87 +
88 X^2
89 RCL 04
90 RCL 05
91 RCL 06
92 +
93 ST+ 01
94 ST+ 02
95 ST+ 03
96 RCL 04
97 +
98 STO 04
99 X^2

100 *
101 RCL 03
102 XEQ 08
103 X>Y?
104 GTO 02
105 DEG
106 RTN
107 LBL "RJZ"
108 XEQ 09
109 SF 07
110 RAD
111 LBL 03
112 RCL 03
113 RCL 02
114 R-P
115 ENTER
116 STO Z
117 LASTX
118 +
119 ST+ X
120 SQRT
121 STO 05
122 RCL 01
123 SQRT
124 ST* T
125 ST+ 05
126 *
127 +
128 ST+ 01
129 ST+ 02
130 X<> 04
131 ST+ 04
132 ENTER

133 X<> 05
134 *
135 +
136 X^2
137 RCL 05
138 RCL 04
139 ENTER
140 *
141 *
142 RCL 02
143 XEQ 08
144 ST- Y
145 /
146 ABS
147 2 E-9
148 X<Y?
149 GTO 03
150 FS?C 07
151 GTO 03
152 RCL 00
153 DEG
154 RTN
155 LBL 08
156 STO 06
157 RDN
158 X=Y?
159 GTO 04
160 -
161 STO 05
162 X<>Y
163 /
164 X<0?
165 GTO 05

166 SQRT
167 ENTER
168 ST+ Y
169 SIGN
170 LASTX
171 -
172 /
173 LN1+X
174 2
175 /
176 GTO 06
177 LBL 04
178 SQRT
179 1/X
180 GTO 07
181 LBL 05
182 CHS
183 SQRT
184 ATAN
185 LBL 06
186 RCL 05
187 ABS
188 SQRT
189 /
190 LBL 07
191 RCL 07
192 ST+ 07
193 *
194 ST+ 08
195 RCL 00
196 RCL 07
197 RCL 01
198 RCL 02

199 +
200 RCL 04
201 ST+ X
202 +
203 RCL 06
204 +
205 5
206 /
207 ENTER
208 SQRT
209 *
210 /
211 RCL 08
212 3
213 *
214 +
215 STO 00
216 RTN
217 LBL 09
218 STO 01
219 RDN
220 STO 02
221 RDN
222 STO 03
223 X<>Y
224 STO 04
225 CLX
226 STO 00
227 STO 08
228 SIGN
229 STO 07
230 END

( 307 bytes / SIZE 009 )

STACK	INPUTS	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	R_F or R_J

p is only useful to compute R_J

Examples:

4 ENTER^
2 ENTER^
1 XEQ "RF" >>>> R_F( 1 , 2 , 4 ) = 0.685085816 ---Exectution time = 9s---

4 ENTER^
3 ENTER^
2 XEQ "RFZ" >>>>   R_F( 2 , 3+4.i , 3-4.i ) = 0.543421944    ---Exectution time = 12s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJ" >>>>   R_J( 1 , 2 , 3 , 4 ) = 0.239848100 ---Exectution time = 21s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJZ" >>>>   R_J( 1 , 2+3.i , 2-3.i , 4 ) = 0.205644141 ---Exectution time = 23s---

Note:

-This program corresponds to the program listed in §1°)
-Likewise, we can modify the routine listed in paragraph 2°) with the same modified formulas:

b) Program#2

Data Registers: R00 to R11: temp
Flag: F07
Subroutines: /

29 STO 01
30 RDN
31 STO 02
32 X<>Y
33 STO 03
34 2
35 STO 00
36 CLX
37 LBL 01
38 RCL 02
39 SQRT
40 RCL X
41 RCL 03
42 SQRT
43 ST* Z
44 +
45 RCL 01
46 SQRT
47 *
48 +
49 ST+ 01
50 ST+ 02
51 ST+ 03
52 R^
53 3
54 RCL 01
55 RCL 02
56 +

57 RCL 03
58 +
59 /
60 SQRT
61 RCL 00
62 ST+ 00
63 *
64 X>Y?
65 GTO 01
66 RTN
67 LBL "RJ"
68 XEQ 14
69 LBL 02
70 RCL 04
71 SQRT
72 ENTER
73 STO 01
74 RCL 05
75 SQRT
76 ST* Z
77 STO 02
78 RCL 06
79 SQRT
80 ST* T
81 ST* 02
82 +
83 ST* 01
84 +

85 RCL 07
86 *
87 +
88 X^2
89 RCL 07
90 RCL 01
91 RCL 02
92 +
93 ST+ 04
94 ST+ 05
95 ST+ 06
96 RCL 07
97 +
98 STO 07
99 RCL 06
100 XEQ 12
101 X>Y?
102 GTO 02
103 RTN
104 LBL "RJZ"
105 XEQ 14
106 SF 07
107 LBL 03
108 RCL 06
109 RCL 05
110 R-P
111 ENTER
112 STO Z

113 LASTX
114 +
115 ST+ X
116 SQRT
117 STO 01
118 RCL 04
119 SQRT
120 ST* T
121 ST+ 01
122 *
123 +
124 ST+ 04
125 ST+ 05
126 X<> 07
127 ST+ 07
128 ENTER
129 X<> 01
130 *
131 +
132 X^2
133 RCL 01
134 RCL 07
135 RCL 05
136 XEQ 12
137 ST- Y
138 /
139 ABS

140 2 E-9
141 X<Y?
142 GTO 03
143 FS?C 07
144 GTO 03
145 RCL 11
146 RTN
147 LBL 12
148 STO 10
149 RDN
150 X^2
151 *
152 ENTER
153 XEQ 13
154 RCL 09
155 ST+ 09
156 *
157 ST+ 08
158 RCL 11
159 RCL 09
160 RCL 04
161 RCL 05
162 +
163 RCL 07
164 ST+ X
165 +
166 RCL 10

167 +
168 5
169 /
170 ENTER
171 SQRT
172 *
173 /
174 RCL 08
175 3
176 *
177 +
178 STO 11
179 RTN
180 LBL 14
181 STO 04
182 RDN
183 STO 05
184 RDN
185 STO 06
186 X<>Y
187 STO 07
188 CLX
189 STO 08
190 STO 11
191 SIGN
192 STO 09
193 END

( 267 bytes / SIZE 012 )

STACK	INPUTS	OUTPUT
T	p > 0	/
Z	z	/
Y	y	/
X	x	R_F or R_J

p is only useful to compute R_J

Examples:

4 ENTER^
2 ENTER^
1 XEQ "RF" >>>> R_F( 1 , 2 , 4 ) = 0.685085816 ---Exectution time = 9s---

4 ENTER^
3 ENTER^
2 XEQ "RFZ" >>>>   R_F( 2 , 3+4.i , 3-4.i ) = 0.543421944    ---Exectution time = 12s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJ" >>>>   R_J( 1 , 2 , 3 , 4 ) = 0.239848100 ---Exectution time = 40s---

4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "RJZ" >>>>   R_J( 1 , 2+3.i , 2-3.i , 4 ) = 0.205644141 ---Exectution time = 47s---

Note:

-The program listed in §3°) could be modified similarly...

5°) Shorter Program ( Exp-Sinh Quadrature )

-Here is the shortest program I've found to compute R_F & R_J
-The Exp-Sinh quadrature is used to compute the integrals.
-Unfortunately, this program is very slow...

  The precision is controlled by the display format.

Data Registers: R00 to R12: temp R08 = x R09 = y R10 = z R11 = p > 0

Flag: F01 CF 01 to calculate R_F(x;y;z)
SF 01 to calculate R_J(x;y;z,p)
Subroutines: /

01 LBL "RFJ"
02 STO 08
03 RDN
04 STO 09
05 RDN
06 STO 10
07 X<>Y
08 STO 11
09 CLX
10 STO 03
11 STO 07
12 SIGN
13 STO 02
14 STO 06
15 XEQ 10
16 RCL 02
17 *
18 STO 01
19 GTO 01
20 LBL 00
21 2

22 ST/ 01
23 ST/ 02
24 STO 06
25 SIGN
26 CHS
27 STO 03
28 LBL 01
29 RCL 06
30 ST+ 03
31 RCL 02
32 RCL 03
33 *
34 E^X-1
35 LASTX
36 CHS
37 E^X-1
38 -
39 2
40 /
41 E^X
42 STO 04

43 STO 05
44 XEQ 10
45 X<> 04
46 1/X
47 XEQ 10
48 RCL 04
49 RCL 05
50 ST/ Z
51 *
52 +
53 RCL 02
54 ST* Y
55 RCL 03
56 *
57 E^X
58 ENTER
59 1/X
60 +
61 *
62 2
63 /

64 RCL 01
65 +
66 STO 01
67 LASTX
68 X#Y?
69 GTO 01
70 VIEW X
71 X<> 07
72 RND
73 X<>Y
74 RND
75 X#Y?
76 GTO 00
77 RCL 07
78 .5
79 FS? 01
80 ISG X
81 *
82 CLD
83 RTN
84 LBL 10

85 STO 12
86 RCL 08
87 +
88 1/X
89 RCL 09
90 RCL 12
91 +
92 /
93 RCL 10
94 RCL 12
95 +
96 /
97 SQRT
98 FC? 01
99 RTN
100 RCL 11
101 RCL 12
102 +
103 /
104 END

( 133 bytes / SIZE 013 )

STACK	INPUTS	OUTPUTS
T	p > 0	/
Z	z	/
Y	y	/
X	x	R_F(x;y;z) or R_J(x;y;z,p)

Example1: CF 01 FIX 8

4 ENTER
2 ENTER
1 XEQ "RFJ" >>>> R_F(1;2;4) = 0.685085817

Example2: SF 01 FIX 9

4 ENTER
3 ENTER
2 ENTER
1 XEQ "RFJ" >>>> R_J(1;2;3;4) = 0.239848100

Notes:

-The execution time = more than 2 minutes !
-We can slightly reduce it if

   -we add RTN after line 103
-we replace lines 62-39-21 by RCL 13
-we add 2 STO 13 after line 08

6°) M-Code Routines

a) RF+RFZ

-"RFZ" employs the method given in paragraph 1
-"RF" uses the formula given in paragraph 4.

-Synthetic register Q is used.

Warning:

-They don't check for alpha data.
-If there are too much "zeros", it could make infinte loops !

-If you have a "newest" HP41, press ENTER ON to stop the infinite loop.

-Otherwise, replace the 2 lines ( near the end )

2FE ?C<0
263 JNC-52d by

3CC   ?KEY
360    ?C RTN
2FE ?C<0?
253 JNC-54d

-Thus, the infinite loop will be stopped by pressing any key for 1 second.

09A "Z"
006 "F"
012 "R"
2A0 SETDEC                  The first executable word of RFZ
078 C=Z
10E A=C ALL
135 C=
060 A*C
070 N=C ALL
0B8 C=Y
10E A=C ALL
135 C=
060 A*C
0B0 C=N ALL
025 C=
060 AB+C
305 C=
060 sqrt(AB)
268 Q=C
0B8 C=Y
025 C=
060 AB+C
04E C
35C =
090 2
269 C=
060 AB/C
128 L=C
0F8 C=X
10E A=C ALL
278 C=Q
01D C=
060 A+C
268 Q=C
0B8 C=Y
2BE C=-C
10E A=C ALL
0F8 C=X
01D C=
060 A+C
13D C=
060 AB*C
0B0 C=N ALL
025 C=
060 AB+C
305 C=
060 sqrt(AB)
070 N=C ALL
278 C=Q
025 C=
060 AB+C
04E C
35C =
090 2
269 C=
060 AB/C
0A8 Y=C
0B0 C=N ALL
2BE C=-C
10E A=C ALL
278 C=Q
01D C=
060 A+C
04E C
35C =
090 2
269 C=
060 AB/C
068 Z=C
138 C=L
0E8 X=C
023 JNC+04
086 "F"
012 "R"
2A0 SETDEC                    The first executable word of RF
04E C=0 ALL
268 Q=C
35C C=
090 2
028 T=C
0F8 C=X LOOP
10E A=C ALL
0B8 C=Y
135 C=
060 A*C
305 C=
060 sqrt(AB)
070 N=C ALL
0F8 C=X
10E A=C ALL
078 C=Z
135 C=
060 A*C
305 C=
060 sqrt(AB)
128 L=C
0B8 C=Y
10E A=C ALL
078 C=Z
135 C=
060 A*C
305 C=
060 sqrt(AB)
0B0 C=N ALL
025 C=
060 AB+C
138 C=L
025 C=
060 AB+C
070 N=C ALL
0F8 C=X
025 C=
060 AB+C
0E8 X=C
0B8 C=Y
10E A=C ALL
0B0 C=N ALL
01D C=
060 A+C
0A8 Y=C
078 C=Z
10E A=C ALL
0B0 C=N ALL
01D C=
060 A+C
068 Z=C
013 JNC+02
28B JNC-47d ( GTO LOOP )
0F8 C=X
025 C=
060 AB+C
0B8 C=Y
025 C=
060 AB+C
10E A=C
04E C
35C =
0D0 3
0AE A<>C ALL
261 C=
060 A/C
305 C=
060 sqrt(AB)
046 C
270 =
038 T
070 N=C ALL
13D C=
060 AB*C
128 L=C
0B0 C=N ALL
10E A=C ALL
01D C=
060 A+C
028 T=C
138 C=L
10E A=C ALL
278 C=Q
0AE A<>C ALL
268 Q=C
2BE C=-C
01D C=
061 A+C
278 C=Q
269 C=
060 AB/C
05E C=|C|
10E A=C ALL
04E C
35C
2BE
090
21C =
250
250
050 -2 10^(-9)
01D C= if you want to stop a possible infinte loop with an old HP41,
060 A+C
2FE ?C<0    replace 3CC   ?KEY 2FE ?C<0?
263 JNC-52d ( GTO LOOP ) these 2 words by 360    ?C RTN 253 JNC-54d
278 C=Q
0E8 X=C
3E0 RTN

( 183 words )

RF STACK	INPUTS	OUTPUTS
Z	z	/
Y	y	/
X	x	R_F(x,y,z)

Example:

    4 ENTER^
    2 ENTER^
    1 XEQ "RF" yields R_F(1,2,4) = 0.6850858167 ( in 4.5 seconds )

RFZ STACK	INPUTS	OUTPUTS
Z	z	/
Y	y	/
X	x	R_F(x,y+iz,y-iz)

Example:

    4 ENTER^
    3 ENTER^
    2 XEQ "RFZ" gives R_F ( 2 , 3+4i , 3-4i ) = 0.5434219446    ( in 5.2 seconds )

b) RJ+RJZ ( 2 Versions )

-I've remarked that if R_C is computed directly, it's faster but the results may be less accurate.
-So, these M-Code routines calculate R_C with R_F

-I've not found a method employing the formulae for R_J & R_F given in paragraph 4 without data register, so the 1st routine employs the usual formulae.

-All the alpha "registers" are used and cleared at the end.

p must be positive

Warning:

-They don't check for alpha data.
-If there are too much "zeros", it could make infinte loops !

-If you have a "newest" HP41, press ENTER ON to stop the infinite loop.
-Otherwise, modify this routine as suggested below.

09A "Z"
00A "J"
012 "R"
108 SETF 8                  The first executable word of RJZ
248 SETF 9
02B JNC+05
08A "J"
012 "R"
104 CLRF 8                  The first executable word of RJ
244 CLRF 9
2A0 SETDEC
046 C
270 =
038 T
228 P=C
04E C=0
028 T=C
268 Q=C
35C C=
050 1
1E8 O=C
0F8 C=X LOOP1
2F9 C=
060 sqrt(C)
168 M=C
0B8 C=Y
10C ?FSET8
0F7 JC+30d
2F9 C=
060 sqrt(C)
1A8 N=C
078 C=Z
2F9 C=
060 sqrt(C)
128 L=C
1B8 C=N
025 C=
060 AB+C
070 N=C ALL
178 C=M
13D C=
060 AB*C
0F0 C<>N ALL
10E A=C ALL
178 C=M
01D C=
060 A+C
10E A=C ALL
1B8 C=N
0AE A<>C ALL
1A8 N=C
138 C=L
135 C=
060 A*C
128 L=C
0B0 C=N ALL
11B JNC+35d
10E A=C ALL
135 C=
060 A*C
070 N=C ALL
078 C=Z
10E A=C ALL
135 C=
060 A*C
0B0 C=N ALL
025 C=
060 AB+C
305 C=
060 sqrt(AB)
128 L=C
0B8 C=Y
025 C=
060 AB+C
025 C=
060 AB+C
305 C=
060 sqrt(AB)
070 N=C ALL
178 C=M
025 C=
060 AB+C
1A8 N=C
178 C=M
013 JNC+02
203 JNC-64d
10E A=C ALL
0B0 C=N ALL
135 C=
060 A*C
138 C=L
025 C=
060 AB+C
070 N=C ALL
138 C=L
10E A=C ALL
178 C=M
135 C=
060 A*C
168 M=C
0B0 C=N ALL
10E A=C ALL
0F8 C=X
01D C=
060 A+C
04E C
35C =
110 4
128 L=C
269 C=
060 AB/C
0E8 X=C
0B0 C=N ALL
10E A=C ALL
0B8 C=Y
01D C=
060 A+C
138 C=L
269 C=
060 AB/C
0A8 Y=C
0B0 C=N ALL
10E A=C ALL
078 C=Z
10C ?FSET8
01D C=
060 A+C
10E A=C ALL
138 C=L
261 C=
060 A/C
068 Z=C
1B8 C=N
10E A=C ALL
238 C=P
135 C=
060 A*C
178 C=M
025 C=
060 AB+C
13D C=
060 AB*C
168 M=C
0B0 C=N ALL
013 JNC+02
223 JNC-60d GTO LOOP1
10E A=C ALL
238 C=P
01D C=
060 A+C
070 N=C ALL
13D C=
060 AB*C
238 C=P
13D C=
060 AB*C
1A8 N=C
0B0 C=N ALL
10E A=C ALL
138 C=L
261 C=
060 A/C
228 P=C
04E C=0
128 L=C
1B8 C=N LOOP2 this loop computes R_F(a,b,b) with M = a & N = b
10E A=C ALL
178 C=M
135 C=
060 A*C
305 C=
060 sqrt(AB)
025 C=
060 AB+C
1B8 C=N
025 C=
060 AB+C
070 N=C ALL
178 C=M
025 C=
060 AB+C
04E C
35C =
110 4
269 C=
060 AB/C
168 M=C
0B0 C=N ALL
10E A=C ALL
1B8 C=N
01D C=
060 A+C
04E C
35C =
110 4
269 C=
060 AB/C
1A8 N=C
025 C=
060 AB+C
178 C=M
025 C=
060 AB+C
10E A=C ALL
04E C
35C =
0D0 3
0AE A<>C ALL
013 JNC+02
203 JNC-64d GTO LOOP1
261 C=
060 A/C
305 C=
060 sqrt(AB)
10E A=C ALL
138 C=L
0AE A<>C ALL
128 L=C
2BE C=-C
01D C=
061 A+C
138 C=L
013 JNC+02
233 JNC-58d GTO LOOP2
269 C=
060 AB/C
05E C=|C|
10E A=C ALL
04E C
35C
2BE
090
21C =
250
250
050 -2 10^(-9) If you want to stop a possible infinte loop with an "old" HP41
01D C=
060 A+C    add    3CC   ?KEY 360    ?C RTN after this line
2FE ?C<0
383 JNC-16d GTO LOOP2 replace 383 JNC-16d by 373 JNC-18d and replace the next GTO LOOP1 by ( see below )
138 C=L C = R_F(a,b,b)
10E A=C ALL
04E C
35C =
0D0 3
135 C=
060 A*C
1F8 C=O
13D C=
060 AB*C
278 C=Q
025 C=
060 AB+C
268 Q=C
1F8 C=O
10E A=C ALL
04E C
35C =
110 4
261 C=
060 A/C
1E8 O=C
0F8 C=X
10E A=C ALL
0B8 C=Y
013 JNC+02
23B JNC-57d GTO LOOP1 or 22B JNC-59d if you have added 3CC   ?KEY    360    ?C RTN before the previous GTO LOOP2
01D C=
060 A+C
078 C=Z
10C ?FSET8
013 JNC+02
0B8 C=Y
025 C=
060 AB+C
238 C=P
025 C=
060 AB+C
238 C=P
025 C=
060 AB+C
04E C
35C =
150 5
269 C=
060 AB/C
128 L=C
305 C=
060 sqrt(AB)
138 C=L
13D C=
060 AB*C
10E A=C
1F8 C=O
0AE A<>C ALL
261 C=
060 A/C
278 C=Q
025 C=
060 AB+C
10E A=C ALL
046 C
270 =
038 T
0AE A<>C ALL
028 T=C
128 L=C
2BE C=-C
01D C=
061 A+C
138 C=L
269 C=
060 AB/C
05E C=|C|
10E A=C ALL
04E C
35C
2BE
090
21C =
250
250
050 -2 10^(-9)
01D C=
060 A+C
2FE ?C<0
223 JNC-60d GTO LOOP1
24C ?FSET9
01B JNC+03
244 CLRF9
203 JNC-64d GTO LOOP1
3B5 ?NCXQ
050 R^
345 ?NCGO
042 CLA

( 335 words )

RJ STACK	INPUTS	OUTPUTS
T	p>0	/
Z	z	/
Y	y	/
X	x	R_J(x;y;z;p)

Example:

     7 ENTER^
     4 ENTER^
     2 ENTER^
     1 XEQ "RJ" yields   R_J(1,2,4,7) = 0.1478544450 ( in 15 seconds )

RJZ STACK	INPUTS	OUTPUTS
T	p>0	/
Z	z	/
Y	y	/
X	x	R_J(x,y+iz,y-iz,p)

Example:

     7 ENTER^
     4 ENTER^
     2 ENTER^
     1 XEQ "RJZ" >>>> R_J ( 1 , 2+4i , 2-4i , 7 ) = 0.1260390651 ( in 17 seconds )

Note:

-The following routine employs the formulae given in paragraph 4°) for R_J , but not in the 2nd loop for R_F to calculate R_C

09A "Z"
00A "J"
012 "R"
108 SETF 8                  The first executable word of RJZ
248 SETF 9
02B JNC+05
08A "J"
012 "R"
104 CLRF 8                  The first executable word of RJ
244 CLRF 9
2A0 SETDEC
046 C
270 =
038 T
228 P=C
04E C=0
028 T=C
268 Q=C
35C C=
050 1
1E8 O=C
0F8 C=X LOOP1
2F9 C=
060 sqrt(C)
168 M=C
0B8 C=Y
10C ?FSET8
0F7 JC+30d
2F9 C=
060 sqrt(C)
1A8 N=C
078 C=Z
2F9 C=
060 sqrt(C)
128 L=C
1B8 C=N
025 C=
060 AB+C
070 N=C ALL
178 C=M
13D C=
060 AB*C
0F0 C<>N ALL
10E A=C ALL
178 C=M
01D C=
060 A+C
10E A=C ALL
1B8 C=N
0AE A<>C ALL
1A8 N=C
138 C=L
135 C=
060 A*C
128 L=C
0B0 C=N ALL
11B JNC+35d
10E A=C ALL
135 C=
060 A*C
070 N=C ALL
078 C=Z
10E A=C ALL
135 C=
060 A*C
0B0 C=N ALL
025 C=
060 AB+C
305 C=
060 sqrt(AB)
128 L=C
0B8 C=Y
025 C=
060 AB+C
025 C=
060 AB+C
305 C=
060 sqrt(AB)
070 N=C ALL
178 C=M
025 C=
060 AB+C
1A8 N=C
178 C=M
013 JNC+02
203 JNC-64d
10E A=C ALL
0B0 C=N ALL
135 C=
060 A*C
138 C=L
025 C=
060 AB+C
070 N=C ALL
138 C=L
10E A=C ALL
178 C=M
135 C=
060 A*C
168 M=C
0B0 C=N ALL
10E A=C ALL
0F8 C=X
01D C=
060 A+C
0E8 X=C
0B0 C=N ALL
10E A=C ALL
0B8 C=Y
01D C=
060 A+C
0A8 Y=C
0B0 C=N ALL
10E A=C ALL
078 C=Z
10C ?FSET8
01D C=
060 A+C
068 Z=C
1B8 C=N
10E A=C ALL
238 C=P
135 C=
060 A*C
178 C=M
025 C=
060 AB+C
13D C=
060 AB*C
168 M=C
0B0 C=N ALL
013 JNC+02
28B JNC-47d GTO LOOP1
10E A=C ALL
238 C=P
01D C=
060 A+C
070 N=C ALL
13D C=
060 AB*C
238 C=P
13D C=
060 AB*C
1A8 N=C
0B0 C=N ALL
228 P=C
04E C=0
128 L=C
1B8 C=N LOOP2 this loop computes R_F(a,b,b) with M = a & N = b
10E A=C ALL
178 C=M
135 C=
060 A*C
305 C=
060 sqrt(AB)
025 C=
060 AB+C
1B8 C=N
025 C=
060 AB+C
070 N=C ALL
178 C=M
025 C=
060 AB+C
04E C
35C =
110 4
269 C=
060 AB/C
168 M=C
0B0 C=N ALL
10E A=C ALL
1B8 C=N
01D C=
060 A+C
04E C
35C =
110 4
269 C=
060 AB/C
1A8 N=C
025 C=
060 AB+C
178 C=M
025 C=
060 AB+C
10E A=C ALL
04E C
35C =
0D0 3
0AE A<>C ALL
013 JNC+02
223 JNC-60d GTO LOOP1
261 C=
060 A/C
305 C=
060 sqrt(AB)
10E A=C ALL
138 C=L
0AE A<>C ALL
128 L=C
2BE C=-C
01D C=
061 A+C
138 C=L
013 JNC+02
233 JNC-58d GTO LOOP2
269 C=
060 AB/C
05E C=|C|
10E A=C ALL
04E C
35C
2BE
090
21C =
250
250
050 -2 10^(-9)
01D C=
060 A+C
2FE ?C<0
383 JNC-16d GTO LOOP2
138 C=L C = R_F(a,b,b)
10E A=C ALL
04E C
35C =
0D0 3
135 C=
060 A*C
1F8 C=O
13D C=
060 AB*C
278 C=Q
025 C=
060 AB+C
268 Q=C
1F8 C=O
10E A=C ALL
01D C=
060 A+C
1E8 O=C
0F8 C=X
10E A=C ALL
0B8 C=Y
013 JNC+02
253 JNC-54d GTO LOOP1
01D C=
060 A+C
078 C=Z
10C ?FSET8
013 JNC+02
0B8 C=Y
025 C=
060 AB+C
238 C=P
025 C=
060 AB+C
238 C=P
025 C=
060 AB+C
04E C
35C =
150 5
269 C=
060 AB/C
128 L=C
305 C=
060 sqrt(AB)
138 C=L
13D C=
060 AB*C
10E A=C
1F8 C=O
0AE A<>C ALL
261 C=
060 A/C
278 C=Q
025 C=
060 AB+C
10E A=C ALL
046 C
270 =
038 T
0AE A<>C ALL
028 T=C
128 L=C
2BE C=-C
01D C=
061 A+C
138 C=L
269 C=
060 AB/C
05E C=|C|
10E A=C ALL
04E C
35C
2BE
090
21C =
250
250
050 -2 10^(-9)
01D C=
060 A+C
2FE ?C<0
223 JNC-60d GTO LOOP1
24C ?FSET9
01B JNC+03
244 CLRF9
203 JNC-64d GTO LOOP1
3B5 ?NCXQ
050 R^
345 ?NCGO
042 CLA

( 315 words )

RJ STACK	INPUTS	OUTPUTS
T	p>0	/
Z	z	/
Y	y	/
X	x	R_J(x;y;z;p)

Example:

     7 ENTER^
     4 ENTER^
     2 ENTER^
     1 XEQ "RJ" yields   R_J(1,2,4,7) = 0.1478544450 ( in 14 seconds )

RJZ STACK	INPUTS	OUTPUTS
T	p>0	/
Z	z	/
Y	y	/
X	x	R_J(x,y+iz,y-iz,p)

Example:

     7 ENTER^
     4 ENTER^
     2 ENTER^
     1 XEQ "RJZ" >>>> R_J ( 1 , 2+4i , 2-4i , 7 ) = 0.1260390651 ( in 16 seconds )

Notes:

-We have saved 20 words and about 1 second.
-Make similar modifications ( adding 3CC   ?KEY 360    ?C RTN and modifying the JNC- ... ) if you don't have a "newest" HP41.

References:

[1] B.C. Carlson - "Numerical Computation of Real or Complex Elliptic Integrals"
[2} W. J. Nellis & B.C. Carlson - "Reduction and Evaluation of Elliptic Integrals"
[3] D. G. Zill & B. C. Carlson - "Symmetric Elliptic Integrals of the Third Kind"
[4] https://en.wikipedia.org/w/index.php?title=Carlson_symmetric_form&oldid=1037291384"
[5] B. C. Carlson - "A Table of Elliptic Integrals of the Third Kind"

hp41programs

Carlson Elliptic Integrals for the HP-41