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^b  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 =  §₀^v ( 1 - m sin² 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 sin² 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  a₀ = 1 ; b₀ = (1-m)^1/2 ; c₀ =  m^1/2
  we determine       a₁ , b₁ , c₁ ; ....... ; a_n , b_n , c_n   with   a_k= (a_k-1+ b_k-1)/2  ;  b_k= (a_k-1.b_k-1)^1/2  ;  c_k = (a_k-1- b_k-1)/2
  and we stop when  c_n = 0    ( in fact, when a_n-1= a_n )
 

-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   v_n = 2ⁿa_nu*180/pi    ( in degrees )
 and then,  v_n-1, ... , v₀  from the recurrence relation:    sin (2v_k-1-v_k) = ( c_k/a_k ) sin v_k
 then,  sn ( u | m ) = sin v₀  ;   cn ( u | m ) = cos v₀   and   dn ( u | m ) = ( 1 - m sin² v₀ )^1/2

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

      sn ( u | m ) = sn ( u m^1/2 | 1/m ) / m^1/2
      cn ( u | m ) = dn ( u m^1/2 | 1/m )
      dn ( u | m ) = cn ( u m^1/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 = c₁/a₁  , ..... ,   Rnn = c_n/a_n
 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 ) =  §₀^pi/2 ( 1 - m sin² t )^-1/2 dt
 - the complete elliptic integrals of the 2nd kind:         E ( m ) =  §₀^pi/2  ( 1 - m sin² t )^1/2 dt
 - the incomplete elliptic integrals of the 1st kind:   F ( v | m ) =  §₀^v ( 1 - m sin² t )^-1/2 dt       ( 0 < = v < = 90° )
 - the incomplete elliptic integrals of the 2nd kind:  E ( v | m ) =  §₀^v  ( 1 - m sin² 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/2a_n
    E(m) = K(m) - K(m).(c₀²+2¹c₁²+ ... +2ⁿc_n²)/2

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

    tan ( v_n+1 - v_n ) = ( b_n/a_n ) tan v_n  and  v₀ = 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 ( v_n+1 - 2v_n ) = (( b_n/a_n ) tan v_n - tan v_n) / ( 1 + ( b_n/a_n ) tan² v_n )

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

   tan ( v_n+1 - 2v_n ) = (( b_n/a_n )  - 1 ) / ( 1 / tan v_n + ( b_n/a_n ) tan v_n )         ( without forgetting lines 45-48-50 )

Then,    F ( v | m ) = v_n / ( 2ⁿa_n )
and       E ( v | m ) = Z ( v | m ) + ( E / K ) F ( v | m )

where    Z( v | m ) = c₁ sin v₁ + ... + c_n sin v_n   is the Jacobian Zeta function.

Notes:

 1-R00 = c₀²+2¹c₁²+ ... +2ⁿc_n²        R01 = 1 , 2 , 4 , 8 , .....      R02 = a_k
    R03 = b_k and  F ( v | m )               R04 = v_k                            R05 = Z( v | m )          R06 = a_n-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:

       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

-This program uses the following property:

      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

-This transformation is performed until  x , y , z are close enough to apply  R_F(x;y;z) = µ ^-1/2   with  µ = (x+y+z)/3     ( we have  R_F(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"  >>>>  R_F(2;3;4) = 0.584082842  ( in 10 seconds )
 

  Degenerate Form
 

-Carson also defines   R_C(x;y) = R_F(x;y;y)  if y > 0
  and  R_C(x;y) = (x/(x-y))^1/2 R_C(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"  >>>>  R_C(1;3) = 0.675510859   ( 11 seconds )

        -3  ENTER^
         1     R/S       >>>>   R_C(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   R_D(x;y;z) = R_J(x;y;z;z)    ( cf c) for R_J )

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"  >>>>  R_D(2;3;4) = 0.165105273    ( 30 seconds )

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

  R_G(x;y;z) =  (1/4)  §₀^infinity  ( ( t + x ).( t + y ).( t + z ) ) ^-1/2 .( x/(t+x) + y/(t+y) + z/(t+z) )  t.dt

And we have:  2.R_G(x;y;z) = z.R_F(x;y;z) - (x-z)(y-z)/3 R_D(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"  >>>>  R_G(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

   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

-We have  R_J(x,x,x,x) = x ^-3/2
-The following program applies  R_J(x;y;z;p) = 3 Sum_{n=0,1,2,....,k}  R_F(a_n,b_n,b_n)/4ⁿ + 1/(4^k+1µ ^3/2)

  where   x₀ = x , y₀ = y , z₀ = z , p₀ = p        ;      x_n+1 = ( x_n+ L_n )/4 , y_n+1 = ( y_n+ L_n )/4 , z_n+1 = ( z_n + L_n )/4 , p_n+1 = ( p_n +L_n )/4

    with    L_n = x_n^1/2y_n^1/2 + x_n^1/2z_n^1/2 + y_n^1/2z_n^1/2
               a_n = ( p_n( x_n^1/2 + y_n^1/2 + z_n^1/2 ) + ( x_ny_nz_n )^1/2 )²    ;   b_n = p_n ( p_n + L_n )²

    and    µ = (x_k+1+y_k+1+z_k+1+2p_k+1)/5

-The iterations stop when  x_n , y_n , z_n , p_n  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"  >>>>  R_J(1;2;3;4) = 0.239848100   ( in 42 seconds )
 

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

Note:   The Cauchy principal value of  R_J(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,  R_J(1;2;3;p) = 0  for p = -0.775227.......

  and this program gives  R_J(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  R_F ( 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"  >>>>  R_F ( 2 , 3+4i , 3-4i ) = 0.543421944    ( in 19 seconds )

             d-2) Elliptic Integrals of the third kind
 

-This program computes  R_J ( 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"  >>>>  R_J ( 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 ) = §₀^phi ( 1 - m sin² t )^-1/2 dt = sin (phi) . R_F ( cos²(phi) ; 1 - m.sin²(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 | k² )   add  X^2 after line 09 and place k in Y-register ( instead of m )
  but m may be negative wheras k² 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 ) = §₀^phi ( 1 - m sin² t )^1/2 dt
                                      = sin (phi) . R_F ( cos²(phi) ; 1 - m.sin²(phi) ; 1 ) - (m/3) sin³(phi) . R_D ( cos²(phi) ; 1 - m.sin²(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 | k² )   add  X^2 after line 10 and place k in Y-register ( instead of m )
  but m may be negative wheras k² 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 ) = §₀^phi ( 1 + n sin² t ) ^-1 ( 1 - m sin² t ) ^-1/2 dt
                                            = sin (phi) . R_F ( cos²(phi) ; 1 - m.sin²(phi) ; 1 ) - (n/3) sin³(phi) . R_J ( cos²(phi) ; 1 - m.sin²(phi) ; 1 ; 1 + n.sin²(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 | k² )   add  X^2 after line 14 and place k in Y-register ( instead of m )
  but m may be negative wheras k² may not.
-You must have  -90° < phi <= 90°
-This program works in all angular mode.
-If  1 + n.sin²(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"
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall