Lauricella Functions for the HP-41

Overview
 

 1°)  Lauricella Function A

   a)  3 variables
   b)  n variables

 2°)  Lauricella Function B

   a)  3 variables
   b)  n variables

 3°)  Lauricella Function C

   a)  3 variables
   b)  n variables

 4°)  Lauricella Function D

   a)  3 variables: Focal Program
   b)  3 variables: M-Code Routine
   c)  Legendre Elliptic Integral of the 3rd kind
   d)  n variables
 

-These functions are defined by hypergeometric series of n variables.
-If n = 2, they correspond to the Appell hypergeometric functions F₂ F₃ F₄ F₁.
-If n = 1, all these 4 functions are equal to the ( Gauss ) hypergeometric series ₂F₁
 

1°)  Lauricella Function A
 

     a)  3 Variables
 

-This function is defined by

  F_A( a ; b₁ , b₂ , b₃ ; c₁ , c₂ , c₃ ; x , y , z ) = SUM_{m,n,p=0 to infinity}   (a)_m+n+p (b₁)_m (b₂)_n (b₃)_p / ((c₁)_m (c₂)_n (c₃)_p)  x^myⁿz^p / (m! n! p!)

  where   (a)_k = a(a+1)(a+2) ....... (a+k-1)   if  k > 0   and   (a)₀ =  1    ( Pochhammer symbol )
 

Data Registers:               R00 , R08 to R10 , R14: temp               ( Registers R01 thru R07 are to be initialized before executing "LCFA3" )

                                      •  R01 = a       •  R02 = b₁       •  R05 = c₁      R11 = x
                                                            •  R03 = b₂       •  R06 = c₂      R12 = y
                                                            •  R04 = b₃       •  R07 = c₃      R13 = z
Flags: /
Subroutines: /
 
 

01  LBL "LCFA3" 02  STO 11          03  RDN 04  STO 12 05  X<>Y 06  STO 13 07  CLX 08  STO 08 09  STO 09 10  STO 10 11  SIGN 12  STO 00 13  STO 14 14  SIGN 15  LBL 01 16  LASTX 17  RCL 01 18  RCL 08 19  + 20  RCL 09

21  + 22  RCL 10          23  + 24  * 25  RCL 02 26  RCL 08 27  + 28  * 29  RCL 11 30  * 31  RCL 05 32  RCL 08 33  ST+ Y 34  ISG X 35  CLX 36  STO 08 37  * 38  / 39  + 40  X#Y?

41  GTO 01 42  CLX 43  STO 08  44  X<> 00          45  RCL 01 46  RCL 09 47  + 48  RCL 10 49  + 50  * 51  RCL 03 52  RCL 09 53  + 54  * 55  RCL 12 56  * 57  RCL 06 58  RCL 09 59  ST+ Y 60  ISG X

61  CLX 62  STO 09 63  * 64  / 65  STO 00          66  + 67  X#Y? 68  GTO 01 69  CLX 70  STO 08 71  STO 09 72  X<> 14 73  RCL 01 74  RCL 10 75  + 76  * 77  RCL 04 78  RCL 10 79  + 80  *

81  RCL 13 82  * 83  RCL 07          84  RCL 10  85  ST+ Y 86  ISG X 87  CLX 88  STO 10 89  * 90  / 91  STO 00 92  STO 14 93  + 94  X#Y? 95  GTO 01 96  END

 
  ( 117 bytes / SIZE 015 )
 
 

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

 
Example:    With   a = 0.7 ,  b₁ = 0.1  b₂ = 0.2  b₃ = 0.3 ,  c₁ = 2.4  c₂ = 2.5  c₃ = 2.6       ( Store these 7 numbers into R01 thru R07 )

   0.13   ENTER^
   0.12   ENTER^
   0.11   XEQ "LCFA3"  >>>>>   F_A( 0.11 , 0.12 , 0.13 ) = 1.021571587                ---Execution time = 104s---

   0.15   ENTER^
   0.25   ENTER^
   0.35      R/S        >>>>>   F_A( 0.35 , 0.25 , 0.15 ) = 1.040641882                ---Execution time = 4mn28s---
 

Note:

-The series is convergent if  | x | + | y | + | z |  <  1
 

     b)  n Variables
 

-This function may be generalized to n variables as follows:
 

  F_A( a ; b₁ , ..... , b_n ; c₁ , ..... , c_n ; x₁ , ..... , x_n ) = SUM_{i1,...,in=0 to infinity}   (a)_i1+...+in (b₁)_i1 ...... (b_n)_in / ((c₁)_i1 ... (c_n)_in)  x₁ⁱ¹ ... x_nⁱⁿ / ( i₁! .... i_n !)
 
 

Data Registers:           •  R00 = n               ( Registers R00 thru R3n+1 are to be initialized before executing "LCFA" )

                                      •  R01 = x₁     •  Rn+1 = a       •  Rn+2 = b₁       •  R2n+2 = c₁               R3n+2 to R5n+1: temp
                                         .............                                ............................................

                                      •  Rnn = x_n                              •  R2n+1 = b_n     •  R3n+1 = c_n

Flag:   F10  is cleared at the end
Subroutines: /
 

 2 M-code routines are used  X+1   to add 1 to register X
                                    and   X/E3  to divide X by 1000

  X+1  may be replaced by  ISG X  CLX
  X/E3 may be replaced by   PI   INT   10^X   /    ( we cannot use  E3  /   because  synthetic registers  P & Q are used )

-Of course, synthetic registers  M  N  O  P  Q  may be replaced by standard registers - for instance R95  R96  R97  R98  R99 - provided n is not too large
-The CX function  CLRGX  is also used.

-Line 113 is a three-byte GTO 12
 
 

01  LBL "LCFA"   02  RCL 00           03  STO N   04  5   05  *   06  X+1    07  STO O   08  .1   09  %   10  +   11  RCL 00   12  -   13  CLRGX   14  INT   15  STO P   16  RCL 00   17  -   18  STO Q   19  RCL P   20  X/E3   21  ST+ O   22  +   23  X+1   24  SIGN

25  STO M   26  LBL 00   27  STO IND L   28  ISG L   29  GTO 00   30  LBL 12   31  SF 10   32  LBL 01   33  RCL 00           34  X+1   35  RCL O   36  RCL IND Y   37  LBL 02   38  RCL IND Y   39  +   40  DSE Y   41  GTO 02   42  RCL Q   43  RCL 00   44  -   45  X<>Y   46  RCL IND Y   47  RCL IND O   48  +

49  *   50  RCL IND Q   51  RCL IND O   52  ST+ Y   53  X+1   54  STO IND O   55  *   56  /   57  RCL IND N   58  *   59  ST* IND P   60  RCL IND P   61  RCL M   62  +   63  STO M   64  LASTX   65  FC? 10   66  GTO 04           67  X#Y?   68  GTO 01   69  LBL 03   70  RCL 00   71  ENTER^   72  ST+ X

73  ST+ X   74  +   75  X+1   76  X/E3   77  RCL O   78  INT   79  +   80  CLRGX    81  DSE N   82  X=0?   83  GTO 06   84  DSE O   85  DSE P   86  DSE Q   87  CF 10   88  GTO 01   89  LBL 04   90  RCL 00           91  ST+ X   92  ST+ X   93  X+1   94  X/E3   95  RCL P   96  +

97  RCL IND X   98  LBL 05   99  STO IND Y 100  ISG Y 101  GTO 05 102  R^ 103  R^ 104  X=Y? 105  GTO 03 106  RCL 00         107  ENTER^ 108  X<> N 109  - 110  ST+ O 111  ST+ P 112  ST+ Q 113  GTO 12 114  LBL 06 115  X<> M 116  CLA 117  END

 
 
    ( 194 bytes / SIZE 5n+2 )
 
 

STACK	INPUTS	OUTPUTS
X	/	FA(x1 , ... , xn)

 
Example:   the same one:    n = 3  STO 00

     x₁ = 0.11   STO 01        a = 0.7   STO 04      b₁ = 0.1   STO 05        c₁ = 2.4    STO 08
     x₂ = 0.12   STO 02                                        b₂ = 0.2   STO 06        c₂ = 2.5    STO 09
     x₃ = 0.13   STO 03                                        b₃ = 0.3   STO 07        c₃ = 2.6    STO 10

  XEQ "LCFA"  >>>>   F_A( 0.11 , 0.12 , 0.13 ) = 1.021571587                ---Execution time = 3m53s---

Note:

-The series is convergent if   | x₁ | + ...... + | x_n |  <  1
 

2°)  Lauricella Function B
 

     a)  3 Variables
 

-This function is defined by

  F_B( a₁ , a₂ , a₃ ; b₁ , b₂ , b₃ ; c ; x , y , z ) = SUM_{m,n,p=0 to infinity}    (a₁)_m (a₂)_n (a₃)_p (b₁)_m (b₂)_n (b₃)_p / (c)_m+n+p   x^myⁿz^p / (m! n! p!)
 

Data Registers:               R00 , R08 to R10 , R14: temp               ( Registers R01 thru R07 are to be initialized before executing "LCFB3" )

                                      •  R01 = a₁       •  R04 = b₁       •  R07 = c         R11 = x
                                      •  R02 = a₂       •  R05 = b₂                                R12 = y
                                      •  R03 = a₃       •  R06 = b₃                                R13 = z
Flags: /
Subroutines: /
 
 

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

21  RCL 04 22  RCL 08          23  + 24  * 25  RCL 11 26  * 27  RCL 07 28  RCL 09 29  + 30  RCL 10  31  + 32  RCL 08 33  ST+ Y 34  ISG X 35  CLX 36  STO 08 37  * 38  / 39  + 40  X#Y?

41  GTO 01 42  CLX 43  STO 08          44  X<> 00 45  RCL 02 46  RCL 09 47  + 48  * 49  RCL 05 50  RCL 09  51  + 52  * 53  RCL 12 54  * 55  RCL 07 56  RCL 10 57  + 58  RCL 09 59  ST+ Y 60  ISG X

61  CLX 62  STO 09          63  * 64  / 65  STO 00 66  + 67  X#Y? 68  GTO 01 69  CLX 70  STO 08 71  STO 09  72  X<> 14 73  RCL 03 74  RCL 10 75  + 76  * 77  RCL 06 78  RCL 10 79  + 80  *

81  RCL 13          82  * 83  RCL 07 84  RCL 10 85  ST+ Y 86  ISG X 87  CLX 88  STO 10  89  * 90  / 91  STO 00 92  STO 14 93  + 94  X#Y? 95  GTO 01 96  END

 
   ( 117 bytes / SIZE 015 )
 
 

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

 
Example:    With   a₁ = 0.1  a₂ = 0.2  a₃ = 0.3 ,  b₁ = 0.4  b₂ = 0.5  b₃ = 0.6 ,  c = 2.7      ( Store these 7 numbers into R01 thru R07 )

   0.43   ENTER^
   0.42   ENTER^
   0.41   XEQ "LCFB3"  >>>>>   F_B( 0.41 , 0.42 , 0.43 ) = 1.057748024                ---Execution time = 2m38s---

   0.73   ENTER^
   0.72   ENTER^
   0.71      R/S        >>>>>   F_B( 0.71 , 0.72 , 0.73 ) = 1.112251899                ---Execution time = 6mn35s---
 

Note:

-The series is convergent if   Max { | x | , | y | , | z |  }  <  1
 

     b)  n Variables
 

F_B( a₁ , ..... , a_n ; b₁ , ..... , b_n ; c ; x₁ , ..... , x_n ) = SUM_{i1,...,in=0 to infinity}  (a₁)_i1 ...... (a_n)_in (b₁)_i1 ... (b_n)_in / (c)_i1+...+in     x₁ⁱ¹ ... x_nⁱⁿ / ( i₁! .... i_n!)
 
 

Data Registers:           •  R00 = n               ( Registers R00 thru R3n+1 are to be initialized before executing "LCFB" )

                                      •  R01 = x₁     •  Rn+1 = a₁       •  R2n+1 = b₁         •  R3n+1 = c             R3n+2 to R5n+1: temp
                                         ..................................................................

                                      •  Rnn = x_n      •  R2n = a_n        •  R3n = b_n

Flag:   F10  is cleared at the end
Subroutines: /
 

 2 M-code routines are used  X+1   to add 1 to register X
                                    and   X/E3  to divide X by 1000

  X+1  may be replaced by  ISG X  CLX
  X/E3 may be replaced by   PI   INT   10^X   /

-Line 35 may be replaced by  ENTER^   ST+ X   +
-Line 115 is a three-byte GTO 12
 
 

01  LBL "LCFB"   02  RCL 00           03  STO N   04  5   05  *   06  X+1    07  STO O   08  .1   09  %   10  +   11  RCL 00   12  -   13  CLRGX   14  INT   15  STO P   16  RCL 00   17  -   18  STO Q   19  RCL P   20  X/E3   21  ST+ O   22  +   23  X+1   24  SIGN

25  STO M   26  ST- Q   27  LBL 00   28  STO IND L   29  ISG L   30  GTO 00           31  LBL 12   32  SF 10   33  LBL 01   34  RCL 00   35  3X    36  X+1   37  RCL O   38  RCL IND Y   39  LBL 02   40  RCL IND Y   41  +   42  DSE Y   43  GTO 02   44  RCL Q   45  RCL 00   46  -   47  RCL IND X   48  RCL IND O

49  +   50  R^   51  /   52  RCL IND Q   53  RCL IND O   54  ST+ Y   55  X+1   56  STO IND O   57  /   58  *   59  RCL IND N   60  *   61  ST* IND P   62  RCL IND P   63  RCL M           64  +   65  STO M   66  LASTX   67  FC? 10   68  GTO 04   69  X#Y?   70  GTO 01   71  LBL 03   72  RCL 00

73  ENTER^   74  ST+ X   75  ST+ X   76  +   77  X+1   78  X/E3   79  RCL O           80  INT   81  +   82  CLRGX    83  DSE N   84  X=0?   85  GTO 06   86  DSE O   87  DSE P   88  DSE Q   89  CF 10   90  GTO 01   91  LBL 04   92  RCL 00   93  ST+ X   94  ST+ X   95  X+1   96  X/E3

97  RCL P   98  +   99  RCL IND X 100  LBL 05 101  STO IND Y 102  ISG Y 103  GTO 05         104  R^ 105  R^ 106  X=Y? 107  GTO 03 108  RCL 00 109  ENTER^ 110  X<> N 111  - 112  ST+ O 113  ST+ P 114  ST+ Q 115  GTO 12 116  LBL 06 117  X<> M 118  CLA 119  END

 
   ( 198 bytes / SIZE 5n+2 )
 
 

STACK	INPUTS	OUTPUTS
X	/	FB(x1 , ... , xn)

 
Example:                 n = 3  STO 00

     x₁ = 0.41   STO 01           a₁ = 0.1   STO 04        b₁ = 0.4    STO 07       c = 2.7   STO 10
     x₂ = 0.42   STO 02           a₂ = 0.2   STO 05        b₂ = 0.5    STO 08
     x₃ = 0.43   STO 03           a₃ = 0.3   STO 06        b₃ = 0.6    STO 09

  XEQ "LCFB"  >>>>   F_B( 0.41 , 0.42 , 0.43 ) = 1.057748024                ---Execution time = 6m14s---

Note:

-The series is convergent if   Max { | x₁ | ,......, | x_n | }  <  1
 

3°)  Lauricella Function C
 

     a)  3 Variables
 

-This function is defined by

  F_C( a ; b ; c₁ , c₂ , c₃ ;  x , y , z ) = SUM_{m,n,p=0 to infinity}    (a)_m+n+p (b)_m+n+p  / ( (c₁)_m(c₂)_n(c₃)_p )  x^myⁿz^p / (m! n! p!)
 

Data Registers:               R00 , R06 to R10: temp               ( Registers R01 thru R05 are to be initialized before executing "LCFC3" )

                                      •  R01 = a        •  R03 = b₁                    R11 = x
                                      •  R02 = b        •  R04 = b₂                    R12 = y
                                                              •  R05 = b₃                    R13 = z
Flags: /
Subroutines: /
 
 

01  LBL "LCFC3" 02  STO 11          03  RDN 04  STO 12 05  X<>Y 06  STO 13 07  CLX 08  STO 08 09  STO 09 10  STO 10 11  SIGN 12  STO 00 13  STO 06 14  SIGN 15  LBL 01 16  LASTX 17  RCL 08 18  RCL 09 19  + 20  RCL 10

21  + 22  STO 07          23  RCL 01 24  + 25  * 26  RCL 02 27  RCL 07  28  + 29  * 30  RCL 11 31  * 32  RCL 03 33  RCL 08 34  ST+ Y 35  ISG X 36  CLX 37  STO 08 38  * 39  / 40  +

41  X#Y? 42  GTO 01 43  CLX 44  STO 08          45  X<> 00 46  RCL 09 47  RCL 10  48  + 49  STO 07 50  RCL 01 51  + 52  * 53  RCL 02 54  RCL 07 55  + 56  * 57  RCL 12 58  * 59  RCL 04 60  RCL 09

61  ST+ Y 62  ISG X 63  CLX 64  STO 09          65  * 66  / 67  STO 00  68  + 69  X#Y? 70  GTO 01 71  CLX 72  STO 08 73  STO 09 74  X<> 06 75  RCL 01 76  RCL 10 77  + 78  * 79  RCL 02 80  RCL 10

81  + 82  * 83  RCL 13          84  * 85  RCL 05 86  RCL 10  87  ST+ Y 88  ISG X 89  CLX 90  STO 10 91  * 92  / 93  STO 00 94  STO 06 95  + 96  X#Y? 97  GTO 01 98  END

 
   ( 119 bytes / SIZE 014 )
 
 

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

 
Example:    With   a = 0.1 , b = 0.2 ,  c₁ = 2.7   c₂ = 2.8   c₃ = 2.9         ( Store these 5 numbers into R01 thru R05 )

   0.12   ENTER^
   0.11   ENTER^
   0.10   XEQ "LCFC3"  >>>>>   F_C( 0.10 , 0.11 , 0.12 ) = 1.002554253                ---Execution time = 2m34s---

Notes:

-The series is convergent if     | x |^1/2 + | y |^1/2 + | z |^1/2   <  1
-This program may return a meaningless result if the above condition is not satisfied.
 

     b)  n Variables
 

F_C( a ; b ; c₁ , ..... , c_n ; x₁ , ..... , x_n ) = SUM_{i1,...,in=0 to infinity}  (a)_i1+...+in (b)_i1+...+in) / ( (c₁)_i1 ...... (c_n)_in)    x₁ⁱ¹ ... x_nⁱⁿ / ( i₁! .... i_n!)
 
 

Data Registers:           •  R00 = n               ( Registers R00 thru R3n+1 are to be initialized before executing "LCFC" )

                                      •  R01 = x₁     •  Rn+1 = a       •  Rn+2 = b         •  Rn+3 = c₁             R2n+3 to R4n+2: temp
                                         .............                                                            .................

                                      •  Rnn = x_n                                                          •  R2n+2 = c_n

Flag:   F10  is cleared at the end
Subroutines: /
 

 2 M-code routines are used  X+1   to add 1 to register X
                                    and   X/E3  to divide X by 1000

  X+1  may be replaced by  ISG X  CLX
  X/E3 may be replaced by  PI   INT   10^X   /

-Line 118 is a three-byte GTO 12
 
 

01  LBL "LCFC"   02  RCL 00           03  STO N   04  4   05  *   06  2   07  +   08  STO O   09  .1   10  %   11  +   12  RCL 00   13  -   14  CLRGX   15  INT   16  STO P   17  RCL 00   18  -   19  STO Q   20  RCL P   21  X/E3   22  ST+ O   23  +   24  X+1   25  SIGN

26  STO M   27  LBL 00   28  STO IND L   29  ISG L   30  GTO 00   31  LBL 12   32  SF 10   33  LBL 01           34  RCL O   35  ENTER^   36  CLX   37  LBL 02   38  RCL IND Y   39  +   40  DSE Y   41  GTO 02   42  RCL 00   43  SIGN   44  ST+ L   45  RCL IND L   46  X<>Y   47  ST+ L   48  RCL IND L   49  R^   50  ST+ T

51  +   52  R^   53  *   54  RCL IND Q   55  RCL IND O   56  ST+ Y   57  X+1   58  STO IND O   59  *   60  /   61  RCL IND N   62  *   63  ST* IND P   64  RCL IND P   65  RCL M   66  +   67  STO M   68  LASTX   69  FC? 10   70  GTO 04   71  X#Y?   72  GTO 01           73  LBL 03   74  RCL 00   75  ENTER^

76  ST+ Y   77  X+1   78  ST+ X   79  +   80  X/E3   81  RCL O   82  INT   83  +   84  CLRGX    85  DSE N   86  X=0?   87  GTO 06   88  DSE O   89  DSE P   90  DSE Q   91  CF 10   92  GTO 01           93  LBL 04   94  RCL 00   95  ENTER^   96  X+1   97  ST+ X   98  +   99  X/E3 100  RCL P

101  + 102  RCL IND X 103  LBL 05 104  STO IND Y 105  ISG Y 106  GTO 05         107  R^ 108  R^ 109  X=Y? 110  GTO 03 111  RCL 00 112  ENTER^ 113  X<> N 114  - 115  ST+ O 116  ST+ P 117  ST+ Q 118  GTO 12 119  LBL 06 120  X<> M 121  CLA 122  END

 
   ( 197 bytes / SIZE 4n+3 )
 
 

STACK	INPUTS	OUTPUTS
X	/	FC(x1 , ... , xn)

 
Example:               n = 3  STO 00

     x₁ = 0.10   STO 01           a = 0.1   STO 04        b = 0.2    STO 05       c₁ = 2.7   STO 06
     x₂ = 0.11   STO 02                                                                                c₂ = 2.8    STO 07
     x₃ = 0.12   STO 03                                                                                c₃ = 2.9    STO 08

  XEQ "LCFC"  >>>>   F_C( 0.10 , 0.11 , 0.12 ) = 1.002554252                ---Execution time = 5m40s---

Notes:

-The series is convergent if     | x₁ |^1/2 + ........  + | x_n |^1/2   <  1
-This program may return a meaningless result if the above condition is not satisfied.
 

4°)  Lauricella Function D
 

     a)  3 Variables: Focal Program
 

-This function is defined by

  F_D( a ; b₁ , b₂ , b₃ ; c ;  x , y , z ) = SUM_{m,n,p=0 to infinity}    (a)_m+n+p (b₁)_m (b₂)_n (b₃)_p  /  (c)_m+n+p   x^myⁿz^p / (m! n! p!)
 

Data Registers:               R00 , R06 to R10: temp                    ( Registers R01 thru R05 are to be initialized before executing "LCFD3" )

                                      •  R01 = a        •  R02 = b₁       •  R05 = c            R11 = x
                                                              •  R03 = b₂                                   R12 = y
                                                              •  R04 = b₃                                   R13 = z
Flags: /
Subroutines: /
 
 

01  LBL "LCFD3" 02  STO 11          03  RDN 04  STO 12 05  X<>Y 06  STO 13 07  CLX 08  STO 08 09  STO 09 10  STO 10 11  SIGN 12  STO 00 13  STO 06 14  SIGN 15  LBL 01 16  LASTX 17  RCL 08 18  RCL 09 19  + 20  RCL 10

21  + 22  STO 07          23  RCL 01 24  + 25  * 26  RCL 05 27  RCL 07 28  + 29  / 30  RCL 11 31  * 32  RCL 02 33  RCL 08 34  ST+ Y 35  ISG X 36  CLX 37  STO 08  38  / 39  * 40  +

41  X#Y? 42  GTO 01 43  CLX 44  STO 08          45  X<> 00 46  RCL 09 47  RCL 10 48  + 49  STO 07 50  RCL 01 51  + 52  * 53  RCL 05 54  RCL 07 55  + 56  / 57  RCL 12  58  * 59  RCL 03 60  RCL 09

61  ST+ Y 62  ISG X 63  CLX 64  STO 09          65  / 66  * 67  STO 00 68  + 69  X#Y? 70  GTO 01 71  CLX 72  STO 08 73  STO 09 74  X<> 06 75  RCL 01 76  RCL 10  77  + 78  * 79  RCL 05 80  RCL 10

81  + 82  / 83  RCL 13          84  * 85  RCL 04 86  RCL 10 87  ST+ Y 88  ISG X 89  CLX 90  STO 10 91  / 92  * 93  STO 00 94  STO 06  95  + 96  X#Y? 97  GTO 01 98  END

 
  ( 119 bytes / SIZE 014 )
 
 

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

 
Example:      With   a = 0.1 ,  b₁ = 0.2   b₂ = 0.3   b₃ = 0.4 ,  c = 2.7             ( Store these 5 numbers into R01 thru R05 )

   0.41   ENTER^
   0.31   ENTER^
   0.21   XEQ "LCFD3"  >>>>>   F_D( 0.21 , 0.31 , 0.41 ) = 1.012345372                ---Execution time = 2m53s---

    0.7    ENTER^
    0.6    ENTER^
    0.5       R/S        >>>>>   F_D( 0.5 , 0.6 , 0.7 ) = 1.025992608                ---Execution time = 11m10s---
 

Notes:

-The series is convergent if    Max { | x | , | y | , | z |  }  <  1

-There are also 10 other Saran-Lauricella functions of 3 variables which may be programmed in a similar way.
-All these triple series are special cases of a generalized hypergeometric series of 3 variables.
-See for instance "Srivastava's Functions of Three Variables for the HP-41"
 

     b)  3 Variables: M-Code Routine
 

-We'll get faster - or at least less slow - results with M-Code.
-Registers R01 thru R05 are replaced by R00 thru R04.
-Their existence is tested but there is no check for alpha data.

-The synthetic registers M & N are denoted m & n in the simplified mnemonics below.
-The alpha "register" is cleared.
 

Data Registers:                                  ( Registers R00 thru R04 are to be initialized before executing "LCFD3" )

                                      •  R00 = a        •  R01 = b₁       •  R04 = c
                                                              •  R02 = b₂
                                                              •  R03 = b₃
Flags: /
Subroutines: /
 

0B3  "3"
004  "D"
006  "F"
003  "C"
00C  "L"
3C8  CLRKEY        first executable word
378   C=c
03C  RCR 3
106  A=C S&X
130  LDI S&X
1FC 508d               correct value for an HP-41CV/CX or C with a quad memory module
306  ?A<C S&X
381  ?NCGO
00A  02E0              ( NONEXISTENT )
2A0  SETDEC
345   ?NCXQ
040   CLA
0F8   C=X
128   L=C
078   C=Z
028   T=C
04E   C
35C   =
050   1
0E8   X=C
068   Z=C
228   P=C
268   Q=C
178   C=M                LOOP  at the address  E48D  in my ROM - Change the 3 ?CGO E48D  written in red below according to your own ROM
10E  A=C ALL
1B8  C=n
01D  C=
060   A+C
1F8  C=O
025  C=
060  AB+C
070  N=C ALL
378  C=c
03C  RCR 3
270  RAMSLCT
038  READATA
025  C=
060  AB+C
046  C=0 S&X
270  RAMSLCT
078  C=Z
13D  C=
060   AB*C
068   Z=C
378   C=c
03C  RCR 3
260  SETHEX
226  C=C+1 S&X
226  C=C+1 S&X
226  C=C+1 S&X
226  C=C+1 S&X
270  RAMSLCT
038  READATA
10E  A=C ALL
0F0  C<>N ALL
2A0  SETDEC
01D  C=
060   A+C
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
078   C=Z
0AE  A<>C ALL
261   C=
060   A/C
138   C=L
13D  C=
060   AB*C
068   Z=C
178   C=m
10E   A=C ALL
378   C=c
03C  RCR 3
260   SETHEX
226   C=C+1 S&X
270  RAMSLCT
038  READATA
2A0  SETDEC
01D  C=
060   A+C
046   C=0 S&X
270   RAMSLCT
078   C=Z
13D  C=
060   AB*C
068   Z=C
178   C=m
00E   A
35C  =
162   1
01D  C=
060   A+C
168   m=C
10E   A=C ALL
078   C=Z
0AE  A<>C ALL
261   C=
060   A/C
068   Z=C
3CC  ?KEY                               These 2 words allow you to stop the routine
360   ?C RTN                             if it's too long: press any key to stop it.
0F8   C=X
025   C=
060   AB+C
10E   A=C ALL
0F8   C=X
0AE  A<>C ALL
0E8   X=C
36E   ?A#C ALL
235   ?CGO                  ?CGO
393    E48D                  LOOP
1B8   C=n
10E   A=C ALL
1F8   C=O
01D  C=
060   A+C
168   m=C
378   C=c
03C  RCR 3
270  RAMSLCT
038  READATA
025   C=
060   AB+C
046   C=0 S&X
270   RAMSLCT
238   C=P
13D  C=
060   AB*C
228   P=C
178   C=m
10E   A=C ALL
0B0   C=N ALL
01D  C=
060   A+C
10E   A=C ALL
238   C=P
0AE  A<>C ALL
261   C=
060   A/C
0B8  C=Y
13D  C=
060   AB*C
228   P=C
378   C=c
03C  RCR 3
260  SETHEX
226  C=C+1 S&X
226  C=C+1 S&X
270  RAMSLCT
038  READATA
10E  A=C ALL
04E  C=0 ALL
270   RAMSLCT
168   m=C
2A0  SETDEC
1B8  C=n
01D  C=
060   A+C
238   C=P
13D  C=
060   AB*C
228   P=C
1B8  C=n
00E   A
35C  =
162   1
01D  C=
060   A+C
1A8   n=C
10E   A=C ALL
238   C=P
0AE  A<>C ALL
261   C=
060   A/C
228   P=C
068   Z=C
0F8   C=X
025   C=
060   AB+C
10E   A=C ALL
0F8   C=X
0AE  A<>C ALL
0E8   X=C
36E   ?A#C ALL
235   ?CGO                  ?CGO
393    E48D                  LOOP
1F8   C=O
10E   A=C ALL
378   C=c
03C  RCR 3
270  RAMSLCT
038  READATA
01D  C=
060   A+C
046   C=0 S&X
270   RAMSLCT
278   C=Q
13D  C=
060   AB*C
268   Q=C
0B0  C=N ALL
10E   A=C ALL
1F8   C=O
01D  C=
060   A+C
10E   A=C ALL
278   C=Q
0AE  A<>C ALL
261   C=
060   A/C
04E   C=0 ALL
1A8  n=C
270   RAMSLCT
038   C=T
13D  C=
060   AB*C
268   Q=C
378   C=c
03C  RCR 3
260  SETHEX
226  C=C+1 S&X
226  C=C+1 S&X
226  C=C+1 S&X
270  RAMSLCT
038  READATA
10E  A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
1F8   C=O
01D  C=
060   A+C
278   C=Q
13D  C=
060   AB*C
268   Q=C
1F8   C=O
00E   A
35C  =
162   1
01D  C=
060   A+C
1E8   O=C
10E   A=C ALL
278   C=Q
0AE  A<>C ALL
261   C=
060   A/C
268   Q=C
228   C=P
068   Z=C
0F8   C=X
025   C=
060   AB+C
10E   A=C ALL
0F8   C=X
0AE  A<>C ALL
0E8   X=C
36E   ?A#C ALL
235   ?CGO                  ?CGO
393    E48D                  LOOP
046    C
270    =
038    T
068    Z=C
345    ?NCGO
042    CLA

( 271 words / SIZE 005 )
 
 

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

 
Example:      With   a = 0.1 ,  b₁ = 0.2   b₂ = 0.3   b₃ = 0.4 ,  c = 2.7                 ( Store these 5 numbers into R00 thru R04 )

   0.41   ENTER^
   0.31   ENTER^
   0.21   XEQ "LCFD3"  >>>>>   F_D( 0.21 , 0.31 , 0.41 ) = 1.012345372                ---Execution time = 1m06s---

Note:

-The M-Code routine is almost 3 times as fast as the focal program
  but it remains relatively slow...
 

     c)  Legendre Elliptic Integral of the 3rd Kind
 

-Legendre elliptic integrals of the 3rd kind may be expressed in terms of Lauricella functions F_D
-More explicitly:

      ¶ ( n ; phi | m ) = §₀^phi ( 1 + n sin² t ) ^-1 ( 1 - m sin² t ) ^-1/2 dt
                             = sin (phi) . F_D ( 1/2 ; 1 , 1/2 , 1/2 ; 3/2 ;  - n sin² (phi) , sin² (phi) , m sin² (phi) )
 
 

01  LBL "LEI3"  02  SIN  03  STO 05  04  X^2  05  ST* Y  06  ST* Z  07  .5  08  STO 00  09  STO 02  10  STO 03  11  STO 04  12  SIGN  13  STO 01  14  ST+ 04  15  X<> T  16  CHS  17  LCFD3  18  RCL 05  19  *  20  END

 
( 35 bytes / SIZE 006 )
 
 

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

 
Example:     in DEG mode:

       0.3  ENTER^
       0.7  ENTER^
       41   XEQ "LEI3"  >>>>  ¶ ( 0.3 ; 41° | 0.7 ) =  0.726768062                        ---Execution time = 3mn14s---

Notes:

-Obviously, this version remains very slow, even with the M-Code LCFD3.
-So, it cannot compete with Carlson elliptic integrals.
-It's just given here for completeness.

-Like in "Numerical Recipes", the sign convention for n is opposite that of Abramowitz & Stegun
-Otherwise, simply delete line 16 CHS
-The angle phi must be between  -90° and +90°
 

     d)  n Variables
 

 F_D( a ; b₁ , ..... , b_n ; c ; x₁ , ..... , x_n ) = SUM_{i1,...,in=0 to infinity}  (a)_i1+....+in (b₁)_i1 ... (b_n)_in / (c)_i1+...+in     x₁ⁱ¹ ... x_nⁱⁿ / ( i₁! .... i_n!)
 
 

Data Registers:           •  R00 = n               ( Registers R00 thru R2n+2 are to be initialized before executing "LCFD" )

                                      •  R01 = x₁     •  Rn+1 = a        •  Rn+2 = b₁         •  R2n+2 = c             R2n+3 to R4n+2: temp
                                         ...........                                      .................

                                      •  Rnn = x_n                                •  R2n+1 = b_n

Flag:   F10  is cleared at the end
Subroutines: /
 

 2 M-code routines are used  X+1   to add 1 to register X
                                    and   X/E3  to divide X by 1000

  X+1  may be replaced by  ISG X  CLX
  X/E3 -------------------   PI   INT   10^X   /

-Line 119 is a three-byte GTO 12
 
 

01  LBL "LCFD"   02  RCL 00           03  STO N   04  4   05  *   06  2   07  +   08  STO O   09  .1   10  %   11  +   12  RCL 00   13  -   14  CLRGX   15  INT   16  STO P   17  RCL 00   18  -   19  STO Q   20  RCL P   21  X/E3   22  ST+ O   23  +   24  X+1   25  SIGN

26  STO M   27  ST- Q   28  LBL 00   29  STO IND L   30  ISG L   31  GTO 00   32  LBL 12   33  SF 10   34  LBL 01   35  RCL O   36  ENTER^   37  CLX   38  LBL 02   39  RCL IND Y   40  +   41  DSE Y   42  GTO 02   43  RCL 00           44  SIGN   45  ST+ L   46  RCL IND L   47  LASTX   48  ST+ X   49  RDN   50  RCL IND T

51  R^   52  ST+ Z   53  +   54  /   55  RCL IND Q   56  RCL IND O   57  ST+ Y   58  X+1   59  STO IND O   60  /   61  *   62  RCL IND N   63  *   64  ST* IND P   65  RCL IND P   66  RCL M   67  +   68  STO M   69  LASTX   70  FC? 10   71  GTO 04           72  X#Y?   73  GTO 01   74  LBL 03   75  RCL 00

76  ENTER^   77  ST+ Y   78  X+1   79  ST+ X   80  +   81  X/E3   82  RCL O   83  INT   84  +   85  CLRGX   86  DSE N   87  X=0?   88  GTO 06    89  DSE O   90  DSE P   91  DSE Q   92  CF 10   93  GTO 01   94  LBL 04   95  RCL 00           96  ENTER^   97  X+1   98  ST+ X   99  + 100  X/E3

101  RCL P 102  + 103  RCL IND X 104  LBL 05 105  STO IND Y 106  ISG Y 107  GTO 05 108  R^ 109  R^ 110  X=Y? 111  GTO 03 112  RCL 00         113  ENTER^ 114  X<> N 115  - 116  ST+ O 117  ST+ P 118  ST+ Q 119  GTO 12 120  LBL 06 121  X<> M 122  CLA 123  END

 
   ( 199 bytes / SIZE 4n+3 )
 
 

STACK	INPUTS	OUTPUTS
X	/	FD(x1 , ... , xn)

 
Example:               n = 3  STO 00

     x₁ = 0.21   STO 01           a = 0.1   STO 04        b₁ = 0.2   STO 05       c = 2.7   STO 08
     x₂ = 0.31   STO 02                                              b₂ = 0.3    STO 06
     x₃ = 0.41   STO 03                                              b₃ = 0.4    STO 07

  XEQ "LCFD"  >>>>   F_D( 0.21 , 0.31 , 0.41 ) = 1.012345372                ---Execution time = 5m47s---

Note:

-The series is convergent if   Max { | x₁ | ,......, | x_n | }  <  1
 

Reference:

[1]   http://mathworld.wolfram.com/LauricellaFunctions.html