The Museum of HP Calculators

Jacobian Elliptic Functions and Elliptic Integrals for the HP-41

This program is by Jean-Marc Baillard and is used here by permission.

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°) Jacobian Elliptic Functions

2°) Legendre Elliptic Integrals of the first and second kind

3°) Carlson Elliptic Integrals

a) Carlson Elliptic Integrals of the first kind
b) Carlson Elliptic Integrals of the second kind ( degenerate form & symmetric form )
c) Carlson Elliptic Integrals of the third kind
d) Complex arguments
d-1) Elliptic Integrals of the first kind
d-2) Elliptic Integrals of the third kind

4°) Legendre Elliptic Integrals    ( another method via Carlson Integrals )

a) Legendre Elliptic Integrals of the first kind
b) Legendre Elliptic Integrals of the second kind
c) Legendre Elliptic Integrals of the third kind
d) Legendre Elliptic Integrals ( all three )

NB:  Integrals are symbolized by  "§" :      §a f(x) dx  ...

1°) Jacobian Elliptic Functions

-There are 12 Jacobian elliptic functions.They can be defined like this:

Let  m be a number that verifies:  0 < = m < = 1      ( m is called the parameter )
if    u =  §0v ( 1 - m sin2 t ) -1/2 dt    ( 0 < = v < = 90° )
the angle  v  = am u    is called the amplitude
and the 3 Jacobian elliptic functions  sn ; cn ; dn  are defined by:

sn ( u | m ) = sin v  ; cn ( u | m) = cos v ;  dn ( u | m ) = ( 1 - m sin2 v )1/2

-The 9 other Jacobian elliptic functions can be obtained with the relations:

cd = cn / dn  ;  dc = dn / cn  ;  ns = 1 / sn
sd = sn / dn  ;   nc = 1 / cn   ;  ds = dn / sn
nd = 1 / dn   ;   sc = sn / cn  ;  cs = cn / sn

The Arithmetic-Geometric Mean:

-The following program uses the process of the Arithmetic-Geometric Mean (AGM):

Starting with  a0 = 1 ; b0 = (1-m)1/2 ; c0 =  m1/2
we determine       a1 , b1 , c1 ; ....... ; an , bn , cn   with   ak= (ak-1+ bk-1)/2  ;  bk= (ak-1.bk-1)1/2  ;  ck = (ak-1- bk-1)/2
and we stop when  cn = 0    ( in fact, when an-1= an )

-This program calculates only the 3 copolar trio:  sn ( u | m ) ; cn ( u | m ) ; dn ( u | m ).
From these, all the 9 other elliptic functions can be determined by the relations mentioned above.

-When the AGM process stops, the HP-41 calculates   vn = 2nanu*180/pi    ( in degrees )
and then,  vn-1, ... , v0  from the recurrence relation:    sin (2vk-1-vk) = ( ck/ak ) sin vk
then,  sn ( u | m ) = sin v0  ;   cn ( u | m ) = cos v0   and   dn ( u | m ) = ( 1 - m sin2 v0 )1/2

-If m > 1, the program uses the relations:

sn ( u | m ) = sn ( u m1/2 | 1/m ) / m1/2
cn ( u | m ) = dn ( u m1/2 | 1/m )
dn ( u | m ) = cn ( u m1/2 | 1/m )

-If m = 1:      sn ( u | m ) = tanh u ; cn ( u | m ) = dn ( u | m ) = sech u

-If m < 0:

sn ( u | m ) = sd ( u (1-m)1/2 | -m/(1-m) ) / (1-m)1/2
cn ( u | m ) = cd ( u (1-m)1/2 | -m/(1-m) )
dn ( u | m ) = nd ( u (1-m)1/2 | -m/(1-m) )

Notes:

1-R00 = 1 , 2 , ... , n  and then   n-1, ... , 1 , 0   ;    R01 = c1/a1  , ..... ,   Rnn = cn/an
2-It seems that n is never greater than 7.
Therefore, synthetic registers M and N can be replaced by R08 and R09
3-If  m < -9999999999 the program can give wrong results!
4-Flags F01 and F02 are cleared by this program.

01  LBL "JEF"
02  DEG
03  CF 01
04  CF 02
05  X<>Y
06  X#0?
07  X>0?
08  GTO 02
09  CHS
10  STO Z
11  1
12  +
13  ST/ Z
14  SQRT
15  *
16  X<>Y
17  SF 01
18  LBL 02
19  1
20  X>Y?
21  GTO 04
22  X#Y?
23  GTO 03
24  X<> Z
25  ENTER^
26  E^X
27  LASTX
28  CHS
29  E^X
30  +
31  2
32  X<>Y
33  /
34  X<>Y
35  ST+ X
36  E^X-1
37  RCL X
38  2
39  +
40  /
41  GTO 06
42  LBL 03
43  RDN
44  STO Z
45  SQRT
46  *
47  X<>Y
48  1/X
49  1
50  SF 02
51  LBL 04
52  STO 00
53  X<> Z
54  STO M
55  RDN
56  STO N
57  -
58  SQRT
59  SIGN
60  LASTX
61  LBL 01
62  RCL Y
63  RCL Y
64  -
65  STO IND 00
66  X<> T
67  RCL Y
68  +
69  ST/ IND 00
70  2
71  /
72  R^
73  X=Y?
74  GTO 02
75  RCL Z
76  *
77  SQRT
78  ISG 00
79  CLX
80  GTO 01
81  LBL 02
82  2
83  RCL 00
84  Y^X
85  *
86  RCL M
87  *
88  R-D
89  LBL 10
90  ENTER^
91  SIN
92  RCL IND 00
93  *
94  ASIN
95  +
96  2
97  /
98  DSE 00
99  GTO 10
100  COS
101  LAST X
102  SIN
103  ENTER^
104  X^2
105  RCL N
106  CHS
107  *
108  1
109  +
110  SQRT
111  STO T
112  X<> M
113  SIGN
114  X<> N
115  RDN
116  FS?C 01
117  GTO 05
118  FC?C 02
119  GTO 06
120  X<>Y
121  X<> T
122  SQRT
123  *
124  GTO 06
125  LBL 05
126  R^
127  R^
128  ST/ T
129  ST/ Z
130  1/X
131  RDN
132  SIGN
133  ST- L
134  X<> L
135  CHS
136  SQRT
137  *
138  LBL 06
139  CLA
140  END

( 193 bytes )

 STACK INPUTS OUTPUTS Z / dn ( u | m ) Y m cn ( u | m ) X u sn ( u | m )

If  0 < =  m < 1 ,  m is saved in T and u is saved in L

Examples:

1- Evaluate  sn ( 0.7 | 0.3 )     cn ( 0.7 | 0.3 )     dn ( 0.7 | 0.3 )

0.3 ENTER^
0.7 XEQ "JEF"    gives      sn ( 0.7 | 0.3 )  = 0.632304777  ( in 11 seconds )
RDN        cn ( 0.7 | 0.3 ) = 0.774719736
RDN        dn ( 0.7 | 0.3 ) = 0.938113640

2- In the same way                 sn ( 0.7 | 1 ) = 0.604367777          sn ( 0.7 | 2 ) = 0.564297008        sn ( 0.7 | -3 ) = 0.759113422
cn ( 0.7 | 1 ) = 0.796705460          cn ( 0.7 | 2 ) = 0.825571855        cn ( 0.7 | -3 ) = 0.650958381
dn ( 0.7 | 1 ) = 0.796705460          dn ( 0.7 | 2 ) = 0.602609139        dn ( 0.7 | -3 ) =1.651895747

2°) Legendre Elliptic Integrals of the first and second kind

-Legendre elliptic integrals are:

- the complete elliptic integrals of the 1st kind:          K ( m ) =  §0pi/2 ( 1 - m sin2 t )-1/2 dt
- the complete elliptic integrals of the 2nd kind:         E ( m ) =  §0pi/2  ( 1 - m sin2 t )1/2 dt
- the incomplete elliptic integrals of the 1st kind:   F ( v | m ) =  §0v ( 1 - m sin2 t )-1/2 dt       ( 0 < = v < = 90° )
- the incomplete elliptic integrals of the 2nd kind:  E ( v | m ) =  §0v  ( 1 - m sin2 t )1/2 dt

-The elliptic integrals of the third kind are much more complicated and will be calculated differently  ( cf  below )

-Here we assume:  0 < = m < = 1
-The same AGM scheme is used here, after which:

K(m) = pi/2an
E(m) = K(m) - K(m).(c02+21c12+ ... +2ncn2)/2

-To obtain   F ( v | m) , the HP-41 must determine n angles  v1 , ... , vn
In Abramowitz' "Handbook of Mathematical Functions" , we find:

tan ( vn+1 - vn ) = ( bn/an ) tan vn  and  v0 = v

-However, with this recurrence relation, we can't determine the angles v  in the proper quadrant
and so, we obtain wrong results. That's why I 've changed this formula into:

tan ( vn+1 - 2vn ) = (( bn/an ) tan vn - tan vn) / ( 1 + ( bn/an ) tan2 vn )

-But there is still a problem: if one angle is equal to 90°, tan2 v  will produce "OUT OF RANGE"
-Finally, the HP-41 uses the relation:

tan ( vn+1 - 2vn ) = (( bn/an )  - 1 ) / ( 1 / tan vn + ( bn/an ) tan vn )         ( without forgetting lines 45-48-50 )

Then,    F ( v | m ) = vn / ( 2nan )
and       E ( v | m ) = Z ( v | m ) + ( E / K ) F ( v | m )

where    Z( v | m ) = c1 sin v1 + ... + cn sin vn   is the Jacobian Zeta function.

Notes:

1-R00 = c02+21c12+ ... +2ncn2        R01 = 1 , 2 , 4 , 8 , .....      R02 = ak
R03 = bk and  F ( v | m )               R04 = vk                            R05 = Z( v | m )          R06 = an-1
2-K( 1 ) is infinite  ( 9.999999999 E99 in the T-register )                  E( 1 ) = 1
F ( v | 1 ) = ln ( tan ( pi/4 + v/2 ) )                                             E ( v | 1 ) = sin v
3-K( m ) and E( m ) are automatically calculated in the T and Z registers
4-The angle v must be expressed in degrees..

01  LBL "ELI"
02  DEG
03  STO 04
04  X<>Y
05  1
06  X#Y?
07  GTO 00
08  90
09  TAN
10  X<>Y
11  R^
12  LASTX
13  X#Y?
14  GTO 02
15  TAN
16  X<>Y
17  SIGN
18  RTN
19  LBL 02
20  RDN
21  SIN
22  ST+ Y
23  X<>Y
24  LASTX
25  COS
26  /
27  LN
28  1
29  X<> Z
30  RTN
31  LBL 00
32  STO 01
33  STO 02
34  X<>Y
35  STO 00
36  -
37  SQRT
38  STO 03
39  CLX
40  STO 05
41  LBL 01
42  1
43  STO 06
44  RCL 04
45  ENTER^
46  ST+ Y
47  TAN
48  RCL 03
49  RCL 02
50  /
51  ST- 06
52  X<>Y
53  ST* Y
54  X#0?
55  1/X
56  +
57  X=0?
58  STO 06
59  X#0?
60  ST/ 06
61  X<> 06
62  ATAN
63  -
64  STO 04
65  RCL 02
66  STO 06
67  RCL 03
68  -
69  2
70  ST* 01
71  /
72  ENTER^
73  X^2
74  RCL 01
75  *
76  ST+ 00
77  RCL 02
78  RCL 03
79  *
80  SQRT
81  X<> 03
82  RCL 02
83  +
84  2
85  /
86  STO 02
87  RCL 04
88  SIN
89  R^
90  *
91  ST+ 05
92  RCL 02
93  RCL 06
94  X#Y?
95  GTO 01
96  RCL 04
97  RCL 02
98  /
99  RCL 01
100  /
101  D-R
102  STO 03
103  PI
104  RCL 02
105  ST+ X
106  /
107  ENTER^
108  ST* 00
109  RCL 00
110  2
111  /
112  -
113  STO Z
114  X<>Y
115  /
116  RCL 03
117  *
118  RCL 05
119  +
120  RCL 03
121  X<>Y
122  END

( 148 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS T / K ( m ) Z / E ( m ) Y m F ( v | m) X v ( degrees ) E ( v | m) L / Z ( v | m)

Z ( v | m)   is in the L-register only if m < 1

Examples:

1-If   v = 84° and m = 0.7

0.7 ENTER^
84 XEQ "ELI"  yields       E ( 84° | 0.7 ) = 1.184070048        ( in 16 seconds )
RDN       F ( 84° | 0.7 ) = 1.884976271
RDN               E ( 0.7 ) = 1.241670567
RDN               K ( 0.7 ) = 2.075363134
and        LASTX        Z ( 84° | 0.7 ) = 0.056306180

2-If   v = 84°  and m = 1

1 ENTER^
84 XEQ "ELI"  yields     E ( 84° | 1 ) = 0.994521895           ( in 2 seconds )
RDN      F ( 84° | 1 ) = 2.948700239
RDN              E ( 1 ) = 1
RDN              K ( 1 ) = 9.999999999E99    ( in fact infinity )

Complete Elliptic Integrals of the first and second kind:

-If you are only interested by the complete elliptic integrals,
here is a shorter and faster program to compute K ( m ) and E ( m ):

01  LBL "CEI"
02  STO 00
03  SIGN
04  ENTER^
05  STO 01
06  LASTX
07  -
08  SQRT
09  ENTER^
10  LBL 01
11  CLX
12  RCL Z
13  RCL Y
14  -
15  2
16  ST* 01
17  /
18  X^2
19  RCL 01
20  *
21  ST+ 00
22  RDN
23  STO Z
24  STO T
25  X<>Y
26  ST* Z
27  +
28  2
29  /
30  X<>Y
31  SQRT
32  R^
33  X#Y?
34  GTO 01
35  ST+ X
36  PI
37  X<>Y
38  /
39  ENTER^
40  ST* 00
41  RCL 00
42  2
43  /
44  -
45  END

( 63 bytes / SIZE 002 )

 STACK INPUTS OUTPUTS T / K (m) Z / K (m) Y / K (m) X m* E (m)

* 0 < = m < 1

Example:

E ( 0.7 ) = 1.241670567   and
K ( 0.7 ) = 2.075363134   are obtained in 5 seconds.

3°) Carlson Elliptic Integrals

a) Carlson Elliptic Integrals of the first kind

-Carlson has given a new definition of a standard elliptic integral of the first kind:

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

-This program uses the following property:

RF(x;y;z) = RF((x+L)/4;(y+L)/4;(z+L)/4)  where  L = x1/2y1/2 + x1/2z1/2 + y1/2z1/2

-This transformation is performed until  x , y , z are close enough to apply  RF(x;y;z) = µ -1/2   with  µ = (x+y+z)/3     ( we have  RF(x;x;x) = x -1/2 )
-Actually, the iterations may be stopped earlier. Then, the function could be evaluated by a Taylor series.
-But this approach would required more bytes.

Data Registers:   R00 unused , R01 thru R03: temp
Flags: /
Subroutines: /

01  LBL "RF"
02  X<Y?
03  X<>Y
04  RDN
05  X>Y?
06  X<>Y
07  STO 01
08  X<> T
09  X>Y?
10  X<>Y
11  STO 02
12  X<>Y
13  STO 03
14  LBL 01
15  RCL 01
16  SQRT
17  RCL 02
18  SQRT
19  STO Z
20  RCL 03
21  SQRT
22  ST* T
23  +
24  *
25  +
26  ST+ 01
27  ST+ 02
28  ST+ 03
29  4
30  ST/ 01
31  ST/ 02
32  ST/ 03
33  RCL 03
34  RCL 01
35  ST- Y
36  RCL 02
37  RCL 03
38  +
39  +
40  /
41   E-5
42  X<Y?
43  GTO 01
44  3
45  LASTX
46  /
47  SQRT
48  END

( 68 bytes / SIZE 004 )

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

Example:      4  ENTER^
3  ENTER^
2  XEQ "RF"  >>>>  RF(2;3;4) = 0.584082842  ( in 10 seconds )

Degenerate Form

-Carson also defines   RC(x;y) = RF(x;y;y)  if y > 0
and  RC(x;y) = (x/(x-y))1/2 RC(x-y;-y)      if y < 0

Data Registers:   R00  thru R03: temp
Flags: /
Subroutine: "RF"

01  LBL "RC"
02  1
03  STO 00
04  X<> Z
05  X>0?
06  GTO 00
07  X<>Y
08  STO 00
09  X<>Y
10  -
11  ST/ 00
12  LASTX
13  CHS
14  LBL 00
15  STO Z
16  X<>Y
17  XEQ "RF"
18  RCL 00
19  SQRT
20  *
21  END

( 35 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Y y / X x RC(x;y)

Example:      3  ENTER^
1  XEQ "RC"  >>>>  RC(1;3) = 0.675510859   ( 11 seconds )

-3  ENTER^
1     R/S       >>>>   RC(1;-3) = 0.274653072   ( 10  seconds )

b) Carlson Elliptic Integrals of the second kind

-They are only a degenerate form of the integrals of the third kind   RD(x;y;z) = RJ(x;y;z;z)    ( cf c) for RJ )

01  LBL "RD"
02  0
03  +
04  XEQ "RJ"
05  END

( 15 bytes / SIZE 013 )

 STACK INPUTS OUTPUTS Z z / Y y / X x RD(x;y;z)

Example:      4  ENTER^
3  ENTER^
2  XEQ "RD"  >>>>  RD(2;3;4) = 0.165105273    ( 30 seconds )

-However, Carlson has also defined a  symmetric Elliptic Integral of the second kind:

RG(x;y;z) =  (1/4)  §0infinity  ( ( t + x ).( t + y ).( t + z ) ) -1/2 .( x/(t+x) + y/(t+y) + z/(t+z) )  t.dt

And we have:  2.RG(x;y;z) = z.RF(x;y;z) - (x-z)(y-z)/3 RD(x;y;z) + ( x.y/z )1/2

Data Registers:   R00  thru R13: temp
Flags: /
Subroutines: "RF" & "RD"

01  LBL "RG"
02  STO 04
03  RDN
04  STO 05
05  X<>Y
06  STO 06
07  R^
08  XEQ "RF"
09  RCL 04
10  SQRT
11  RCL 05
12  SQRT
13  *
14  RCL 06
15  ST* Z
16  SQRT
17  /
18  +
19  STO 00
20  RCL 04
21  RCL 06
22  ST- Y
23  RCL 05
24  -
25  *
26  STO 13
27  RCL 06
28  RCL 05
29  RCL 04
30  XEQ "RD"        or replace this line by   XEQ "RJ"   and add   RCL 06  after line 26
31  RCL 13
32  *
33  3
34  /
35  RCL 00
36  +
37  2
38  /
39  END

( 54 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS Z z / Y y / X x RG(x;y;z)

Example:      4  ENTER^
3  ENTER^
2  XEQ "RG"  >>>>  RG(2;3;4) = 1.725503029  ( in 42 seconds )

-If  one of the arguments is zero, do not place it in  Z-register ( there would be a DATA ERROR line 17 )
Or add   X#0?   X<>Y  after line 03

c) Carlson Elliptic Integrals of the third kind

-The elliptic integral of the third kind is defined by

RJ(x;y;z;p) =  (3/2)  §0infinity  ( 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

-We have  RJ(x,x,x,x) = x -3/2
-The following program applies  RJ(x;y;z;p) = 3 Sumn=0,1,2,....,k  RF(an,bn,bn)/4n + 1/(4k+1µ 3/2)

where   x0 = x , y0 = y , z0 = z , p0 = p        ;      xn+1 = ( xn+ Ln )/4 , yn+1 = ( yn+ Ln )/4 , zn+1 = ( zn + Ln )/4 , pn+1 = ( pn +Ln )/4

with    Ln = xn1/2yn1/2 + xn1/2zn1/2 + yn1/2zn1/2
an = ( pn( xn1/2 + yn1/2 + zn1/2 ) + ( xnynzn )1/2 )2    ;   bn = pn ( pn + Ln )2

and    µ = (xk+1+yk+1+zk+1+2pk+1)/5

-The iterations stop when  xn , yn , zn , pn  are close enough to produce a good accuracy.
-This program also computes the Cauchy principal value of the integral if  p < 0

Data Registers:   R00  thru R12: temp
Flags: /
Subroutines: "RF" & "RC"

01  LBL "RJ'"
02  X<Y?
03  X<>Y
04  RDN
05  X>Y?
06  X<>Y
07  STO 04
08  X<> T
09  X<Y?
10  X<>Y
11  STO 06
12  X<>Y
13  STO 05
14  ST- Y
15  R^
16  -
17  *
18  RCL 05
19  R^
20  STO 07
21  1
22  STO 08
23  STO 09
24  STO 10
25  CLX
26  STO 11
27  STO 12
28  RDN
29  X>0?
30  GTO 01
31  -
32  STO 09
33  /
34  STO 10
35  RCL 05
36  +
37  ENTER^
38  X<> 07
39  *
40  RCL 04
41  RCL 06
42  *
43  RCL 05
44  ST/ Z
45  /
46  XEQ "RC"
47  STO 11
48  RCL 06
49  RCL 05
50  RCL 04
51  XEQ "RF"
52  ST- 11
53  LBL 01
54  RCL 04
55  SQRT
56  ENTER^
57  STO 01
58  RCL 05
59  SQRT
60  ST* Z
61  STO 02
62  RCL 06
63  SQRT
64  ST* T
65  ST* 02
66  +
67  ST* 01
68  +
69  RCL 07
70  *
71  +
72  X^2
73  RCL 07
74  RCL 01
75  RCL 02
76  +
77  ST+ 04
78  ST+ 05
79  ST+ 06
80  RCL 07
81  +
82  STO 07
83  X^2
84  *
85  ENTER^
86  XEQ "RF"
87  RCL 08
88  *
89  ST+ 12
90  4
91  ST/ 04
92  ST/ 05
93  ST/ 06
94  ST/ 07
95  ST/ 08
96  RCL 07
97  RCL 06
98  -
99  ABS
100  RCL 04
101  ST- Y
102  RCL 05
103  +
104  RCL 06
105  ST+ Z
106  +
107  RCL 07
108  ST+ X
109  +
110  /
111   E-5
112  X<Y?
113  GTO 01
114  RCL 08
115  LASTX
116  5
117  /
118  ENTER^
119  SQRT
120  *
121  /
122  RCL 12
123  3
124  *
125  +
126  RCL 10
127  *
128  RCL 11
129  3
130  *
131  +
132  RCL 09
133  /
134  END

( 175 bytes / SIZE 013 )

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

Examples:     4  ENTER^
3  ENTER^
2  ENTER^
1  XEQ "RJ"  >>>>  RJ(1;2;3;4) = 0.239848100   ( in 42 seconds )

-4  ENTER^
3  ENTER^
2  ENTER^
1     R/S        >>>>  RJ(1;2;3;-4) =  -0.237867695  ( in 58 seconds )

Note:   The Cauchy principal value of  RJ(x;y;z;p)  has a zero for some p < 0  and a loss of significant digits is to be expected near the zero:

-For instance,  RJ(1;2;3;p) = 0  for p = -0.775227.......

and this program gives  RJ(1;2;3;-0.775227) = 8.4498 10-8

Line131 adds  0.1220082998  and  -0.1220080653   Therefore,  the result has probably no more than 2 or 3 significant figures.

d) Complex Arguments

d-1) Elliptic Integrals of the first kind

-This program computes  RF ( x , y+i.z , y-i.z )
-This is not the general case, but it's useful when the radicand has one real root and a pair of complex conjugate roots.

Data Registers:   R00 unused  R01 thru R03: temp
Flags: /
Subroutines: /

01  LBL "RFZ"
02  STO 01
03  RDN
04  STO 02
05  X<>Y
06  ABS
07  STO 03
08  LBL 01
09  RCL 03
10  RCL 02
11  R-P
12  STO Z
13  SQRT
14  ST+ X
15  X<>Y
16  2
17  /
18  COS
19  *
20  RCL 01
21  SQRT
22  *
23  +
24  ST+ 01
25  ST+ 02
26  4
27  ST/ 01
28  ST/ 02
29  ST/ 03
30  RCL 03
31  RCL 02
32  RCL 01
33  -
34  ABS
35  +
36  RCL 01
37  RCL 02
38  ST+ X
39  +
40  /
41   E-5
42  X<Y?
43  GTO 01
44  3
45  LASTX
46  /
47  SQRT
48  END

( 67 bytes / SIZE 004 )

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

Example:      4  ENTER^
3  ENTER^
2  XEQ "RFZ"  >>>>  RF ( 2 , 3+4i , 3-4i ) = 0.543421944    ( in 19 seconds )

d-2) Elliptic Integrals of the third kind

-This program computes  RJ ( x , y+i.z , y-i.z , p )  with  p > 0

Data Registers:   R00 unused  R01 thru R09: temp  ( R09 may be replaced by R00 )
Flags: /
Subroutine:  "RF"

01  LBL "RJZ"
02  STO 04
03  RDN
04  STO 05
05  RDN
06  ABS
07  STO 06
08  X<>Y
09  STO 07
10  CLX
11  STO 08
12  SIGN
13  STO 09
14  LBL 01
15  RCL 06
16  RCL 05
17  R-P
18  STO Z
19  SQRT
20  ST+ X
21  X<>Y
22  2
23  /
24  COS
25  *
26  STO 01
27  RCL 04
28  SQRT
29  ST* T
30  ST+ 01
31  *
32  +
33  ST+ 04
34  ST+ 05
35  X<> 07
36  ST+ 07
37  ENTER^
38  X<> 01
39  *
40  +
41  X^2
42  RCL 01
43  RCL 07
44  X^2
45  *
46  ENTER^
47  XEQ "RF"
48  RCL 09
49  *
50  ST+ 08
51  4
52  ST/ 04
53  ST/ 05
54  ST/ 06
55  ST/ 07
56  ST/ 09
57  RCL 05
58  RCL 07
59  -
60  ABS
61  RCL 06
62  +
63  RCL 04
64  RCL 05
65  -
66  ABS
67  +
68  RCL 04
69  RCL 05
70  RCL 07
71  +
72  ST+ X
73  +
74  /
75   E-5
76  X<Y?
77  GTO 01
78  RCL 09
79  LASTX
80  5
81  /
82  ENTER^
83  SQRT
84  *
85  /
86  RCL 08
87  3
88  *
89  +
90  END

( 120 bytes / SIZE 010 )

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

Example:      4  ENTER^
3  ENTER^
2  ENTER^
1  XEQ "RJZ"  >>>>  RJ ( 1 , 2+3i , 2-3i , 4 ) = 0.205644141   ( in 52 seconds )

4°) Legendre Elliptic Integrals  ( another method via Carlson Integrals )

a) Legendre Elliptic Integrals of the first kind

-We use the formula:    F ( phi | m ) = §0phi ( 1 - m sin2 t )-1/2 dt = sin (phi) . RF ( cos2(phi) ; 1 - m.sin2(phi) ; 1 )

Data Registers:   R00  thru R03: temp
Flags: /
Subroutine: "RF"

01  LBL "LEI1"
02  1
03  P-R
04  X^2
05  X<>Y
06  STO 00
07  X^2
08  1
09  R^
10  -
11  *
12  +
13  X=Y?
14  X#0?
15  GTO 00
16  DEG
17  90
18  TAN
19  RTN
20  LBL 00
21  1
22  XEQ "RF"
23  RCL 00
24  *
25  END

( 39 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Y m / X phi F(phi | m)

Example:     in DEG mode:

0.7  ENTER^
84   XEQ "LEI1"  >>>>  F( 84° | 0.7 ) =  1.884976270    ( in 13 seconds )

-If you prefer  F( phi , k ) = F ( phi | k2 )   add  X^2 after line 09 and place k in Y-register ( instead of m )
but m may be negative wheras k2 may not.
-You must have  -90° < phi <= 90°  If  90° < phi <= 180°  you'll get  F ( 180°- phi | m )
-This program works in all angular mode.

b) Legendre Elliptic Integrals of the second kind

Formula:    E ( phi | m ) = §0phi ( 1 - m sin2 t )1/2 dt
= sin (phi) . RF ( cos2(phi) ; 1 - m.sin2(phi) ; 1 ) - (m/3) sin3(phi) . RD ( cos2(phi) ; 1 - m.sin2(phi) ; 1 )

Data Registers:   R00  thru R14: temp
Flags: /
Subroutines: "RF" & "RD"

01  LBL "LEI2"
02  1
03  P-R
04  X^2
05  STO 04
06  X<>Y
07  STO 00
08  X^2
09  1
10  R^
11  STO 13
12  -
13  *
14  +
15  X=Y?
16  X#0?
17  GTO 00
18  SIGN
19  RTN
20  LBL 00
21  STO 05
22  1
23  XEQ "RF"
24  STO 14
25  1
26  RCL 04
27  RCL 05
28  XEQ "RD"        or replace this line by   0   +   XEQ "RJ"
29  RCL 00
30  X^2
31  RCL 13
32  *
33  *
34  RCL 14
35  X<>Y
36  3
37  /
38  -
39  RCL 00
40  *
41  END

( 57 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS Y m / X phi E (phi | m)

Example:     in DEG mode:

0.7  ENTER^
84   XEQ "LEI2"  >>>>  E ( 84° | 0.7 ) =  1.184070047    ( in 53 seconds )

-If you prefer  E( phi , k ) = E ( phi | k2 )   add  X^2 after line 10 and place k in Y-register ( instead of m )
but m may be negative wheras k2 may not.
-You must have  -90° < phi <= 90°
-This program works in all angular mode.

c) Legendre Elliptic Integrals of the third kind

Formula:     ¶ ( n ; phi | m ) = §0phi ( 1 + n sin2 t ) -1 ( 1 - m sin2 t ) -1/2 dt
= sin (phi) . RF ( cos2(phi) ; 1 - m.sin2(phi) ; 1 ) - (n/3) sin3(phi) . RJ ( cos2(phi) ; 1 - m.sin2(phi) ; 1 ; 1 + n.sin2(phi) )

-This sign convention for n is opposite that of Abramowitz and Steigun

Data Registers:   R00  thru R15: temp
Flags: /
Subroutines: "RF" & "RJ"

01  LBL "LEI3"
02  1
03  P-R
04  X^2
05  STO 04
06  X<>Y
07  STO 15
08  X^2
09  STO 06
10  R^
11  STO 13
12  CLX
13  SIGN
14  R^
15  -
16  *
17  +
18  X=Y?
19  X#0?
20  GTO 00
21  DEG
22  90
23  TAN
24  RTN
25  LBL 00
26  STO 05
27  1
28  XEQ "RF"
29  STO 14
30  RCL 04
31  RCL 06
32  1
33  RCL 13
34  +
35  *
36  +
37  1
38  RCL 05
39  R^
40  XEQ "RJ"
41  RCL 15
42  X^2
43  *
44  RCL 13
45  *
46  RCL 14
47  X<>Y
48  3
49  /
50  -
51  RCL 15
52  *
53  END

( 70 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS Z n / Y m / X phi ¶( n ; phi | m )

Example:     in DEG mode:

0.9  ENTER^
0.7  ENTER^
84   XEQ "LEI3"  >>>>  ¶ ( 0.9 ; 84° | 0.7 ) =  1.336853616    ( in 69 seconds )

-If you prefer  ¶ ( phi , n , k ) = ¶ ( n ; phi | k2 )   add  X^2 after line 14 and place k in Y-register ( instead of m )
but m may be negative wheras k2 may not.
-You must have  -90° < phi <= 90°
-This program works in all angular mode.
-If  1 + n.sin2(phi) is non-negative, register R15 may be replaced by register R00

d) Legendre Elliptic Integrals ( all three )

-Considering their great similarity, the last 3 routines may form a single program, thus saving many bytes:

Data Registers:   R00  thru R15: temp
Flags:   F01    Set flag F01 to compute Legendre elliptic integrals of the first kind
F02     ------- F02 ------------------------------------------- second --          Set only one of these 3 flags
F03     ------- F03 ------------------------------------------- third ----
Subroutines: "RF" & "RJ"

01  LBL "LEI"
02  1
03  STO 06
04  R^
05  STO 13
06  RDN
07  P-R
08  X^2
09  STO 04
10  X<>Y
11  STO 15
12  X^2
13  FS? 03
14  STO 06
15  1
16  R^                 add   X^2   after this line if you prefer  k  instead of  m
17  FC? 03
18  STO 13
19  -
20  *
21  +
22  X=Y?
23  X#0?
24  GTO 00
25  SIGN
26  FS? 02
27  RTN
28  DEG
29  90
30  TAN
31  RTN
32  LBL 00
33  STO 05
34  1
35  XEQ "RF"
36  FS? 01
37  GTO 00
38  STO 14
39  RCL 04
40  RCL 06
41  1
42  RCL 13
43  +
44  *
45  +
46  FC? 03
47  RCL 06
48  1
49  RCL 05
50  RCL 04
51  XEQ "RJ"
52  RCL 15
53  X^2
54  *
55  RCL 13
56  *
57  RCL 14
58  X<>Y
59  3
60  /
61  -
62  LBL 00
63  RCL 15
64  *
65  END

( 87 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS Z n* / Y m / X phi Leg-Ell-Int

*  n is used for the integrals of the 3rd kind only

Leg-Ell-Int  =  E ( phi | m )        if flag F01 is set
With     Leg-Ell-Int  =  F ( phi | m )        if flag F02 is set
Leg-Ell-Int  =  ¶ ( n ; phi | m )    if flag F03 is set

References:   Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
B.C. Carlson, "Numerical Computation of Real or Complex Elliptic Integrals"