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
038   READATA
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
038   READATA
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/