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
 
 

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

  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
 
 

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

  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