Quadratic Forms for the HP-41

Overview
 

 1°)  Easy Cases
 2°)  General Case
 3°)  Jacobi's Method
 

-A quadratic form is a homogeneous multivariate polynomial of degree 2:

   q(x₁,.....,x_n) = SUM _i<=j  a_i,j x_i x_j

-The following programs diagonalize real quadratic forms.
-In paragraphs 1 & 2, we find coefficients  c_i  and  b_i,j  ( j > i )  such that

   q(x₁,.....,x_n) = SUM _i=1,...,n  c_i ( x_i + b_i,i+1 x_i+1 + ...... + b_i,n x_n )²

-Note that the coefficient of the 1st variable inside the parentheses is always 1, so it is not stored in a register.

  c_i is stored, then b_i,i+1 and so on.
 

1°)  Easy Cases
 

-The problem is relatively easy if all the coefficients of  x_i²  are - more exactly become - different from 0:

-For example,  2 x² + 6 x y + 10 x z + ..... = 2 ( x + 3y/2 + 5z/2 + ... )² - 2 ( 3y/2 )² - 2 ( 5z/2 )² - 2 ( 15yz/2 ) - ...... + ............

-The coefficient of y² is modified and - if it becomes different from zero - we can proceed in the same way with the remaining terms.
 
 

Data Registers:           •  R00 = n < 25                          ( Registers R00 thru Rn(n+1)/2 are to be initialized before executing "QF0" )

                                      •  R01 = a_1,1     •  Rn+1 = a_2,2     ........................    •  Rn(n+1)/2-2 = a_n-1,n-1   •  Rn(n+1)/2 = a_n,n
                                      •  R02 = a_1,2     •  Rn+2 = a_2,3     ........................    •  Rn(n+1)/2-1 = a_n-1,n
                                          ........................................................................

                                      •  Rn-1= a_1,n-1  •  R2n-1 = a_2,n
                                      •  Rnn = a_1,n

  >>> When the program stops,  c_i  and  b_i,j  have replaced  a_i,i  and  a_i,j  ( i < j )

Flag:   F00

   If flag F00 is clear at the end, all worked well !
   If flag F00 is set at the end, there is one - or several - c_i = 0
   and the results will be wrong, except if all the corresponding  b_i,j = 0  too.

Subroutines: /
 
 

 
   ( 121 bytes / SIZE 1+n(n+1)/2 )
 
 

  With   eee = n(n+1)/2

Example1:           q( x,y,z,t ) = x² + 4 x y + 6 x z + 8 x t + 24 y z + 8 y t + 16 z² + 44 z t + 18 t²

  Store these 10 coefficients into registers R01 thru R10

  R01 = 1     R05 =  0     R08 = 16    R10 = 18            ( don't forget to store 0 in R05 since there is no term in y² )
  R02 = 4     R06 = 24    R09 = 44
  R03 = 6     R07 =  8
  R04 = 8

  4   STO 00    ( there are 4 variables )

  XEQ "QF0"   >>>>    1.0100000           ---Execution time = 8s---

  •  Flag F00 is clear  and these registers now contain:

  R01 = 1     R05 =  -4      R08 =  16     R10 = 5
  R02 = 2     R06 = -1.5    R09 = 0.25
  R03 = 3     R07 =   1
  R04 = 4

  So, we have  q(x,y,z,t) = 1 ( x + 2 y + 3 z + 4 t )² - 4 ( y - 3z/2 + t )² + 16 ( z + t/4 )² + 5 t²

-The quadratic form is non-degenerate and its signature is  + + + -
 

Example2:           q( x,y,z,t ) = x² + 4 x y + 6 x z + 8 x t + 4 y² + 24 y z + 8 y t + 16 z² + 44 z t + 18 t²

  Store these 10 coefficients into registers R01 thru R10

  R01 = 1     R05 =  4     R08 = 16    R10 = 18
  R02 = 4     R06 = 24    R09 = 44
  R03 = 6     R07 =  8
  R04 = 8

  4   STO 00    ( 4 variables )

  XEQ "QF0"   >>>>    1.0100000           ---Execution time = 8s---

  •  Flag F00 is set and we have:

  R01 = 1     R05 =   0     R08 =  7     R10 = -86/7
  R02 = 2     R06 =  12    R09 = 10/7
  R03 = 3     R07 =  -8
  R04 = 4

-Since R05 = 0 but  R06 & R07 are different from zero, the results are meaningless !

-But if we re-write q( x,y,z,t ) = 18 t² + 8 t x + 8 t y + 44 t z + x² + 4 x y + 6 x z + 4 y² + 24 y z + 16 z²  and if we store

  R01 = 18     R05 =  1     R08 =  4     R10 = 16
  R02 =  8      R06 =  4     R09 = 24
  R03 =  8      R07 =  6
  R04 = 44

 XEQ "QF0"   >>>>    1.0100000

  •  Flag F00 is clear  and the registers contain:

  R01 = 18      R05 =  1/9      R08 =   -8      R10 = 83/2
  R02 = 2/9     R06 =  10       R09 = -13/4
  R03 = 2/9     R07 = -17
  R04 = 11/9

  So, we have  q(x,y,z,t) = 18 ( t + 2x/9 + 2y/9 + 11z/9 )² + (1/9) ( x + 10 y - 17 z )² - 8 ( y - 13z/4 )² + 83 z² /2

-The quadratic form is non-degenerate and its signature is  + + + -

Example3:           q( x,y,z,t ) = x² + 4 x y + 6 x z + 8 x t + 4 y² + 12 y z + 16 y t + 16 z² + 44 z t + 18 t²

  Store these 10 coefficients into registers R01 thru R10

  R01 = 1     R05 =  4     R08 = 16    R10 = 18
  R02 = 4     R06 = 12    R09 = 44
  R03 = 6     R07 = 16
  R04 = 8

  4   STO 00    ( 4 variables )

  XEQ "QF0"   >>>>    1.0100000           ---Execution time = 8s---

  •  Flag F00 is set   Nevertheless:

  R01 = 1     R05 =  0     R08 =  7     R10 = -86/7
  R02 = 2     R06 =  0     R09 = 10/7
  R03 = 3     R07 =  0
  R04 = 4

  So the result is correct and  q(x,y,z,t) = 1 ( x + 2 y + 3 z + 4 t )² + 7 ( z + 10t/7 )² - 86 t² /7

-This quadratic form is degenerate and its signature is  + + -

Notes:

-For n = 10 , execution time = 87 seconds
-For n = 20 , execution time = 9mn49s
-Since there are at most 319 registers, n cannot exceed 24, but you will probably seldom have to diagonalize a quadratic form over a 25-dimension space...

-The routine is executed in FIX 7 format ( line 06 ) because roundoff-errors may produce a very small number instead of 0
-Change this line according to your preferences, or delete line 06 and set the display setting as you like before executing "QF0"

-As shown in example 2, even if it fails sometimes, the program may work well after changing the order of the variables.
-But of course, it will never work if all the a_i,i = 0 ( unless q = 0 )
-Thus, a more general program must be written:
 

2°)  General Case
 

-More difficult is the case  a_i,i = 0. Suppose the array looks like this, with a # 0

   0    *    *    f    *    *
   0    *    e    g    *
   0    d    *    h
   a    *    *
   b    *
   c

-We must deal with the 1st & 4th columns so that the elements d , e , f , g , h  will be replaced by 0 and if the decomposition is:

   ( x + A y + B z + C t + D u + E v )² - ( x + A y + B z + C' t + D' u + E' v )² + .....        the required conditions lead to:

  A = d/a   B = e/a   and so on if there are similar terms ( more zeros before a # 0 )

  C = f/a + a/4
  C' = f/a - a/4

  D = g/a - bf/a² + b/4         E = h/a - cf/a² + c/4          and so on if there are more terms after a , b , c , ....
  D' = g/a - bf/a² - b/4         E' = h/a - cf/a² - c/4

-Practically, roundoff-errors may produce very small numbers instead of 0 for d , e , ...
 

Data Registers:           •  R00 = n < 25                            ( Registers R00 thru Rn(n+1)/2 are to be initialized before executing "QF" )

                                      •  R01 = a_1,1     •  Rn+1 = a_2,2     ........................    •  Rn(n+1)/2-2 = a_n-1,n-1   •  Rn(n+1)/2 = a_n,n
                                      •  R02 = a_1,2     •  Rn+2 = a_2,3     ........................    •  Rn(n+1)/2-1 = a_n-1,n
                                          ........................................................................

                                      •  Rn-1= a_1,n-1  •  R2n-1 = a_2,n
                                      •  Rnn = a_1,n

  >>> When the program stops,  c_i  and  b_i,j  have replaced  a_i,i  and  a_i,j  ( i < j ) and other registers may also contain such coefficients if needed.

Flags:  /
Subroutines: /

-Lines 27 & 191 are three-byte GTOs
 
 

 
   ( 327 bytes / SIZE 1+n(n+1)/2 + ??? )
 
 

  With   eee >= n(n+1)/2

Example1:           q( x,y,z,t ) = 4 x z + x t + 2 y² + 6 y z + 4 y t + z² + 8 z t - t²

  Store the 10 coefficients into registers R01 thru R10 as before

  R01 = 0     R05 = 2     R08 = 1    R10 = -1
  R02 = 0     R06 = 6     R09 = 8
  R03 = 4     R07 = 4
  R04 = 1

  4   STO 00    ( there are 4 variables )

  XEQ "QF"   >>>>    1.0150000           ---Execution time = 15s---

  •  Registers R01 thru R15 contain:

  R01 =  1        R05 =  2      R08 = 0    R10 = -119/32    R11 = an alpha string   R12 =   -1
  R02 = 3/2      R06 =  0      R09 = 0                                                                   R13 =  3/2
  R03 = 5/4      R07 = 5/8                                                                                    R14 = -3/4
  R04 = 35/16                                                                                                      R15 = 27/16

  So,   q(x,y,z,t) = 1 ( x + 3y/2 + 5z/4 + 35t/16 )² + 2 ( y + 5t/8 )² - 119 t² /32 - 1 ( x + 3y/2 - 3z/4 + 27t/16 )²

-The quadratic form is non-degenerate and its signature is  + + - -

Example2:           q( u,v,x,y,z ) = - 6 u y + 6 u z + 3 v² - 24 v x - 18 v y + 24 v z + 48 x² - 6 x y - 12 x z - 3 y² - 18 y z + 3 z²

  Store the 15 coefficients into registers R01 thru R15

  R01 =  0     R06 =   3      R10 =   48    R13 =  -3    R15 = 3
  R02 =  0     R07 = -24     R11 =  -6     R14 = -18
  R03 =  0     R08 = -18     R12 = -12
  R04 = -6    R09 =   24
  R05 =  6

  5   STO 00    ( there are 5 variables )

  XEQ "QF"   >>>>    1.0250000           ---Execution time = 27s---

  •  Registers R01 thru R25 now contain:

  R01 =  1     R06 =  3      R10 =  1    R13 = 0    R15 = 0     R16 = alpha string   R17 = -1    R22 = alpha string    R23 = -1
  R02 =  3     R07 = -4     R11 =  0     R14 = 0                                                    R18 = 3                                     R24 =  0
  R03 =  1     R08 =  0     R12 = -2                                                                      R19 = 1                                     R25 = -5
  R04 = -1    R09 =  1                                                                                          R20 = 2
  R05 =  5                                                                                                            R21 = 2

 Thus,   q( u,v,x,y,z ) = 1 ( u + 3 v + x - y + 5 z )² + 3 ( v - 4 x + z )² + 1 ( x - 2 z )²  - 1 ( u + 3 v + x + 2 y + 2 z )² - 1 ( x - 5 z )²

-The quadratic form is non-degenerate and its signature is  + + + - -

Notes:

-Since synthetic registers P & Q are used, don't interrupt "QF": their contents would be modified.
-The coefficient that just follows an alpha string is always -1

-If there are m coefficients after an alpha string, the corresponding variables are the last m variables.

-These alpha strings depend on the content of the "alpha" register line 97 and they are only useful to separate the different terms of the sum.
-If there is no alpha string - i-e if eee = n(n+1)/2 - the results are the same as those returned by "QF0".

-Line 12 ( FIX 7 ) may be deleted provided you choose the display format according to the magnitude of the coefficients.
-The M-Code routine X/2 may be used to save bytes & time: replace lines 133 to 138 by  X/2  X/2  to divide X-register by 4 without disturbing Y Z T.
-Likewise,  PI   INT   10^X   /   may be replaced by  X/E3

Variant:

-The coefficients c_i may be inserted inside the parentheses so that the diagonalized quadratic form is expressed by

     q(x₁,.....,x_n) = SUM _i=1,...,n  sign(b_i,i) ( b_i,i x_i + b_i,i+1 x_i+1 + ...... + b_i,n x_n )²     provided we choose the proper sign for  b_i,i

-If you prefer this output,

  Add  X<>Y  after line 142
  Add  CHS  after line 119
  And replace lines 25 thru 33 by

  SIGN          X=0?        SQRT                FS? 30       ST+ X
  STO Z        STO L       *                       GTO 08      ABS
  LASTX      X<> L       STO IND Y       X=0?
  RND          ABS          ISG M               GTO 03                                ----->  209 lines / 340 bytes

Example:       q( x,y,z,t ) = -2 x z - 4 x t + 9 y² - 6 y z - 3 z² - 6 z t + 29 t²

  Store the 10 coefficients into registers R01 thru R10

  R01 =  0     R05 =  9     R08 = -3    R10 = 29
  R02 =  0     R06 = -6     R09 = -6
  R03 = -2     R07 = 0
  R04 = -4

  4   STO 00    ( the 4 variables )

  XEQ "QF"   >>>>    1.0150000           ---Execution time = 15s---

  •  Registers R01 thru R15 contain:

  R01 =  1        R05 = 3      R08 = 0    R10 = 5      R11 = an alpha string       R12 =  -1
  R02 =  3        R06 = 0      R09 = 0                                                              R13 =  -3
  R03 =  1        R07 = 2                                                                                 R14 =  -2
  R04 = -1                                                                                                     R15 =  -1

  So,   q(x,y,z,t) =  ( 1 x + 3 y + z - t )² + ( 3 y + 2 t )² + ( 5 t )² -  ( - x - 3 y - 2 z - t )²

-The signs between the last parentheses may of course be all replaced by +
-However, the formulas involve square-roots which are often irrational numbers...
 

3°)  Jacobi's Method
 

-Another method - though not the fastest - may also be used:
  the method of Jacobi to diagonalize the symmetric matrix M defined by m_i,i = a_i,i & m_i,j = a_i,j/2 for i # j
-It returns the eigenvectors & eigenvalues, which is more than what we really need, but it always work !

-For example, let  q(x,y,z) = 2 x y + 4 x z + 6 y z   We have to diagonalize

    0   1   2
    1   0   3
    2   3   0

-Store these 9 numbers into R01 thru R09
-Store 3 into R00 ( the order of this matrix )
-Press SF 02 and XEQ "JCB"  ( cf "Eigenvalues & Eigenvectors for the HP-41" ), we get in 100 seconds

  R01    R04    R07    R10       4.113090583        0.464141852     -0.842521081      -0.273368927
  R02    R05    R08    R11  =  -0.911178808       0.592851303      0.524794206       -0.610834162
  R03    R06    R09    R12      -3.201911776       0.658103088      0.121446574        0.743068675

-And you can check that

     q(x,y,z) = 4.113090583 ( 0.464141852 x + 0.592851303 y + 0.658103088 z )²
                   - 0.911178808 ( -0.842521081 x + 0.524794206 y + 0.121446574 z )²
                   - 3.201911776 ( -0.273368927 x - 0.610834162 y + 0.743068675 z )² = 2 x y + 4 x z + 6 y z

  with small roundoff-errors as usual...

-Using "JCB2" would be less slow but "QF" returns in 8 seconds:

    q(x,y,z) = ( x + y/2 + 4 z )² - 12 z² - ( x - y/2 + 2 z )²

-So "QF" is by far preferable in most cases !