# hp41programs

JacobiFunct Jacobian Elliptic Functions for the HP-41

Overview

1°)  Via the Arithmetic-Geometric Mean
2°)  Using Gauss' Transformations

a)  Program#1
b)  M-Code Routine ( 0 <= m < 1 )
c)  Program#2

3°)  Jacobian Elliptic Function  Sn ( x | m )  0 <= m < 1

a)  Focal Program
b)  M-Code Routine

1°)  Via the Arithmetic-Geometric Mean

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

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

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

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

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

The Arithmetic-Geometric Mean:

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

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

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

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

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

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

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

-If m < 0:

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

Data Registers:    R00 = 1 , 2 , ... , n  and then   n-1, ... , 1 , 0   ;    R01 = c1/a1  , ..... ,   Rnn = cn/an
Flags:   F01 & F02
Subroutines: /

 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 / SIZE 008 )

 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                           --Execution time = 11s---
RDN      cn ( 0.7 | 0.3 ) = 0.774719736
RDN      dn ( 0.7 | 0.3 ) = 0.938113640

2-Likewise,

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

Notes:

-It seems that n is never greater than 7.
-Therefore, synthetic registers M and N can be replaced by R08 and R09
-If  m < -9999999999 the program can give wrong results!

2°)  Using Gauss' Transformations

a)  Program#1

-If  0 <= m < 1 , this variant employs the following formulae:

-With  m' = 1-m ,  let be  µ = [ ( 1-sqrt(m') / ( 1+sqrt(m') ]2   and   v = u / ( 1+sqrt(µ) ]  ,  we have:

sn ( u | m ) = [ ( 1 + sqrt(µ) ) sn ( v | µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
cn ( u | m ) = [ cn ( v | µ ) dn ( v | µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
dn ( u | m ) = [ 1 - sqrt(µ) sn2 ( v | µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]

-These formulas are applied recursively until µ is small enough to use

sn ( v | 0 ) = Sin v
cn ( v | 0 ) = Cos v
dn ( v | 0 ) = 1

Data Registers:    R00 = 1 , 2 , ... , n  and then   n-1, ... , 1 , 0   ;    R01 = sqrt(µ1) , ..... ,   Rnn = sqrt(µn)
Flags:   F06-F07
Subroutines: /

-Line 41 should be a three-byte GTO 04. It may also be replaced by  RTN  thus saving 2 bytes.

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

( 187 bytes / SIZE 008 )

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

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.632304776                           --Execution time = 7s---
RDN      cn ( 0.7 | 0.3 ) = 0.774719736
RDN      dn ( 0.7 | 0.3 ) = 0.938113640

2-Likewise,

sn ( 0.7 | 1 ) = 0.604367777          sn ( 0.7 | 2 ) = 0.564297007        sn ( 0.7 | -3 ) = 0.759113421
cn ( 0.7 | 1 ) = 0.796705460          cn ( 0.7 | 2 ) = 0.825571855        cn ( 0.7 | -3 ) = 0.650958382
dn ( 0.7 | 1 ) = 0.796705460          dn ( 0.7 | 2 ) = 0.602609138        dn ( 0.7 | -3 ) =1.651895746

Notes:

-It also seems that n is never greater than 7.
-Therefore, synthetic registers M and N can be replaced by R08 and R09
-Storing 1 in synthetic register O and replacing 1 by  RCL O in the loops would reduce execution time.
-If  m < -9999999999 the program can give wrong results!

-If m = 0.9999999999 , execution time = 14 seconds.

b)  M-Code Routine ( 0 <= m < 1 )

-Lines 57 to 111 are the slow part of the program above.
-The following routine "AJF" does the same job much faster.

086   "F"
00A  "J"
001  "A"
378   C=c
03C  RCR 3
106  A=C S&X
1BC  RCR 11
130   LDI S&X
1F8   1F8h=504d                         Correct value for  HP-41 CX , HP-41CV or HP-41C with a quadmemory module
306   ?A<C S&X
381   ?NCGO
00A   NEXIST                             The existence of the required registers is checked
0A6  A<>C S&X
228   P=C
0B8  C=Y
2FE  ?C#0 MS
0B5  ?C GO
0A3  DATAERROR                     DATA ERROR  is displayed if  m < 0
028  T=C
10E  A=C ALL
04E  C
35C  =
050   1
36E  ?A#C ALL
0B5  ?NC GO
0A2  DATAERROR                     DATA ERROR  is displayed if  m = 1
068   Z=C
0F8   C=X
128   L=C
2A0  SETDEC                             LOOP1
0B8  C=Y
2BE  C=-C
10E  A=C ALL
078  C=Z
01D  C=
060   A+C
305  C=                                        If  m > 1 , the HP-41 will display  DATA ERROR
060  sqrt(AB)                                when it tries to calculate sqrt(1-m)
1E8  O=C
2BE  C=-C
10E  A=C ALL
078  C=Z
01D  C=
060   A+C
2EE  ?C#0 ALL
133   JNC+38d                           GOTO  LBL 1
070   N=C ALL
1F8   C=O
10E   A=C ALL
078   C=Z
01D  C=
060   A+C
10E   A=C ALL
0B0   C=N ALL
0AE  A<>C ALL
261   C=
060   A/C
1E8  O=C
001   C=
060   AB+1
10E   A=C ALL
0F8   C=X
0AE  A<>C ALL
261   C=
060   A/C
0E8  X=C
1F8   C=O
10E   A=C ALL
135   C=
060   A*C
0A8  Y=C
1F8   C=O
10E   A=C ALL
238   C=P
260   SETHEX
226   C=C+1 S&X
228   P=C
270   RAMSLCT
0AE  A<>C ALL
2F0   WRITDATA
046   C=0 S&X
270   RAMSLCT
25B   JNC-53d                          GOTO  LOOP1
0F8   C=X                                  LBL 1
231   performs
064   R-D
0A8  Y=C
10E   A=C
078   C=Z
0E8   X=C
070   C&N
351
084   =
1D5
078   P-R(C,A)
0A8  Y=C
0B0   C=N
0E8   X=C
238   C=P
106   A=C S&X                          LOOP2
03C  RCR 3
366   ?A#C S&X
3A0   ?NC RTN                          The routine stops here
078   C=Z
10E   A=C ALL
0B8   C=Y
135   C=
060   A*C
0A8  Y=C
0F8   C=X
10E   A=C ALL
135   C=
060   A*C
238   C=P
270   RAMSLCT
070   N=C ALL
046   C=0 S&X
270   RAMSLCT
0B0   C=N ALL
13D   C=
060   AB*C
2BE  C=-C
1E8  O=C
001  C=
060  AB+1
1A8  N=C
04E  C
35C  =
050   1
268   Q=C
10E   A=C ALL
1F8   C=O
01D   C=
060   A+C
1B8   C=N
269   C=
060   AB/C
068   Z=C
0B8  C=Y
10E  A=C ALL
1B8  C=N
261   C=
060   A/C
0A8  Y=C
278   C=Q
10E   A=C ALL
0B0   C=N ALL
01D  C=
060   A+C
0F8   C=X
13D  C=
060  AB*C
1B8  C=N
269   C=
060   AB/C
0E8   X=C
260   SETHEX
238   C=P
266   C=C-1 S&X
228   P=C
2A0   SETDEC
20B   JNC-63d                                 GOTO  LOOP2

( 163 words / SIZE 008 )

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

with  0 <= m < 1

Example:

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

0.3 ENTER^
0.7 XEQ "AJF"    gives      sn ( 0.7 | 0.3 )  = 0.632304777                           --Execution time = 3.5s---
RDN      cn ( 0.7 | 0.3 ) = 0.774719736
RDN      dn ( 0.7 | 0.3 ) = 0.938113640

Notes:

-If m = 0.9999999999 , execution time = 5.5 seconds instead of 14 seconds.
-The alpha register is not cleared and synthetic register Q is used.
-There is no check for alpha data, underflow or underflow.

c)  Program#2

-The program listed in paragraph 2-a) may now be simplified:

Data Registers:    R00 is unused   ,    R01 = sqrt(µ1) , ..... ,   Rnn = sqrt(µn)
Flags:   F06-F07
Subroutines: /

 01  LBL "JEF"  02  DEG  03  CF 06  04  CF 07  05  X<>Y  06  X#0?  07  X>0?  08  GTO 00         09  CHS  10  STO Z  11  1  12  +  13  ST/ Z  14  SQRT  15  *  16  X<>Y  17  SF 06 18  LBL 00  19  1  20  X>Y?  21  GTO 01  22  X#Y?  23  GTO 00         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 02         42  LBL 00  43  RDN  44  STO Z  45  SQRT  46  *  47  X<>Y  48  1/X  49  ENTER^  50  SF 07  51  LBL 01 52  RDN  53  STO M  54  X<>Y  55  AJF  56  FS?C 06  57  GTO 00         58  FC?C 07  59  GTO 04  60  RDN  61  X<>Y  62  R^  63  RCL M  64  SQRT  65  *  66  GTO 02  67  LBL 00  68  X<> Z 69  ST/ Y  70  ST/ Z  71  1/X  72  X<> Z  73  1  74  ST- M  75  X<> M  76  CHS  77  SQRT  78  *  79  LBL 02         80  CLA  81  END

( 117 bytes / SIZE 008 )

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

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                           --Execution time = 4s---
RDN      cn ( 0.7 | 0.3 ) = 0.774719736
RDN      dn ( 0.7 | 0.3 ) = 0.938113640

2-Likewise,

sn ( 0.7 | 1 ) = 0.604367777          sn ( 0.7 | 2 ) = 0.564297007        sn ( 0.7 | -3 ) = 0.759113420
cn ( 0.7 | 1 ) = 0.796705460          cn ( 0.7 | 2 ) = 0.825571855        cn ( 0.7 | -3 ) = 0.650958382
dn ( 0.7 | 1 ) = 0.796705460          dn ( 0.7 | 2 ) = 0.602609138        dn ( 0.7 | -3 ) =1.651895745

Notes:

-Since R00 is unused by the M-code routine AJF, you can replace synthetic register M by R00 in the listing above.

-If  m < -9999999999 the program can give wrong results !
-If m = 0.9999999999 , execution time = 6 seconds instead of 14 seconds.

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

3°)  Jacobian Elliptic Function  Sn ( x | m )  0 <= m < 1

a)  Focal Program

-If you only have to compute sn ( x | m )  for  0 <= m < 1 - for instance to solve the Ellipsoidal Wave differential equation -
the following program is faster than "JEF"
-Gauss' transformation is used

Data Registers:    R00 = 1 , 2 , ... , n  and then   n-1, ... , 1 , 0   ;    R01 = sqrt(µ1) , ..... ,   Rnn = sqrt(µn)
Flags:  /
Subroutines: /

 01  LBL "SNX"  02  DEG  03  0  04  STO 00          05  RDN  06  LBL 01  07  1  08  RCL Z  09  -  10  SQRT 11  1  12  X<>Y  13  -  14  X=0?  15  GTO 00          16  1  17  ST+ 00  18  LASTX  19  +  20  / 21  STO IND 00  22  X^2  23  X<>Y  24  1  25  LASTX  26  +  27  /  28  GTO 01          29  LBL 00  30  X<>Y 31  R-D  32  SIN  33  LBL 02  34  ENTER^  35  X^2  36  RCL IND 00  37  *  38  1  39  +  40  X<>Y 41  RCL IND 00  42  LASTX  43  +  44  *  45  X<>Y  46  /  47  DSE 00  48  GTO 02  49  END

( 66 bytes / SIZE 008 )

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

Example:            Evaluate  sn ( 0.7 | 0.8 )

0.8  ENTER^
0.7  XEQ "SNX"   >>>>   sn ( 0.7 | 0.8 ) = 0.612365484          ---Execution time = 7s---

Notes:

-There will be an infinte loop if  m = 1
-If m = 0.9999999999 , execution time = 11 seconds instead of 14 seconds.

b)  M-Code Routine

-Likewise, the M-Code routine may be simplified to return sn( u | m ) only:

098  "X"
00E  "N"
013  "S"
378   C=c
03C  RCR 3
106  A=C S&X
1BC  RCR 11
130   LDI S&X
1F8   1F8h=504d                         Correct value for  HP-41 CX , HP-41CV or HP-41C with a quadmemory module
306   ?A<C S&X
381   ?NCGO
00A   NEXIST                             The existence of the required registers is checked
0A6  A<>C S&X
228   P=C
0B8  C=Y
2FE  ?C#0 MS
0B5  ?C GO
0A3  DATAERROR                     DATA ERROR  is displayed if  m < 0
168  M=C
10E  A=C ALL
04E  C
35C  =
050   1
36E  ?A#C ALL
0B5  ?NC GO
0A2  DATAERROR                     DATA ERROR  is displayed if  m = 1
1A8  N=C
0F8   C=X
128   L=C
2A0  SETDEC                             LOOP1
0B8  C=Y
2BE  C=-C
10E  A=C ALL
1B8  C=N
01D  C=
060   A+C
305  C=                                        If  m > 1 , the HP-41 will display  DATA ERROR
060  sqrt(AB)                                when it tries to calculate sqrt(1-m)
1E8  O=C
2BE  C=-C
10E  A=C ALL
1B8  C=N
01D  C=
060   A+C
2EE  ?C#0 ALL
133   JNC+38d                           GOTO  LBL 1
070   N=C ALL
1F8   C=O
10E   A=C ALL
1B8   C=N
01D  C=
060   A+C
10E   A=C ALL
0B0   C=N ALL
0AE  A<>C ALL
261   C=
060   A/C
1E8  O=C
001   C=
060   AB+1
10E   A=C ALL
0F8   C=X
0AE  A<>C ALL
261   C=
060   A/C
0E8  X=C
1F8   C=O
10E   A=C ALL
135   C=
060   A*C
0A8  Y=C
1F8   C=O
10E   A=C ALL
238   C=P
260   SETHEX
226   C=C+1 S&X
228   P=C
270   RAMSLCT
0AE  A<>C ALL
2F0   WRITDATA
046   C=0 S&X
270   RAMSLCT
25B   JNC-53d                          GOTO  LOOP1
130   LDI S&X                           LBL 1
090   090h                                  these instructions
358   ST=C XP                          set the required
144   CLRF 6                             CPU flags
284   CLRF 7                             so that Sin(C)
0F8   C=X                                  is executed in radians
070   N=C ALL                          without changing the current angular mode
229   C=
048   Sin(C)
0E8   X=C
238   C=P
106   A=C S&X                         LOOP2
03C  RCR 3
366   ?A#C S&X
11B   JNC+35d                          GTO LBL2
0F8   C=X
10E   A=C ALL
135   C=
060   A*C
238   C=P
270   RAMSLCT
070   N=C ALL
046   C=0 S&X
270   RAMSLCT
0B0  C=N ALL
13D  C=
060   AB*C
001   C=
060   AB+1
0A8  Y=C
1B8  C=N
10E  A=C ALL
0B0  C=N ALL
01D  C=
060   A+C
0F8   C=X
13D  C=
060   AB*C
0B8  C=N
269   C=
060  AB/C
0E8  X=C
260  SETHEX
238   C=P
266   C=C-1 S&X
228   P=C
2A0  SETDEC
2DB  JNC-37d                        GTO LOOP2
178   C=M                               LBL2
0A8  Y=C
345  ?NCGO
042   CLA

( 136 words / SIZE 008 )

 STACK INPUTS OUTPUTS T T T Z Z Z Y m m X u sn ( u | m ) L L u

with  0 <= m < 1

Example:

-Evaluate  sn ( 0.7 | 0.3 )

0.3 ENTER^
0.7 XEQ "SNX"    gives      sn ( 0.7 | 0.3 )  = 0.632304777                           --Execution time = 2.2s---

Notes:

-If m = 0.8 , execution time = 2.5 seconds.
-If m = 0.9999999999 , execution time = 4 seconds.
-Thus, SNX is significantly faster than AJF.
-Moreover, registers Z & T are saved.

-The alpha register is cleared and synthetic register Q is used.
-There is no check for alpha data, underflow or underflow.

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/