hp41programs

Elliptic Integrals for the HP-41
 

Overview
 

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

  a) Carlson Elliptic Integrals of the first kind

    a-1) Program#1
    a-2) Program#2  ( faster )

  b) Carlson Elliptic Integrals of the second kind ( degenerate form & symmetric form )
  c) Carlson Elliptic Integrals of the third kind

    c-1) Program#1
    c-2) Program#2  ( faster )
    c-3) Program#3  ( even faster )

  d) Complex arguments

    d-1) Elliptic Integrals of the first kind

      d-1a)  Program#1
      d-1b)  Program#2  ( faster )

    d-2) Elliptic Integrals of the third kind

      d-2a)  Program#1
      d-2b)  Program#2  ( faster )
      d-2c)  Program#3  ( even faster )

  e) Micro-Code Functions

     e-1) Integrals of the first kind - real & complex
     e-2) Integrals of the third kind - real & complex

 3°) Legendre Elliptic Integrals    ( 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°) 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 arithmetic-geometric mean ( AGM ) scheme is used here ( cf "Jacobian Elliptic Functions" ), 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.
 

2°) Carlson Elliptic Integrals
 

       a) Carlson Elliptic Integrals of the first kind
 

            a-1) Program#1
 

-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 requires more bytes ( cf next paragraph ).
 

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)

 
Examples:

      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 )
 

            a-2) Program#2
 

-Now, "RF2" calculates  R_F(x;y;z)  by the same method as "RF" but we stop the iterations when  Max { | X | , | Y | , | Z | } is smaller than 1/40
             where  X = 1 - x/µ  ,  Y = 1 - y/µ  ,   Z = 1 - z/µ

-Then,   R_F(x;y;z) =  ( 1 - E/10 + F/14 + E²/24 - 3 E.F/44 ) µ ^-1/2    with  E = XY - Z²   and  F = X.Y.Z
 

Data Registers:   R00 thru R05: temp
Flags: /
Subroutines: /
 
 

01  LBL "RF2" 02  STO 01 03  RDN 04  STO 02 05  X<>Y 06  STO 03 07  40 08  1/X 09  STO 04 10  LBL 01 11  RCL 01        12  SQRT 13  RCL 02 14  SQRT 15  STO Z 16  RCL 03 17  SQRT 18  ST* T 19  +

20  * 21  + 22  ST+ 01 23  ST+ 02 24  ST+ 03 25  4 26  ST/ 01 27  ST/ 02 28  ST/ 03 29  RCL 01        30  RCL 02 31  RCL 03 32  + 33  + 34  3 35  / 36  STO 00 37  RCL 03 38  -

39  RCL 00 40  / 41  STO 05 42  ABS 43  RCL 04 44  X<Y? 45  GTO 01 46  RCL 00        47  RCL 02 48  - 49  RCL 00 50  / 51  ST* 05 52  ABS 53  RCL 04 54  X<Y? 55  GTO 01 56  RCL 00 57  RCL 01

58  - 59  RCL 00 60  / 61  ABS 62  RCL 04 63  X<Y? 64  GTO 01 65  RCL 05        66  LASTX 67  ST* 05 68  X^2 69  - 70  STO 04 71  24 72  / 73  .1 74  - 75  RCL 05 76  3

77  * 78  44 79  / 80  - 81  RCL 04 82  * 83  RCL 05        84  14 85  / 86  + 87  1 88  + 89  RCL 00 90  SQRT 91  / 92  END

 
   ( 118 bytes / SIZE 006 )
 
 

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

 
Example:

    4  ENTER^
    2  ENTER^
    1  XEQ "RF2"  >>>>  R_F(1,2,4) = 0.685085817  ( in 6 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 R04: temp
Flags: /
Subroutines:  The M-code functions RF & RJ  or "RF" & "RJ" ...  In this case, use other registers to avoid a register usage conflict.
 
 

01  LBL "RG" 02  STO 01 03  RDN 04  STO 02 05  X<>Y 06  STO 03        07  R^ 08  RF

09  RCL 01 10  SQRT 11  RCL 02        12  SQRT 13  * 14  RCL 03 15  ST* Z 16  SQRT

17  / 18  + 19  STO 00        20  RCL 03 21  RCL 03 22  RCL 02 23  RCL 01 24  RJ

25  RCL 03 26  RCL 01        27  - 28  * 29  RCL 02 30  RCL 03 31  - 32  *

33  3 34  / 35  RCL 00        36  + 37  2 38  / 39  END

 
   ( 49 bytes / SIZE 004 )
 
 

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.725503028  ( in 15 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
 

            c-1) Program#1
 

-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.
 

            c-2) Program#2
 

-This new program computes  R_J(x;y;z;p)  for p>0  only.
-Instead of calling "RF" as a subroutine, it calculates R_F(a,b,b)  directly by 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

-So, it is much faster.
 

Data Registers:   R00 thru R08: temp
Flags: /
Subroutines: /

-Line 117 is a three-byte GTO 01
 
 

01  LBL "RJ2"   02  STO 01   03  RDN   04  STO 02   05  RDN   06  STO 03   07  X<>Y   08  STO 04   09  CLX   10  STO 00          11  STO 08   12  SIGN   13  STO 07   14  RAD   15  LBL 01   16  RCL 01   17  SQRT   18  ENTER^   19  STO 05   20  RCL 02   21  SQRT   22  ST* Z   23  STO 06   24  RCL 03

25  SQRT   26  ST* T   27  ST* 06   28  +   29  ST* 05   30  +   31  RCL 04          32  *   33  +   34  X^2   35  RCL 04   36  RCL 05   37  RCL 06   38  +   39  ST+ 01   40  ST+ 02   41  ST+ 03   42  RCL 04   43  +   44  STO 04   45  X^2   46  *   47  X=Y?   48  GTO 02

49  -   50  X<0?   51  GTO 03   52  STO 05          53  X<>Y   54  /   55  SQRT   56  ENTER^   57  ST+ Y   58  SIGN   59  LASTX   60  -   61  /   62  LN1+X   63  2   64  /   65  GTO 04   66  LBL 02   67  SQRT   68  1/X   69  GTO 05   70  LBL 03   71  CHS   72  STO 05

73  X<>Y   74  /   75  SQRT   76  ATAN   77  LBL 04   78  RCL 05          79  SQRT   80  /   81  LBL 05   82  RCL 07   83  *   84  ST+ 00   85  4   86  ST/ 01   87  ST/ 02   88  ST/ 03   89  ST/ 04   90  ST/ 07   91  RCL 08   92  RCL 07   93  RCL 01   94  RCL 02   95  +   96  RCL 03

97  +   98  RCL 04          99  ST+ X 100  + 101  5 102  / 103  ENTER^ 104  SQRT 105  * 106  / 107  RCL 00 108  3 109  * 110  + 111  STO 08 112  ST- Y 113  / 114  ABS 115  2 E-9 116  X<Y? 117  GTO 01 118  RCL 08 119  DEG 120  END

 
    ( 153 bytes / SIZE 009 )
 
 

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

 
Example:

    4  ENTER^
    3  ENTER^
    2  ENTER^
    1  XEQ "RJ2"  >>>>  R_J(1;2;3;4) = 0.239848100   ( in 22 seconds - instead of 42 with "RJ" )

-The iterations stop when 2 consecutive results are nearly equal.
 

            c-3) Program#3
 

-Let    X = 1 - x/µ  ,  Y = 1 - y/µ  ,  Z = 1 - z/µ  ,  P = 1 - p/µ

-"RJ3" stops the iterations when   Max { | X | , | Y | , | Z | , | P | } is smaller than 1/67

-Then,  R_J(x;y;z;p) = 3 Sum_{n=0,1,2,....,k}  R_F(a_n,b_n,b_n)/4ⁿ + (1+eps)/(4^k+1µ ^3/2)

  where  eps = -3.E/14 + F/6 + 9.E²/88 - 3.G/22 - 9.E.F/52 + 3.H/26

  with        E = X.(Y+Z) + Y.Z - 3.P²          G = ( 2.X.Y.Z + E.P + 3.P³ ).P
  and         F = X.Y.Z + 2.E.P + 4.P³           H = X.Y.Z.P²
 

Data Registers:   R00 thru R11: temp
Flags: /
Subroutines: /

-Lines 112-122-132-142 are three-byte GTO 01
 
 

01  LBL "RJ3"   02  STO 01   03  RDN   04  STO 02   05  RDN   06  STO 03   07  X<>Y   08  STO 04   09  CLX   10  STO 00   11  SIGN   12  STO 07   13  67   14  1/X   15  STO 08           16  RAD   17  LBL 01   18  RCL 01   19  SQRT   20  ENTER^   21  STO 05   22  RCL 02   23  SQRT   24  ST* Z   25  STO 06   26  RCL 03   27  SQRT   28  ST* T   29  ST* 06   30  +   31  ST* 05   32  +   33  RCL 04   34  *   35  +   36  X^2   37  RCL 04

38  RCL 05   39  RCL 06   40  +   41  ST+ 01   42  ST+ 02   43  ST+ 03   44  RCL 04   45  +   46  STO 04   47  X^2   48  *   49  X=Y?   50  GTO 02   51  -   52  X<0?   53  GTO 03   54  STO 05           55  X<>Y   56  /   57  SQRT   58  ENTER^   59  ST+ Y   60  SIGN   61  LASTX   62  -   63  /   64  LN1+X   65  2   66  /   67  GTO 04   68  LBL 02   69  SQRT   70  1/X   71  GTO 05   72  LBL 03   73  CHS   74  STO 05

75  X<>Y   76  /   77  SQRT   78  ATAN   79  LBL 04   80  RCL 05   81  SQRT   82  /   83  LBL 05   84  RCL 07   85  *   86  ST+ 00   87  4   88  ST/ 01   89  ST/ 02   90  ST/ 03   91  ST/ 04   92  ST/ 07   93  RCL 01           94  RCL 02   95  +   96  RCL 03   97  +   98  RCL 04   99  ST+ X 100  + 101  5 102  / 103  STO 05 104  RCL 01 105  - 106  RCL 05 107  / 108  STO 06 109  ABS 110  RCL 08 111  X<Y?

112  GTO 01 113  RCL 05 114  RCL 02 115  - 116  RCL 05 117  / 118  STO 09 119  ABS 120  RCL 08 121  X<Y? 122  GTO 01 123  RCL 05 124  RCL 03 125  - 126  RCL 05 127  / 128  STO 10 129  ABS 130  RCL 08         131  X<Y? 132  GTO 01 133  RCL 05 134  RCL 04 135  - 136  RCL 05 137  / 138  STO 11 139  ABS 140  RCL 08 141  X<Y? 142  GTO 01 143  DEG 144  RCL 09 145  RCL 10 146  + 147  X<> 06 148  ST* 06

149  RCL 09 150  LASTX 151  * 152  ST+ 06 153  * 154  STO 09 155  RCL 11 156  X^2 157  STO 10 158  3 159  * 160  ST- 06 161  RCL 10 162  ST+ X 163  RCL 06 164  + 165  RCL 11         166  * 167  ST+ X 168  RCL 09 169  + 170  STO 08 171  LASTX 172  + 173  RCL 06 174  RCL 10 175  ST* 09 176  + 177  RCL 11 178  * 179  - 180  ST* 11 181  RCL 06 182  88 183  / 184  RCL 08 185  52

186  / 187  - 188  3 189  * 190  14 191  1/X 192  - 193  RCL 06 194  * 195  RCL 09 196  26 197  / 198  + 199  RCL 11         200  22 201  / 202  - 203  3 204  * 205  RCL 08 206  6 207  / 208  + 209  RCL 07 210  ST* Y 211  + 212  RCL 05 213  ENTER^ 214  SQRT 215  * 216  / 217  RCL 00 218  3 219  * 220  + 221  END

 
  ( 271 bytes / SIZE 012 )
 
 

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

 
Example:

     4  ENTER^
     3  ENTER^
     2  ENTER^
     1  XEQ "RJ3"  >>>>  R_J(1,2,3,4) = 0.239848100  ( in 14 seconds )

  R_J(1,2,4,7) = 0.147854445  is obtained in 17 seconds.

-The execution time increases when x , y , z , p  have huge ratios.
 

       d) Complex Arguments
 

            d-1) Elliptic Integrals of the first kind
 

         d-1a) Program#1
 
 

-This program computes  R_F ( x , y+i.z , y-i.z )
-This is not most 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 Y 13  LASTX 14  + 15  ST+ X 16  SQRT 17  RCL 01 18  SQRT

19  * 20  + 21  ST+ 01 22  ST+ 02 23  4 24  ST/ 01 25  ST/ 02 26  ST/ 03 27  RCL 03

28  RCL 02 29  RCL 01 30  - 31  ABS 32  + 33  RCL 01        34  RCL 02 35  ST+ X 36  +

37  / 38   E-5 39  X<Y? 40  GTO 01        41  3 42  LASTX 43  / 44  SQRT 45  END

 
    ( 64 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 14 seconds )
 

         d-1b) Program#2
 

-This routine uses the same series expansion as "RF2"
 

Data Registers:   R00 thru R05: temp
Flags: /
Subroutines: /
 
 

01  LBL "RFZ2" 02  STO 01  03  RDN 04  STO 02        05  X<>Y 06  STO 03 07  40 08  1/X 09  STO 04 10  LBL 01 11  RCL 03 12  RCL 02 13  R-P 14  STO Y 15  LASTX 16  + 17  ST+ X

18  SQRT 19  RCL 01        20  SQRT 21  * 22  + 23  ST+ 01 24  ST+ 02 25  4 26  ST/ 01 27  ST/ 02 28  ST/ 03 29  RCL 01 30  RCL 02 31  ST+ X 32  + 33  3 34  /

35  STO 00 36  RCL 01        37  - 38  RCL 00 39  / 40  STO 05 41  ABS 42  RCL 04 43  X<Y? 44  GTO 01 45  RCL 00 46  RCL 02 47  - 48  RCL 03 49  R-P 50  RCL 00 51  /

52  RCL 04 53  X<Y? 54  GTO 01 55  RCL 05        56  X^2 57  RCL Z 58  X^2 59  ST* 05 60  X<>Y 61  - 62  STO 04 63  24 64  / 65  .1 66  - 67  RCL 05 68  3

69  * 70  44 71  / 72  - 73  RCL 04        74  * 75  RCL 05 76  14 77  / 78  + 79  1 80  + 81  RCL 00 82  SQRT 83  / 84  END

 
   ( 110 bytes / SIZE 006 )
 
 

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

 
Example:

     4  ENTER^
     2  ENTER^
     1  XEQ "RFZ2"  >>>>  R_F ( 1 , 2+4i , 2-4i ) = 0.631488738    ( in 11 seconds )
 

            d-2) Elliptic Integrals of the third kind
 

-The 3 following programs compute  R_J ( x , y+i.z , y-i.z , p )  with  p > 0
 

         d-2a) Program#1
 
 

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  ENTER^

19  STO Z 20  LASTX 21  + 22  ST+ X 23  SQRT 24  STO 01  25  RCL 04        26  SQRT 27  ST* T 28  ST+ 01 29  * 30  + 31  ST+ 04 32  ST+ 05 33  X<> 07 34  ST+ 07 35  ENTER^ 36  X<> 01

37  * 38  + 39  X^2 40  RCL 01 41  RCL 07        42  X^2 43  * 44  ENTER^ 45  XEQ "RF" 46  RCL 09 47  * 48  ST+ 08 49  4 50  ST/ 04 51  ST/ 05 52  ST/ 06 53  ST/ 07 54  ST/ 09

55  RCL 05 56  RCL 07 57  - 58  ABS 59  RCL 06        60  + 61  RCL 04  62  RCL 05 63  - 64  ABS 65  + 66  RCL 04 67  RCL 05 68  RCL 07 69  + 70  ST+ X 71  + 72  /

73   E-5 74  X<Y? 75  GTO 01 76  RCL 09        77  LASTX 78  5 79  / 80  ENTER^ 81  SQRT 82  * 83  / 84  RCL 08 85  3 86  * 87  + 88  END

 
    ( 118 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 46 seconds )

         d-2b) Program#2
 

-Like "RJ2" this routine calculates R_F(a,b,b)  directly, without calling "RF"
 

Data Registers:   R00 thru R07: temp
Flag:   F10
Subroutines:  /

-Lines 115 & 117 are three-byte GTO 01
 
 

01  LBL "RJZ2"   02  STO 01    03  RDN   04  STO 02   05  RDN   06  STO 03   07  X<>Y   08  STO 04   09  CLX   10  STO 00   11  STO 06          12  SIGN   13  STO 07   14  RAD   15  SF 10   16  LBL 01   17  RCL 03   18  RCL 02   19  R-P   20  ENTER^   21  STO Z   22  LASTX   23  +   24  ST+ X

25  SQRT   26  STO 05    27  RCL 01   28  SQRT   29  ST* T   30  ST+ 05   31  *   32  +   33  ST+ 01   34  ST+ 02   35  X<> 04   36  ST+ 04   37  ENTER^   38  X<> 05   39  *   40  +   41  X^2   42  RCL 05          43  RCL 04   44  ENTER^   45  *   46  *   47  X=Y?   48  GTO 02

49  -   50  X<0?   51  GTO 03   52  STO 05          53  X<>Y   54  /   55  SQRT   56  ENTER^   57  ST+ Y   58  SIGN   59  LASTX   60  -   61  /   62  LN1+X   63  2   64  /   65  GTO 04   66  LBL 02   67  SQRT   68  1/X   69  GTO 05   70  LBL 03   71  CHS   72  STO 05

73  X<>Y   74  /   75  SQRT   76  ATAN   77  LBL 04   78  RCL 05          79  SQRT   80  /   81  LBL 05   82  RCL 07   83  *   84  ST+ 00   85  4   86  ST/ 01   87  ST/ 02   88  ST/ 03   89  ST/ 04   90  ST/ 07   91  RCL 06   92  RCL 07   93  RCL 02   94  RCL 04   95  +   96  ST+ X

97  RCL 01   98  +   99  5 100  / 101  ENTER^ 102  SQRT 103  * 104  / 105  RCL 00        106  3 107  * 108  + 109  STO 06 110  ST- Y 111  / 112  ABS 113  2 E-9 114  X<Y? 115  GTO 01 116  FS?C 10 117  GTO 01 118  RCL 06 119  DEG 120  END

 
  ( 160 bytes / SIZE 008 )
 
 

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 "RJZ2"  >>>>  R_J ( 1 , 2+3i , 2-3i , 4 ) = 0.205644141   ( in 26 seconds - instead of 46 seconds with "RJZ" )

-It may happen that 2 consecutive approximations are nearly equal too early,
 so we use a flag to prevent a premature stop!
 

         d-3b) Program#3
 

-"RJZ3" uses the same series expansion and the same termination criterion as "RJ3"
 

Data Registers:   R00 thru R11: temp
Flags: /
Subroutines: /

-Lines 109-121-131 are three-byte GTO 01
 
 

01  LBL "RJZ3"   02  STO 01    03  RDN   04  STO 02   05  RDN   06  STO 03   07  X<>Y   08  STO 04   09  CLX   10  STO 00          11  SIGN   12  STO 07   13  67   14  1/X   15  STO 08   16  RAD   17  LBL 01   18  RCL 03   19  RCL 02   20  R-P   21  ENTER^   22  STO Z   23  LASTX   24  +   25  ST+ X   26  SQRT   27  STO 05   28  RCL 01   29  SQRT   30  ST* T   31  ST+ 05   32  *   33  +   34  ST+ 01   35  ST+ 02   36  X<> 04

37  ST+ 04   38  ENTER^   39  X<> 05    40  *   41  +   42  X^2   43  RCL 05   44  RCL 04          45  ENTER^   46  *   47  *   48  X=Y?   49  GTO 02   50  -   51  X<0?   52  GTO 03   53  STO 05   54  X<>Y   55  /   56  SQRT   57  ENTER^   58  ST+ Y   59  SIGN   60  LASTX   61  -   62  /   63  LN1+X   64  2   65  /   66  GTO 04   67  LBL 02   68  SQRT   69  1/X   70  GTO 05   71  LBL 03   72  CHS

73  STO 05    74  X<>Y   75  /   76  SQRT   77  ATAN   78  LBL 04   79  RCL 05          80  SQRT   81  /   82  LBL 05   83  RCL 07   84  *   85  ST+ 00   86  4   87  ST/ 01   88  ST/ 02   89  ST/ 03   90  ST/ 04   91  ST/ 07   92  RCL 01   93  RCL 02   94  RCL 04   95  +   96  ST+ X   97  +   98  5   99  / 100  STO 05 101  RCL 01 102  - 103  RCL 05 104  / 105  STO 06 106  ABS 107  RCL 08 108  X<Y?

109  GTO 01 110  RCL 05  111  RCL 02 112  - 113  RCL 03 114  R-P 115  RCL 05        116  / 117  STO 09 118  ABS 119  RCL 08 120  X<Y? 121  GTO 01 122  RCL 05 123  RCL 04 124  - 125  RCL 05 126  / 127  STO 11 128  ABS 129  RCL 08 130  X<Y? 131  GTO 01 132  DEG 133  RCL 05 134  RCL 02 135  - 136  ST+ X 137  RCL 05 138  / 139  X<> 06 140  ST* 06 141  RCL 09 142  X^2 143  ST+ 06 144  *

145  STO 09 146  RCL 11  147  X^2 148  STO 10        149  3 150  * 151  ST- 06 152  RCL 10 153  ST+ X 154  RCL 06 155  + 156  RCL 11 157  * 158  ST+ X 159  RCL 09 160  + 161  STO 08 162  LASTX 163  + 164  RCL 06 165  RCL 10 166  ST* 09 167  + 168  RCL 11 169  * 170  - 171  ST* 11 172  RCL 06 173  88 174  / 175  RCL 08 176  52 177  / 178  - 179  3 180  *

181  14 182  1/X 183  - 184  RCL 06        185  * 186  RCL 09 187  26 188  / 189  + 190  RCL 11 191  22 192  / 193  - 194  3 195  * 196  RCL 08 197  6 198  / 199  + 200  RCL 07 201  ST* Y 202  + 203  RCL 05 204  ENTER^ 205  SQRT 206  * 207  / 208  RCL 00 209  3 210  * 211  + 212  END

 
   ( 264 bytes / SIZE 012 )
 
 

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

 
Example:

      7  ENTER^
      4  ENTER^
      2  ENTER^
      1  XEQ "RJZ3"  >>>>  R_J ( 1 , 2+4i , 2-4i , 7 ) = 0.126039065   ( in 18 seconds )
 

       e) Micro-Code Functions
 

            e-1) Integrals of the first kind - real & complex
 

Code  Mnemonics            Comments

086    "F"
012    "R"
104    CLRF8                  The first executable word of RF
244    CLRF9
00E    A=0  ALL
0F8    READ3(X)
2EE    ?C#0?  ALL
04F    JC+09
166    A=A+1  S&X
03B    JNC+07
09A   "Z"
006    "F"
012    "R"
108    SETF8                  The first executable word of RFZ
248    SETF9
00E    A=0  ALL            We first check that no more than 1 argument  ( among  x , y , z  for RF or  x , x+iy , x-iy  for RFZ )  equals zero
0B8    READ2(Y)
2EE    ?C#0?  ALL
017    JC+02
166    A=A+1  S&X
078    READ1(Z)
2EE    ?C#0?  ALL
017    JC+02
166    A=A+1  S&X
130    LDI S&X
002    002
306    ?A<C  S&X
0B5   ?NCGO               otherwise we get DATA ERROR
0A2    282D
0F8    READ3(X)          these 8 words check for alpha data in registers X , Y , Z and execute SETDEC
361    ?NCXQ              If you don't want to check for alpha data, replace them ( 0F8  361  050  0B8  10E  078  355  050  )
050    14D8                   by  2A0    SETDEC
0B8    READ2(Y)
10E    A=C  ALL
078    READ1(Z)
355    ?NCXQ
050    14D5
04E    C=0  ALL
028    WRIT0(T)
0F8    READ3(X)      loop       This loop begins at the address F1A7  in my ROM , you'll have to change the ?NCGOs written in red ( near the end )
2F9    ?NCXQ            C=                                                                                       if you store this word at another address in your own ROM
060    18BE                 sqrt(C)
128    WRIT4(L)
0B8   READ2(Y)
10C   ?FSET8
097    JC+18d
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
070    N=C  ALL
078    READ1(Z)
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
268    WRIT9(Q)
10E    A=C  ALL
0B0    C=N  ALL
135    ?NCXQ           C=
060    184D               A*C
0F0    C<>N  ALL
10E    A=C  ALL
278    READ9(Q)
01D    ?NCXQ          C=
060     1807               A+C
0CB    JNC+25d
10E    A=C  ALL
135    ?NCXQ           C=
060    184D               A*C
070    N=C  ALL
078    READ1(Z)
10E    A=C  ALL
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
0B0    C=N  ALL
01D    ?NCXQ          C=
060     1807               A+C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
070    N=C  ALL
10E    A=C  ALL
0B8   READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
01D    ?NCXQ          C=
060     1807               A+C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
10E    A=C  ALL
138    READ4(L)
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
0B0    C=N  ALL
01D    ?NCXQ          C=
060     1807               A+C
070    N=C  ALL
10E    A=C  ALL
0B8   READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
04E    C=0  ALL        C
35C    PT=12             =
110    LD@PT-4        4
128    WRIT4(L)
261    ?NCXQ           C=
060    1898                A/C
0A8    WRIT2(Y)
0B0    C=N  ALL
10E    A=C  ALL
078    READ1(Z)
10C    ?FSET8
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
138    READ4(L)
261    ?NCXQ           C=
060    1898                A/C
068    WRIT1(Z)
0B0    C=N  ALL
10E    A=C  ALL
0F8    READ3(X)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
138    READ4(L)
261    ?NCXQ           C=
060    1898                A/C
0E8    WRIT3(X)
10E    A=C  ALL
0B8   READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
078    READ1(Z)
10C    ?FSET8
013    JNC+02
0B8   READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
046    C=0 S&X        C
270    RAMSLCT      =
038    READDATA   T
2BE    C=-C-1 MS    C=-C
0AE    A<>C
028    WRIT0(T)
070    N=C  ALL
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
0B0    C=N  ALL
261    ?NCXQ           C=
060    1898                A/C
05E    C=0 MS          C= | C |
10E    A=C  ALL
04E    C=0                 |
2BE    C=-C-1 MS    |
35C    PT=12            |
090    LD@PT-2       |      C = -2 E-9
21C    PT=2              |
250    LD@PT-9      |
250    LD@PT-9      |
050    LD@PT-1      |
01D    ?NCXQ          C=
060     1807               A+C
2FE    ?C#0 MS
29D    ?NCGO         if not carry goto
3C6    F1A7             loop
24C    ?FSET9
023    JNC+04
244    CLRF9
29D    ?NCGO         if not carry goto
3C6    F1A7             loop
04E    C=0               C
35C    PT=12           =
0D0    LD@PT-3     3
10E    A=C  ALL
270    RAMSLCT      C=
038    READDATA    T
261    ?NCXQ           C=
060    1898                A/C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
0E8    WRIT3(X)
3E0    RTN
 
 

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

 
Example:

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

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

 
Example:

    4  ENTER^
    2  ENTER^
    1  XEQ "RFZ"  gives  R_F ( 1 , 2+4i , 2-4i ) = 0.631488738    ( in 5.5 seconds )

-The execution time may be longer if the arguments have huge ratios.
-For instance:  R_F(0,10 ^-70,10⁸⁰) =  1.740801762 10 ^-38   is obtained in 9 seconds

-There will be a DATA ERROR if  x , y , z < 0 ( with RF )  or  x < 0  ( with RFZ )
-Space can be saved if you don't check for zero-arguments...
-In this case, there will be an infinite loop if you execute RF with x = y = 0 and z # 0  or if you execute RFZ with y = z = 0
  press  ENTER^  ON to stop the loop with a "newest" HP-41
-On older HP-41s, insert

3CC   ?KEY
360    ?C RTN    at the beginning of the loop: any infinite loop will be stopped by pressing any key for 1 second.

-RF & RFZ  do not use any data register and the alpha registers are undisturbed
-However, RF uses the "scratch" register Q

-The routine stops when 2 consecutive approximations are nearly equal ( CPU-flag9 is also used by RFZ to avoid a premature stop ).
-You can certainly write a faster function if you use the method of "RF2" & "RFZ2":
  the number of iterations will be reduced, but the series expansion needs several extra words.
 
 

            e-2) Integrals of the third kind - real & complex
 
 

Code  Mnemonics            Comments

08A    "J"
012    "R"
104    CLRF8                  The first executable word of RJ
244    CLRF9
00E    A=0  ALL
0F8    READ3(X)
2EE    ?C#0?  ALL
04F    JC+09
166    A=A+1  S&X
03B    JNC+07
09A   "Z"
00A    "J"
012    "R"
108    SETF8                  The first executable word of RJZ
248    SETF9
00E    A=0  ALL            We check that no more than 1 argument  ( among  x , y , z  for RJ or  x , x+iy , x-iy  for RJZ )  equals zero
0B8    READ2(Y)
2EE    ?C#0?  ALL
017    JC+02
166    A=A+1  S&X
078    READ1(Z)
2EE    ?C#0?  ALL
017    JC+02
166    A=A+1  S&X
130    LDI S&X
002    002
306    ?A<C  S&X
0B5   ?NCGO               otherwise we get DATA ERROR
0A2    282D
0F8    READ3(X)          these 12 words check for alpha data in registers X , Y , Z , T  execute SETDEC and return T in C
10E    A=C  ALL          If you don't want to check for alpha data, replace them ( 0F8  10E  0B8  355  050  046  270  038  10E  078  355  050  )
0B8    READ2(Y)         by  2A0    SETDEC
355    ?NCXQ                    046    C=0 S&X
050    14D5                        270    RAMSLCT
046    C=0 S&X                 038    READDATA
270    RAMSLCT
038    READDATA
10E    A=C  ALL
078    READ1(Z)
355    ?NCXQ
050    14D5
228    WRIT8(P)
04E    C=0
028    WRIT0(T)
268    WRIT9(Q)
35C    PT=12             C=
050    LD@PT-1        1
1E8    WRIT7(O)
0F8    READ3(X)      loop       This loop begins at the address FAD4  in my ROM , you'll have to change the ?NCGOs written in red ( near the end )
2F9    ?NCXQ           C=                                                                                               if you store this word at another address in your own ROM
060    18BE                sqrt(C)
168    WRIT5(M)
0B8    READ2(Y)
10C    ?FSET8
10F    JC+33d
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
1A8    WRIT6(N)
078    READ1(Z)
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
128    WRIT4(L)
10E    A=C  ALL
1B8    READ6(N)
01D    ?NCXQ          C=
060     1807               A+C
070    N=C  ALL
10E    A=C  ALL
178    READ5(M)
135    ?NCXQ           C=
060    184D               A*C
0F0    C<>N  ALL
10E    A=C  ALL
178    READ5(M)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
1B8    READ6(N)
0AE    A<>C  ALL
1A8    WRIT6(N)
138    READ4(L)
135    ?NCXQ           C=
060    184D               A*C
128    WRIT4(L)
10E    A=C  ALL
0B0    C=N  ALL
133    JNC+38d
10E    A=C  ALL
135    ?NCXQ           C=
060    184D               A*C
070    N=C  ALL
078    READ1(Z)
10E    A=C  ALL
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
0B0    C=N  ALL
01D    ?NCXQ          C=
060     1807               A+C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
128    WRIT4(L)
10E    A=C  ALL
0B8    READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
01D    ?NCXQ          C=
060     1807               A+C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
070    N=C  ALL
10E    A=C  ALL
178    READ5(M)
01D    ?NCXQ          C=
060     1807               A+C
1A8    WRIT6(N)
178    READ5(M)
10E    A=C  ALL
0B0    C=N  ALL
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
138    READ4(L)
01D    ?NCXQ          C=
060     1807               A+C
070    N=C  ALL
138    READ4(L)
10E    A=C  ALL
178    READ5(M)
135    ?NCXQ           C=
060    184D               A*C
168    WRIT5(M)
0B0    C=N  ALL
10E    A=C  ALL
0F8    READ3(X)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
04E    C=0  ALL        C
35C    PT=12             =
110    LD@PT-4        4
128    WRIT4(L)
261    ?NCXQ           C=
060    1898                A/C
0E8    WRIT3(X)
0B0    C=N  ALL
10E    A=C  ALL
0B8    READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
138    READ4(L)
261    ?NCXQ           C=
060    1898                A/C
0A8    WRIT2(Y)
0B0    C=N  ALL
10E    A=C  ALL
078    READ1(Z)
10C    ?FSET8
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
138    READ4(L)
261    ?NCXQ           C=
060    1898                A/C
068    WRIT1(Z)
1B8    READ6(N)
10E    A=C  ALL
238    READ8(P)
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
178    READ5(M)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
135    ?NCXQ           C=
060    184D               A*C
168    WRIT5(M)
0B0    C=N  ALL
10E    A=C  ALL
238    READ8(P)
01D    ?NCXQ          C=
060     1807               A+C
070    N=C  ALL
10E    A=C  ALL
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
238    READ8(P)
135    ?NCXQ           C=
060    184D               A*C
1A8    WRIT6(N)
0B0    C=N  ALL
10E    A=C  ALL
138    READ4(L)
261    ?NCXQ           C=
060    1898                A/C
228    WRIT8(P)
1B8    READ6(N)
2BE    C=-C-1 MS    C=-C
10E    A=C  ALL
178    READ5(M)
01D    ?NCXQ          C=
060     1807               A+C
1A8    WRIT6(N)       then the HP41 calculates R_F(a,b,b) according to the sign of  a-b
2EE    ?C#0  ALL
1F3    JNC+62d
2FE    ?C#0  MS
09B    JNC+19d
05E    C=0 MS          C= | C |
10E    A=C  ALL
1A8    WRIT6(N)
178    READ5(M)
261    ?NCXQ           C=
060    1898                A/C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
070    N=C  ALL
130    LDI S&X         The 8 following words calculate ATAN in radians without disturbing the current angular mode
090    090
358    ST=C XP
144    CLRF6
284    CLRF7
0B0    C=N  ALL
205    ?NCXQ           C=
040    1081                Atan(C)
103    JNC+32d
10E    A=C  ALL
178    READ5(M)
261    ?NCXQ           C=
060    1898                A/C
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
168    WRIT5(M)
2BE    C=-C-1 MS    C=-C
00E    A=0  ALL        A
35C    PT=12             =
162    A=A+1@PT    1
01D    ?NCXQ          C=
060     1807               A+C
070    N=C  ALL
178    READ5(M)
10E    A=C  ALL
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
0B0    C=N  ALL
261    ?NCXQ           C=
060    1898                A/C
084    CLRF5            C
1CD   ?NCXQ          =
06C    1B73              LN1+C
10E    A=C  ALL
04E    C=0  ALL        C
35C    PT=12             =
090    LD@PT-2        2
261    ?NCXQ           C=
060    1898                A/C               now, we have  atanh ((a-b)/a)^(1/2)
070    N=C  ALL
1B8    READ6(N)
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
10E    A=C  ALL
0B0    C=N  ALL
0AE   A<>C  ALL
261    ?NCXQ           C=
060    1898                A/C
033    JNC+06
178    READ5(M)
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
22D   ?NCXQ           C=
060    188B                1/C
10E    A=C  ALL
04E    C=0  ALL        C
35C    PT=12             =
0D0    LD@PT-3       3
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
1F8    READ7(O)
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
278    READ9(Q)
01D    ?NCXQ          C=
060     1807               A+C
268    WRIT9(Q)
1F8    READ7(O)
10E    A=C  ALL
138    READ4(L)
261    ?NCXQ           C=
060    1898                A/C
1E8    WRIT7(O)
0F8    READ3(X)
10E    A=C  ALL
0B8   READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
078    READ1(Z)
10C    ?FSET8
013    JNC+02
0B8   READ2(Y)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
238    READ8(P)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
238    READ8(P)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
04E    C=0  ALL        C
35C    PT=12             =
150    LD@PT-5        5
261    ?NCXQ           C=
060    1898                A/C
128    WRIT4(L)
2F9    ?NCXQ           C=
060    18BE                sqrt(C)
10E    A=C  ALL
138    READ4(L)
135    ?NCXQ           C=
060    184D               A*C
10E    A=C  ALL
1F8    READ7(O)
0AE   A<>C  ALL
261    ?NCXQ           C=
060    1898                A/C
10E    A=C  ALL
278    READ9(Q)
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
046    C=0 S&X        C
270    RAMSLCT      =
038    READDATA   T
2BE    C=-C-1 MS    C=-C
0AE    A<>C
028    WRIT0(T)
070    N=C  ALL
01D    ?NCXQ          C=
060     1807               A+C
10E    A=C  ALL
0B0    C=N  ALL
261    ?NCXQ           C=
060    1898                A/C
05E    C=0 MS          C= | C |
10E    A=C  ALL
04E    C=0                 |
2BE    C=-C-1 MS    |
35C    PT=12            |
090    LD@PT-2       |      C = -2 E-9
21C    PT=2              |
250    LD@PT-9      |
250    LD@PT-9      |
050    LD@PT-1      |
01D    ?NCXQ          C=
060     1807               A+C
2FE    ?C#0 MS
351    ?NCGO          if not carry goto
3EA    FAD4             loop
24C    ?FSET9
023    JNC+04
244    CLRF9
351    ?NCGO          if not carry goto
3EA    FAD4             loop
3B5    ?NCXQ         if not carry
050    14ED              R^
345    ?NCXQ         if not carry
040    10D1             CLA ( clears the alpha registers M N O P )
0F8    READ3(X)
331    ?NCGO        checks for
002    00CC            underflow/overflow
 
 

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

 
Example:

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

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

 
Example:

     7  ENTER^
     4  ENTER^
     2  ENTER^
     1  XEQ "RJZ" produces  R_J ( 1 , 2+4i , 2-4i , 7 ) = 0.126039065   ( in 11.8 seconds )
 

-If the arguments have huge ratios, the execution time may be huge too!
-For instance:  R_J(0,10 ^-50,10 ^-10,10⁵⁰) =  1.423139884 10 ^-43   is obtained in almost 3 minutes!
-Note that the other programs would lead to OUT OF RANGE though the final result is neither underflowed nor overflowed.

-If the result is really underflowed, it is replaced by zero.
-If the result is really overflowed, you'll get OUT OF RANGE, but the content of X-register will remain the same
-For example, R_J(0,10 ^-80,10 ^-80,10 ^-80) produces an OUT OF RANGE  with X-register = 2.356194491 E-80  apparently but if you press 1/X you'll get 0.

-However, if you've set flag F24,     R_J(0,10 ^-80,10 ^-80,10 ^-80)   gives   9.999999999 E99
-In fact, the correct result is 2.35619449 E+120  so the mantissa was good but the exponent was not ( and couldn't be ) exact.

-There will be a DATA ERROR if  x or y or z < 0 or p <= 0 ( with RJ )  or  x < 0  or p <= 0 ( with RJZ )
-Space can be saved if you don't check for zero-arguments...
-In this case, there will be an infinite loop if you execute RJ with x = y = 0 and z # 0  or if you execute RJZ with y = z = 0
  press  ENTER^  ON to stop the loop with a "newest" HP-41
-On older HP-41s, insert

3CC   ?KEY
360    ?C RTN    at the beginning of the loop: the infinite loop will be stopped by pressing any key for 2 seconds.

-RJ & RJZ  do not use any data register but the alpha registers M N O P are cleared and Q is also used.

-The routine stops when 2 consecutive approximations are nearly equal ( CPU-flag9 is also used by RJZ to avoid a premature stop ).

-You can certainly write a faster function if you employ the method of "RJ3" & "RJZ3":
  the number of iterations will be reduced, but the series expansion needs many extra words.

-Another improvement is to modify RJ so that it returns the Cauchy principal value of the integral if p < 0
-All depends on the place you want to allocate on your HP-41 for Carlson Elliptic Integrals ...
 

3°) 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: temp
Flags: /
Subroutine:  The M-code function RF or the focal program "RF" or "RF2".  In this case, use another register to avoid a register usage conflict.
 
 

01  LBL "LEI1" 02  1 03  P-R 04  X^2 05  X<>Y

06  STO 00        07  X^2 08  SIGN  09  LASTX 10  R^

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  RF 23  RCL 00        24  * 25  END

 
   ( 37 bytes / SIZE 001 )
 
 

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 6.5 seconds )

-If you prefer  F( phi , k ) = F ( 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°  If  90° < phi <= 180°  you'll get  F ( 180°- phi | m )
-This program works in all angular modes.
 

       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 R04: temp
Flags: /
Subroutines:  The M-code routines RF & RJ  or "RF" & "RJ" ...  In this case, use other registers to avoid a register usage conflict.
 
 

01  LBL "LEI2" 02  1 03  P-R 04  X^2 05  STO 00        06  X<>Y 07  STO 01 08  X^2 09  SIGN

10  LASTX 11  R^ 12  STO 02        13  * 14  - 15  X=Y? 16  X#0? 17  GTO 00 18  SIGN

19  RTN 20  LBL 00 21  STO 03        22  1 23  RF 24  STO 04 25  1 26  STO Y 27  RCL 00

28  RCL 03 29  RJ 30  RCL 04        31  X<>Y 32  RCL 01 33  X^2 34  * 35  RCL 02 36  *

37  3 38  / 39  - 40  RCL 01        41  * 42  END

 
   ( 55 bytes / SIZE 005 )
 
 

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 18 seconds )

-If you prefer  E( phi , k ) = E ( phi | k² )   add  X^2 after line 11 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 modes.
 

       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 R04: temp
Flags: /
Subroutines:  The M-code routines RF & RJ  or "RF" & "RJ" ...  In this case, use other registers to avoid a register usage conflict.
 
 

01  LBL "LEI3" 02  1 03  P-R 04  X^2 05  STO 00        06  X<>Y 07  STO 01 08  X^2 09  STO 02 10  R^

11  STO 03        12  X<> T 13  * 14  1 15  X<>Y 16  - 17  X=Y? 18  X#0? 19  GTO 00 20  DEG

21  90 22  TAN 23  RTN 24  LBL 00 25  STO 04        26  1 27  RF 28  X<> 02 29  RCL 03 30  *

31  1 32  ST+ Y 33  RCL 04 34  RCL 00        35  RJ 36  RCL 02  37  X<>Y 38  RCL 01 39  X^2 40  *

41  RCL 03        42  * 43  3 44  / 45  - 46  RCL 01 47  * 48  END

 
   ( 64 bytes / SIZE 005 )
 
 

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.336853615    ( in 19 seconds )

-If you prefer  ¶ ( phi , n , k ) = ¶ ( n ; phi | k² )   add  X^2 after line 12 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 modes.
->   1 + n.sin²(phi) is supposed to be positive ( otherwise, use "RJ" instead of RJ )
 

       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 R04: temp

Flags:   F01    Set flag F01 to compute Legendre elliptic integrals of the first kind
             F02     Set flag F02 to compute Legendre elliptic integrals of the second kind         Set only one of these 2 flags
             Clear   F01 & F02 to compute Legendre elliptic integrals of the third kind

Subroutines: The M-code routines RF & RJ   or "RF" & "RJ" ...  In this case, use other registers to avoid a register usage conflict.

-Add  X^2  after line 11  if you prefer  k  instead of  m
 
 

01  LBL "LEI" 02  1 03  P-R 04  X^2 05  STO 01 06  X<>Y 07  STO 02        08  X^2 09  R^ 10  STO 00  11  X<> T 12  FS? 02

13  STO 00  14  * 15  1 16  X<>Y 17  - 18  X=Y? 19  X#0? 20  GTO 00        21  SIGN 22  FS? 02 23  RTN 24  DEG

25  90 26  TAN 27  RTN 28  LBL 00 29  STO 03  30  1 31  RF 32  FS? 01 33  GTO 00 34  STO 04        35  RCL 02 36  X^2

37  RCL 00  38  * 39  1 40  ST+ Y 41  FS? 02 42  STO Y 43  RCL 03 44  RCL 01        45  RJ 46  RCL 02 47  X^2 48  *

49  RCL 00  50  * 51  3 52  / 53  RCL 04        54  X<>Y 55  - 56  LBL 00 57  RCL 02 58  * 59  END

 
   ( 79 bytes / SIZE 005 )
 
 

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
      Leg-Ell-Int  =  F ( phi | m )        if flag  F02 is set
      Leg-Ell-Int  =  ¶ ( n ; phi | m )    if flags F01 & F02 are cleared.
 
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  B.C. Carlson, "Numerical Computation of Real or Complex Elliptic Integrals"
[3]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4