QuatWEF

Quaternionic Weierstrass Elliptic Functions for the HP-41

Overview

1°) Laurent Series
2°) Duplication formula

1°) Laurent Series

-The Weierstrass Elliptic Function P(q;g₂,g₃) may be calculated by a Laurent series:

P(q;g₂;g₃) = q^-2 + c₂.q² + c₃.q⁴ + ...... + c_n.q^2n-2 + ....

where c₂ = g₂/20 ; c₃ = g₃/28 and c_n = 3 ( c₂. c_n-2 + c₃. c_n-3 + ....... + c_n-2. c₂ ) / (( 2n+1 )( n-3 )) ( n > 3 )

-The successive c_n are computed and stored into registers R17 R18 ......

Data Registers: R00 to R16: temp ( Registers R09-R10 are to be initialized before executing "QWEF" )

• R09 = g₂ R17 = c₂
• R10 = g₃ R18 = c₃ R19 = c₄ R20 = c₅ ..........................

>>>> When the program stops, R11 thru R14 = P(q;g₂;g₃)

Flags: /
Subroutines: "Q*Q" "Q^2" "1/Q" ( cf "Quaternions for the HP-41" )

01 LBL "QWEF"
02 XEQ "Q^2"
03 STO 01
04 STO 05
05 RDN
06 STO 02
07 STO 06
08 RDN
09 STO 03
10 STO 07
11 RDN
12 STO 04
13 STO 08
14 RDN
15 XEQ "1/Q"
16 STO 11
17 RDN
18 STO 12
19 RDN
20 STO 13
21 X<>Y
22 STO 14
23 RCL 01
24 RCL 09
25 20
26 /
27 STO 17
28 *
29 ST+ 11
30 RCL 02
31 LASTX
32 *

33 ST+ 12
34 RCL 03
35 LASTX
36 *
37 ST+ 13
38 RCL 04
39 LASTX
40 *
41 ST+ 14
42 RCL 10
43 28
44 /
45 STO 18
46 17
47 STO 00
48 XEQ "Q*Q"
49 STO 05
50 RDN
51 STO 06
52 RCL 18
53 *
54 ST+ 12
55 RDN
56 STO 07
57 LASTX
58 *
59 ST+ 13
60 X<>Y
61 STO 08
62 LASTX
63 *
64 ST+ 14

65 RCL 05
66 LASTX
67 *
68 ST+ 11
69 LBL 10
70 3
71 STO 16
72 LBL 01
73 RCL 00
74 .1
75 %
76 17
77 +
78 STO 15
79 CLX
80 LBL 02
81 RCL IND Y
82 RCL IND 15
83 *
84 +
85 DSE Y
86 ISG 15
87 GTO 02
88 RCL 00
89 1
90 +
91 STO 00
92 14
93 -
94 ENTER^
95 ST+ Y
96 ISG Y

97 CLX
98 3
99 ST* T
100 -
101 *
102 /
103 ISG 15
104 CLX
105 STO IND 15
106 XEQ "Q*Q"
107 STO 05
108 RDN
109 STO 06
110 RCL IND 15
111 *
112 RCL 12
113 +
114 STO 12
115 LASTX
116 -
117 ABS
118 X<>Y
119 STO 07
120 RCL IND 15
121 *
122 RCL 13
123 +
124 STO 13
125 LASTX
126 -
127 ABS
128 +

129 X<>Y
130 STO 08
131 RCL IND 15
132 *
133 RCL 14
134 +
135 STO 14
136 LASTX
137 -
138 ABS
139 +
140 RCL 05
141 RCL IND 15
142 *
143 RCL 11
144 +
145 STO 11
146 LASTX
147 -
148 ABS
149 +
150 X#0?
151 GTO 10
152 DSE 16
153 GTO 01
154 RCL 14
155 RCL 13
156 RCL 12
157 RCL 11
158 END

( 218 bytes / SIZE 020+??? )

STACK	INPUTS	OUTPUTS
T	t	t'
Z	z	z'
Y	y	y'
X	x	x'

where P( x + y i + z j + t k ; g₂ , g₃) = x' + y' i + z' j + t' k

Example: q = 0.2 + 0.3 i + 0.4 j + 0.5 k ; g₂ = 0.9 , g₃ = 1.4

0.9 STO 09
1.4 STO 10

   0.5   ENTER^
   0.4   ENTER^
   0.3   ENTER^
   0.2   XEQ "QWEF"   >>>>   -1.591637275                      ---Execution time = 55s---
                                     RDN    -0.411614082
                                     RDN    -0.548818776
                                     RDN    -0.686023470

P( 0.2 + 0.3 i + 0.4 j + 0.5 k ; 0.9 , 1.4 ) = -1.591637275 - 0.411614082 i - 0.548818776 j - 0.686023470 k

Notes:

-SIZE 029 is enough for this example.
-Lines 107 to 149 may be replaced by

   STOQ     RCL IND 15    ST* T
   5              ST* L              X<> L
   SIGN       ST* Y             DSQ
   CLX         ST* Z             11

-See "M-Code routines for hypercomplex numbers"

2°) Duplication Formula

-The argument q may be "large" and the Laurent series doesn't converge.
-The duplication formula may be used one or several times in this case.

-"WF2Q" takes P(q) in registers X Y Z T and returns P(2q) in the stack and in R01 thru R04.

P(2q) = -2 P(q) + ( 6 P²(q) - g₂/2 )² / ( 4 ( 4 P³(q) - g₂ P(q) - g₃ ) )

Data Registers: R00 to R18: temp ( Registers R09-R10 are to be initialized before executing "WF2Q" )

• R09 = g₂
• R10 = g₃ R11-R12-R13-R14 = P(q)

>>>> When the program stops, R01 thru R04 = P(2q;g₂;g₃)

Flags: /
Subroutines: "Q*Q" "Q^2" "1/Q" ( cf "Quaternions for the HP-41" )

01 LBL "WF2Q"
02 STO 11
03 RDN
04 STO 12
05 RDN
06 STO 13
07 RDN
08 STO 14
09 RDN
10 XEQ "Q^2"
11 STO 01
12 STO 05
13 RDN
14 STO 02
15 STO 06
16 RDN
17 STO 03
18 STO 07
19 X<>Y
20 STO 04
21 STO 08
22 6

23 ST* 01
24 ST* 02
25 ST* 03
26 ST* 04
27 RCL 09
28 2
29 /
30 ST- 01
31 RCL 04
32 RCL 03
33 RCL 02
34 RCL 01
35 XEQ "Q^2"
36 STO 15
37 RDN
38 STO 16
39 RDN
40 STO 17
41 X<>Y
42 STO 18
43 4
44 ST/ 15

45 ST/ 16
46 ST/ 17
47 ST/ 18
48 ST* 05
49 ST* 06
50 ST* 07
51 ST* 08
52 RCL 09
53 ST- 05
54 RCL 11
55 STO 01
56 RCL 12
57 STO 02
58 RCL 13
59 STO 03
60 RCL 14
61 STO 04
62 XEQ "Q*Q"
63 X<> 10
64 ST- 10
65 X<> 10
66 XEQ "1/Q"

67 STO 01
68 RDN
69 STO 02
70 RDN
71 STO 03
72 X<>Y
73 STO 04
74 RCL 15
75 STO 05
76 RCL 16
77 STO 06
78 RCL 17
79 STO 07
80 RCL 18
81 STO 08
82 XEQ "Q*Q"
83 STO 01
84 CLX
85 RCL 11
86 ST+ X
87 ST- 01
88 CLX

89 RCL 12
90 ST+ X
91 -
92 STO 02
93 CLX
94 RCL 13
95 ST+ X
96 -
97 STO 03
98 CLX
99 RCL 14
100 ST+ X
101 -
102 STO 04
103 RCL 03
104 RCL 02
105 RCL 01
106 END

( 163 bytes / SIZE 019 )

STACK	INPUTS	OUTPUTS
T	t	t'
Z	z	z'
Y	y	y'
X	x	x'

where P(q;g₂;g₃) = x + y i + z j + t k and P(2q;g₂;g₃) = x' + y' i + z' j + t' k

Example: We have found P( 0.2 + 0.3 i + 0.4 j + 0.5 k ; 0.9 , 1.4 ) = -1.591637275 - 0.411614082 i - 0.548818776 j - 0.686023470 k

-If the result is still in the stack ( and g₂ , g₃ in R09-R10 ),

simply press XEQ "WF2Q" >>>>   -0.370892104                            ---Execution time = 6s---
                                              RDN    -0.169841921
                                              RDN    -0.226455895
                                              RDN    -0.283069871

-So, P( 0.4 + 0.6 i + 0.8 j + k ; 0.9 , 1.4 ) = -0.370892104 - 0.169841921 i - 0.226455895 j - 0.283069871 k

Notes:

-In fact, the Laurent series still converges for 2 q and a direct evaluation of P(2q) with "QWEF" yields ( in 2m17s )

P( 0.4 + 0.6 i + 0.8 j + k ; 0.9 , 1.4 ) = -0.370892103 - 0.169841920 i - 0.226455894 j - 0.283069870 k

-But SIZE 040 is required here...

Reference:

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4

hp41programs

Quaternionic Weierstrass Elliptic Functions for the HP-41