# hp41programs

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:

b1 = a1,1 q1 + ................ + a1,n qn

...................................................                          where  bi , qi , ai,j  are quaternions,  qi are the unknowns

bn = an,1 q1 + ................ + an,n qn

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

Data Registers:              R00 = | ProdPiv |2                ( 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 = p1,1               ...................................     •  R21+4n(m-1) to R24+4n(m-1) = p1,m
.....................................................................................................................................

•  R17+4n to R20+4n = pn,1   ....................................   •  R17+4n.m to R20+4n.m = pn,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 |2

where   n = number of rows , m = number of columns
&      | ProdPiv |2  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 |2                       ( Registers R01 thru R4nm are to be initialized before executing "QLS" )

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

•  R4n-3 to R4n = pn,1      ....................................     •  R4n.m-3 to R4n.m = pn,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 |2

where   n = number of rows , m = number of columns
&      | ProdPiv |2  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 |2  , replace lines 138-139 by   SQRT   ST* 00