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 = §₀^v ( 1 - m sin² t ) ^-1/2dt    ( 0 < = v < = 90° )
   the angle v = am u   is called the amplitude
   and the 3 Jacobian elliptic functions sn , cn , dn are defined by:

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

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

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

The Arithmetic-Geometric Mean:

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

Starting with a₀ = 1 ; b₀ = (1-m)^1/2 ; c₀ = m^1/2
we determine a₁ , b₁ , c₁ ; ....... ; a_n , b_n , c_n with a_k= (a_k-1+ b_k-1)/2 ; b_k= (a_k-1.b_k-1)^1/2 ; c_k = (a_k-1- b_k-1)/2
and we stop when c_n = 0 ( in fact, when a_n-1= a_n )

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

-When the AGM process stops, the HP-41 calculates   v_n = 2ⁿa_nu*180/pi    ( in degrees )
and then, v_n-1, ... , v₀ from the recurrence relation:    sin (2v_k-1-v_k) = ( c_k/a_k ) sin v_k
then, sn ( u | m ) = sin v₀ ;   cn ( u | m ) = cos v₀   and   dn ( u | m ) = ( 1 - m sin² v₀ )^1/2

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

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

-If m < 0:

Data Registers: R00 = 1 , 2 , ... , n and then n-1, ... , 1 , 0 ; R01 = c₁/a₁ , ..... , Rnn = c_n/a_n
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') ]² and v = u / ( 1+sqrt(µ) ] , we have:

-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(µ₁) , ..... , 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(µ₁) , ..... , 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(µ₁) , ..... , 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/