Quaternionic Linear Systems

Quaternionic Linear Systems for the HP-41

Overview

-This program solves unilateral systems of the form B = A Q or B = Q A
-Gaussian elimination with partial pivoting is used.
-The product of the pivots is also computed:
-There is a unique solution if and only if this product is different from zero.

-Several unilateral systems with the same matrix A may also be solved at the same time.
-Underdetermined & overdetermined systems too.
-Actually, the program is structured like "LS3" ( cf 'Linear and Non-Linear Systems for the HP-41" )

Program Listing

-For 1 system of n equations in n unknowns, re-write it as:

b₁ = a_1,1 q₁ + ................ + a_1,n q_n

................................................... where b_i , q_i , a_i,j are quaternions, q_i are the unknowns

b_n = a_n,1 q₁ + ................ + a_n,n q_n

-More generally, the coefficients form a nxm matrix P = [ p_i,j ] 1 <= i <= n , 1 <= j <= m , which are to be stored in column order from register R21.

Data Registers: R00 = | ProdPiv |² ( Registers R21 thru R20+4nm are to be initialized before executing "QLS" )

R01 to R09 & R16 to R19: temp R20 is unused

>>>> R12 to R15 = Product of the Pivots.

• R21 to R24 = p_1,1 ................................... • R21+4n(m-1) to R24+4n(m-1) = p_1,m
.....................................................................................................................................

• R17+4n to R20+4n = p_n,1 .................................... • R17+4n.m to R20+4n.m = p_n,m

Flags: F00-F01-F02

CF 00 A pivot is regarded as 0 if it's smaller than E-7 ( line 123 )
SF 00 Only 0 is regarded as 0

CF 01 Partial pivoting
SF 01 No Pivoting ( not recommended in general - a few exceptions however )

CF 02 B = A Q is solved the solutions are usually different
SF 02 B = Q A is solved because the product is not commutative.

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

-Lines 145 and 268 are three-byte GTOs

01 LBL "QLS"
02 X<>Y
03 .1
04 %
05 +
06 STO 10
07 ST* Y
08 FRC
09 -
10 STO 11
11 1
12 STO 12
13 LASTX
14 %
15 RCL Z
16 INT
17 +
18 4
19 ST* 10
20 ST* 11
21 *
22 20.02
23 ST+ 11
24 +
25 3
26 -
27 .001004
28 RCL X
29 .004
30 +
31 FS? 02
32 X<>Y
33 STO 19
34 X<>Y
35 STO 18
36 CLX
37 STO 13
38 STO 14
39 STO 15
40 R^
41 LBL 01
42 STO 09
43 FS? 01
44 GTO 00
45 INT
46 RCL 11
47 FRC

48 +
49 STO 00
50 3
51 ST+ 00
52 CLX
53 LBL 02
54 RCL IND 00
55 X^2
56 DSE 00
57 RCL IND 00
58 X^2
59 +
60 DSE 00
61 RCL IND 00
62 X^2
63 +
64 DSE 00
65 RCL IND 00
66 X^2
67 +
68 X<=Y?
69 GTO 02
70 RCL 00
71 X<>Y
72 ENTER^
73 LBL 02
74 RDN
75 DSE 00
76 GTO 02
77 RCL 09
78 ENTER^
79 STO 00
80 FRC
81 R^
82 INT
83 +
84 X=Y?
85 GTO 00
86 STO 01
87 LBL 03
88 4
89 RCL 00
90 RCL 01
91 LBL 04
92 RCL IND X
93 X<> IND Z
94 STO IND Y

95 CLX
96 SIGN
97 ST+ Z
98 +
99 DSE Z
100 GTO 04
101 DSE 00
102 CLX
103 DSE 01
104 GTO 03
105 LASTX
106 CHS
107 ST* 12
108 ST* 13
109 ST* 14
110 ST* 15
111 LBL 00
112 RCL 18
113 RCL 09
114 INT
115 +
116 REGMOVE
117 RCL 19
118 LASTX
119 +
120 REGMOVE
121 CLX
122 FC? 00
123 E-7
124 RCL 01
125 ABS
126 RCL 02
127 ABS
128 +
129 RCL 03
130 ABS
131 +
132 RCL 04
133 ABS
134 +
135 X<=Y?
136 CLX
137 X#0?
138 GTO 00
139 STO 12
140 STO 13
141 STO 14

142 STO 15
143 LBL 00
144 X=0?
145 GTO 00
146 12
147 RCL 18
148 +
149 REGMOVE
150 XEQ "Q*Q"
151 STO 12
152 RDN
153 STO 13
154 RDN
155 STO 14
156 X<>Y
157 STO 15
158 RCL 18
159 5
160 +
161 REGMOVE
162 RCL 04
163 RCL 03
164 RCL 02
165 RCL 01
166 XEQ "1/Q"
167 STO 01
168 RDN
169 STO 02
170 RDN
171 STO 03
172 X<>Y
173 STO 04
174 RCL 18
175 1
176 +
177 REGMOVE
178 RCL 09
179 STO 00
180 LBL 05
181 RCL 00
182 INT
183 RCL 19
184 +
185 REGMOVE
186 XEQ "Q*Q"
187 STO IND 00
188 CLX

189 SIGN
190 ST+ 00
191 X<>Y
192 STO IND 00
193 RDN
194 ST+ 00
195 X<>Y
196 STO IND 00
197 RDN
198 ST+ 00
199 X<>Y
200 STO IND 00
201 3
202 ST- 00
203 DSE 00
204 GTO 05
205 RCL 11
206 X<>Y
207 -
208 STO 00
209 LBL 06
210 RCL 09
211 STO 16
212 ENTER^
213 FRC
214 RCL 00
215 INT
216 +
217 X=Y?
218 GTO 06
219 STO 17
220 INT
221 RCL 18
222 +
223 REGMOVE
224 LBL 07
225 RCL 16
226 INT
227 RCL 19
228 +
229 REGMOVE
230 XEQ "Q*Q"
231 ST- IND 17
232 CLX
233 SIGN
234 ST+ 17
235 X<>Y

236 ST- IND 17
237 RDN
238 ST+ 17
239 X<>Y
240 ST- IND 17
241 RDN
242 ST+ 17
243 X<>Y
244 ST- IND 17
245 DSE 16
246 CLX
247 3
248 ST- 17
249 DSE 17
250 GTO 07
251 LBL 06
252 3
253 ST- 00
254 DSE 00
255 GTO 06
256 LBL 00
257 3
258 ST- 09
259 RCL 10
260 ST- 11
261 RCL 11
262 RCL 09
263 1
264 -
265 X<=Y?
266 CLX
267 DSE X
268 GTO 01
269 RCL 12
270 X^2
271 RCL 13
272 X^2
273 +
274 RCL 14
275 X^2
276 +
277 RCL 15
278 X^2
279 +
280 STO 00
281 END

( 411 bytes / SIZE 21+4n.m )

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	\| ProdPiv \|²

where n = number of rows , m = number of columns
& | ProdPiv |² is the square of the magnitude of a product of all the pivots.

Example1: Find 3 quaternions u , v , w such that

    1 + 2 i + 3 j + 4 k = ( 2 + 3 i + j + 2 k ) u + ( 1 + 3 i + 4 j + 2 k ) v + ( 3 + 4 i + 2 j + k ) w
    2 + i + 4 j + 3 k = ( 1 + 2 i + j + 3 k ) u + ( 2 + i + 3 j + 2 k ) v + ( 2 + 3 i + 4 j + 5 k ) w
    3 + 4 i + 2 j + k = ( 3 + 4 i + 5 j + 6 k ) u + ( 4 + 2 i + 3 j + 4 k ) v + ( 1 + 2 i + 3 j + 4 k ) w

-Store these 48 integers in the same order as above into

R21 R22 R23 R24     R33 R34 R35 R36    R45 R46 R47 R48    R57 R58 R59 R60
R25 R26 R27 R28     R37 R38 R39 R40    R49 R50 R51 R52    R61 R62 R63 R64
R29 R30 R31 R32     R41 R42 R43 R44    R53 R54 R55 R56    R65 R66 R67 R68

CF 00 CF 01 CF 02

-There are 3 rows & 4 columns, so:

3 ENTER^
4 XEQ "QLS" >>>> 19782 = R00 ---Execution time = 78s---

-Since R00 # 0 , there is a unique solution ( u , v , w ) in registers R21 thru R32:

u = -0.398544131 + 0.660903352 i - 0.558993024 j + 0.824992418 k
v = 1.079162875 - 0.798604791 i + 0.672732787 j - 0.722171672 k ( rounded to 9D )
w = 0.462996664 - 0.098422809 i - 0.377767668 j - 0.104185624 k

-"The" product of the pivots, in registers R12 thru R15 is

ProdPiv = -6.096844948 - 66.21646089 i - 123.5648833 j + 9.587928679 k

-This is not really a determinant, except if all the coefficients of the matrix A are complex numbers - in this case R14 = R15 = 0
-Otherwise, you'll get different results if there is other - or no - permutations of the rows !

-Registers R33 to R68 now contain the Identity matrix, with very small numbers instead of zero sometime...

Example2: If you want to solve the similar system ( same coefficients but different order for the products )

    1 + 2 i + 3 j + 4 k = u ( 2 + 3 i + j + 2 k ) + v ( 1 + 3 i + 4 j + 2 k ) + w ( 3 + 4 i + 2 j + k )
    2 + i + 4 j + 3 k = u ( 1 + 2 i + j + 3 k ) + v ( 2 + i + 3 j + 2 k ) + w ( 2 + 3 i + 4 j + 5 k )
    3 + 4 i + 2 j + k = u ( 3 + 4 i + 5 j + 6 k ) + v ( 4 + 2 i + 3 j + 4 k ) + w ( 1 + 2 i + 3 j + 4 k )

-Store the same coefficients into the same registers ( R21 thru R68 )

SF 02

3 ENTER^
4 R/S >>>> 45542 = R00 # 0

-The solution is also unique and we find in registers R21 thru R32

u = -0.204075359 - 0.409029028 i + 0.185499100 j - 0.192569496 k
v = 0.781696017 + 0.824864960 i - 0.006762988 j - 0.118791445 k ( rounded to 9D )
w = 0.170721532 - 0.247441922 i - 0.104189539 j + 0.300623600 k

-And the product of the pivots, in registers R12 thru R15 is now

ProdPiv = -11.66262622 - 103.0309765 i - 113.7225589 j - 147.8437711 k

Notes:

-The execution times above are obtained if on uses M-Code routines Q*Q and 1/Q to compute the multiplication and the inverse of quaternions.
-The focal routines "Q*Q" & "1/Q" are of course a little slower...

-Since calculating the "quasideterminant" is not essential here, we may compute the square of its module only.
-On the other hand, it's more natural to store the coefficients into R01 thru R4n.m and to get the results in the same registers.
-The following variant does that:

Data Registers: R00 = | ProdPiv |² ( Registers R01 thru R4nm are to be initialized before executing "QLS" )

• R01 to R04 = p_1,1 ................................... • R4n(m-1)+1 to R4n(m-1)+4 = p_1,m
........................................................................................................................

• R4n-3 to R4n = p_n,1 .................................... • R4n.m-3 to R4n.m = p_n,m

Flags: F00-F01-F02

CF 00 A pivot is regarded as 0 if it's smaller than E-7 ( line 126 )
SF 00 Only 0 is regarded as 0

CF 01 Partial pivoting
SF 01 No Pivoting ( not recommended in general, with a few exceptions however )

CF 02 B = A Q is solved
SF 02 B = Q A is solved

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

-Lines 145 and 256 are three-byte GTOs

01 LBL "QLS"
02 RCL Y
03 RCL Y
04 *
05 25 E4
06 /
07 1.018
08 +
09 REGMOVE
10 STO 17
11 X<> Z
12 .1
13 %
14 +
15 STO 10
16 ST* Y
17 FRC
18 -
19 STO 11
20 1
21 STO 00
22 LASTX
23 %
24 RCL Z
25 INT
26 +
27 4
28 ST* 10
29 ST* 11
30 *
31 .017
32 ST- 17
33 17
34 ST+ 17
35 +
36 ST+ 11
37 +
38 3
39 -
40 .001004
41 RCL X
42 .004
43 +
44 FS? 02

45 X<>Y
46 STO 13
47 X<>Y
48 STO 12
49 R^
50 LBL 01
51 STO 09
52 FS? 01
53 GTO 00
54 INT
55 RCL 11
56 FRC
57 +
58 STO 14
59 3
60 ST+ 14
61 CLX
62 LBL 02
63 RCL IND 14
64 X^2
65 DSE 14
66 RCL IND 14
67 X^2
68 +
69 DSE 14
70 RCL IND 14
71 X^2
72 +
73 DSE 14
74 RCL IND 14
75 X^2
76 +
77 X<=Y?
78 GTO 02
79 RCL 14
80 X<>Y
81 ENTER^
82 LBL 02
83 RDN
84 DSE 14
85 GTO 02
86 RCL 09
87 ENTER^
88 STO 01

89 FRC
90 R^
91 INT
92 +
93 X=Y?
94 GTO 00
95 STO 02
96 LBL 03
97 4
98 RCL 01
99 RCL 02
100 LBL 04
101 RCL IND X
102 X<> IND Z
103 STO IND Y
104 CLX
105 SIGN
106 ST+ Z
107 +
108 DSE Z
109 GTO 04
110 DSE 01
111 CLX
112 DSE 02
113 GTO 03
114 LBL 00
115 RCL 12
116 RCL 09
117 INT
118 +
119 REGMOVE
120 RCL 13
121 LASTX
122 +
123 REGMOVE
124 CLX
125 FC? 00
126 E-7
127 RCL 01
128 X^2
129 RCL 02
130 X^2
131 +
132 RCL 03

133 X^2
134 +
135 RCL 04
136 X^2
137 +
138 ST* 00
139 SQRT
140 X<=Y?
141 CLX
142 X=0?
143 STO 00
144 X=0?
145 GTO 00
146 RCL 12
147 5
148 +
149 REGMOVE
150 RCL 04
151 RCL 03
152 RCL 02
153 RCL 01
154 XEQ "1/Q"
155 STO 01
156 RDN
157 STO 02
158 RDN
159 STO 03
160 X<>Y
161 STO 04
162 RCL 12
163 1
164 +
165 REGMOVE
166 RCL 09
167 STO 14
168 LBL 05
169 RCL 14
170 INT
171 RCL 13
172 +
173 REGMOVE
174 XEQ "Q*Q"
175 STO IND 14
176 CLX

177 SIGN
178 ST+ 14
179 X<>Y
180 STO IND 14
181 RDN
182 ST+ 14
183 X<>Y
184 STO IND 14
185 RDN
186 ST+ 14
187 X<>Y
188 STO IND 14
189 3
190 ST- 14
191 DSE 14
192 GTO 05
193 RCL 11
194 X<>Y
195 -
196 STO 14
197 LBL 06
198 RCL 09
199 ENTER^
200 STO 15
201 FRC
202 RCL 14
203 INT
204 +
205 X=Y?
206 GTO 06
207 STO 16
208 INT
209 RCL 12
210 +
211 REGMOVE
212 LBL 07
213 RCL 15
214 INT
215 RCL 13
216 +
217 REGMOVE
218 XEQ "Q*Q"
219 ST- IND 16
220 CLX

221 SIGN
222 ST+ 16
223 X<>Y
224 ST- IND 16
225 RDN
226 ST+ 16
227 X<>Y
228 ST- IND 16
229 RDN
230 ST+ 16
231 X<>Y
232 ST- IND 16
233 DSE 15
234 CLX
235 3
236 ST- 16
237 DSE 16
238 GTO 07
239 LBL 06
240 3
241 ST- 14
242 DSE 14
243 GTO 06
244 LBL 00
245 3
246 ST- 09
247 RCL 10
248 ST- 11
249 RCL 11
250 RCL 09
251 1
252 -
253 X<=Y?
254 CLX
255 DSE X
256 GTO 01
257 RCL 17
258 REGMOVE
259 RCL 00
260 END

( 387 bytes / SIZE 4n.m+18 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	\| ProdPiv \|²

where n = number of rows , m = number of columns
& | ProdPiv |² is the square of the magnitude of a product of all the pivots.

-In example 1, the same results are in obtained in 76 seconds instead of 78.

Notes:

-If all the coefficients of the matrix A are integers, X-output = R00 will be an integer too.
-If you prefer to get | ProdPiv | instead of | ProdPiv |² , replace lines 138-139 by SQRT ST* 00

hp41programs

Quaternionic Linear Systems for the HP-41