# 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(x1,.....,xn) = SUM i<=j  ai,j xi xj

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

q(x1,.....,xn) = SUM i=1,...,n  ci ( xi + bi,i+1 xi+1 + ...... + bi,n xn )2

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

ci is stored, then bi,i+1 and so on.

1°)  Easy Cases

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

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

-The coefficient of y2 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 = a1,1     •  Rn+1 = a2,2     ........................    •  Rn(n+1)/2-2 = an-1,n-1   •  Rn(n+1)/2 = an,n
•  R02 = a1,2     •  Rn+2 = a2,3     ........................    •  Rn(n+1)/2-1 = an-1,n
........................................................................

•  Rn-1= a1,n-1  •  R2n-1 = a2,n
•  Rnn = a1,n

>>> When the program stops,  ci  and  bi,j  have replaced  ai,i  and  ai,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 - ci = 0
and the results will be wrong, except if all the corresponding  bi,j = 0  too.

Subroutines: /

 01  LBL "QF0"  02  CF 00  03  CLA  04  RCL 00  05  STO N  06  FIX 7  07  ISG N  08  LBL 01  09  DSE N  10  RCL N  11   E3  12  /  13  ISG X  14  ST+ M  15  RCL IND M 16  ENTER^  17  ISG M  18  FS? 30  19  GTO 05         20  RCL M  21  X<>Y  22  RND  23  X=0?  24  SF 00  25  X=0?  26  STO L  27  X<> L  28  ST+ X  29  X=0?  30  SIGN 31  LBL 02  32  ST/ IND Y  33  ISG Y  34  GTO 02  35  RDN  36  STO O  37  X<>Y  38  LBL 03  39  RCL M  40  RCL IND X  41  X^2  42  RCL Z  43  *  44  ST- IND O  45  CLX 46  RCL IND Y  47  R^  48  *  49  ST+ X  50  SIGN  51  ISG Y  52  X=0?  53  GTO 01  54  ISG O  55  LBL 04  56  CLX  57  RCL IND Y  58  LASTX  59  *  60  ST- IND O 61  ISG O  62  CLX  63  ISG Y  64  GTO 04  65  R^  66  ISG M  67  GTO 03  68  LBL 05  69  RCL M         70  FRC  71  ISG X  72  CLA  73  END

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

 STACK INPUTS OUTPUTS X / 1.eee

With   eee = n(n+1)/2

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

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 y2 )
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 )2 - 4 ( y - 3z/2 + t )2 + 16 ( z + t/4 )2 + 5 t2

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

Example2:           q( x,y,z,t ) = x2 + 4 x y + 6 x z + 8 x t + 4 y2 + 24 y z + 8 y t + 16 z2 + 44 z t + 18 t2

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 t2 + 8 t x + 8 t y + 44 t z + x2 + 4 x y + 6 x z + 4 y2 + 24 y z + 16 z2  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 )2 + (1/9) ( x + 10 y - 17 z )2 - 8 ( y - 13z/4 )2 + 83 z2 /2

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

Example3:           q( x,y,z,t ) = x2 + 4 x y + 6 x z + 8 x t + 4 y2 + 12 y z + 16 y t + 16 z2 + 44 z t + 18 t2

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 )2 + 7 ( z + 10t/7 )2 - 86 t2 /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 ai,i = 0 ( unless q = 0 )
-Thus, a more general program must be written:

2°)  General Case

-More difficult is the case  ai,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 )2 - ( x + A y + B z + C' t + D' u + E' v )2 + .....        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/a2 + b/4         E = h/a - cf/a2 + c/4          and so on if there are more terms after a , b , c , ....
D' = g/a - bf/a2 - b/4         E' = h/a - cf/a2 - 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 = a1,1     •  Rn+1 = a2,2     ........................    •  Rn(n+1)/2-2 = an-1,n-1   •  Rn(n+1)/2 = an,n
•  R02 = a1,2     •  Rn+2 = a2,3     ........................    •  Rn(n+1)/2-1 = an-1,n
........................................................................

•  Rn-1= a1,n-1  •  R2n-1 = a2,n
•  Rnn = a1,n

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

Flags:  /
Subroutines: /

-Lines 27 & 191 are three-byte GTOs

 01  LBL "QF"   02  CLA   03  RCL 00   04  ENTER^   05  STO N   06  1   07  +   08  *   09  2   10  /   11  STO P   12  FIX 7   13  ISG N   14  LBL 01   15  DSE N   16  RCL N   17  PI   18  INT   19  10^X   20  /   21  ISG X   22  ST+ M   23  RCL M   24  RCL IND X   25  ISG M   26  FS? 30   27  GTO 08   28  STO Z   29  RND   30  X=0?   31  GTO 03   32  X<> L   33  ST+ X 34  ISG Y   35  LBL 02   36  ST/ IND Y   37  ISG Y   38  GTO 02   39  RDN   40  STO O   41  X<>Y   42  XEQ 04   43  GTO 01   44  LBL 03   45  CLX   46  RCL IND Y   47  RND   48  X#0?   49  GTO 00   50  STO IND Y   51  ISG Y   52  GTO 03   53  SIGN   54  -   55  STO M   56  GTO 01   57  LBL 04   58  RCL M   59  RCL IND X   60  X^2   61  RCL Z   62  *   63  ST- IND O   64  CLX   65  RCL IND Y   66  R^ 67  *   68  ST+ X   69  SIGN   70  ISG Y   71  X=0?   72  RTN   73  ISG O   74  LBL 05   75  CLX   76  RCL IND Y   77  LASTX   78  *   79  ST- IND O   80  ISG O   81  CLX   82  ISG Y   83  GTO 05   84  R^   85  ISG M   86  GTO 04   87  LBL 00   88  X<> L   89  X<>Y   90  STO O   91  RCL M   92  -   93  STO Q   94  SIGN   95  ST+ P   96  ST+ Q   97  ASTO IND P   98  ST+ P   99  ST- M 100  STO IND M 101  ST+ M 102  CHS 103  STO IND P 104  ST- P 105  RCL N 106  RCL Z 107  RCL O 108  LBL 06 109  DSE Z 110  RCL Z 111  + 112  RCL IND X 113  DSE Q 114  FS? 30 115  GTO 00 116  X<> Z 117  ST/ Z 118  X<> Z 119  STO IND M 120  STO IND P 121  CLX 122  SIGN 123  ST+ M 124  ST+ P 125  RDN 126  GTO 06 127  LBL 00 128  RCL Z 129  / 130  STO O 131  R^ 132  LBL 07 133  SIGN 134  ST* X 135  ST+ X 136  ST+ X 137  ST/ L 138  X<> L 139  + 140  STO IND M 141  LASTX 142  ST+ X 143  - 144  STO IND P 145  CLX 146  SIGN 147  ST+ P 148  + 149  ISG M 150  X=0? 151  GTO 00 152  RCL IND M 153  RCL O 154  * 155  RCL IND Y 156  X<>Y 157  - 158  R^ 159  / 160  RCL IND M 161  GTO 07 162  LBL 00 163  RCL M 164  STO O 165  STO Q 166  RCL N 167  ST- M 168  SIGN 169  ST+ M 170  XEQ 04 171  RCL Q 172  STO O 173  RCL P 174  RCL N         175  - 176  DSE P 177  RCL P 178  PI 179  INT 180  10^X 181  / 182  + 183  X<> M 184  STO Q 185  SIGN 186  ST+ M 187  CHS 188  XEQ 04 189  RCL Q 190  STO M 191  GTO 01 192  LBL 08 193  RCL P 194   E3 195  / 196  ISG X 197  CLA 198  END

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

 STACK INPUTS OUTPUTS X / 1.eee

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

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

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 + 2 ( y + 5t/8 )2 - 119 t2 /32 - 1 ( x + 3y/2 - 3z/4 + 27t/16 )2

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

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

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 )2 + 3 ( v - 4 x + z )2 + 1 ( x - 2 z )2  - 1 ( u + 3 v + x + 2 y + 2 z )2 - 1 ( x - 5 z )2

-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 ci may be inserted inside the parentheses so that the diagonalized quadratic form is expressed by

q(x1,.....,xn) = SUM i=1,...,n  sign(bi,i) ( bi,i xi + bi,i+1 xi+1 + ...... + bi,n xn )2     provided we choose the proper sign for  bi,i

-If you prefer this output,

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 y2 - 6 y z - 3 z2 - 6 z t + 29 t2

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 )2 + ( 3 y + 2 t )2 + ( 5 t )2 -  ( - x - 3 y - 2 z - t )2

-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 mi,i = ai,i & mi,j = ai,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 )2
- 0.911178808 ( -0.842521081 x + 0.524794206 y + 0.121446574 z )2
- 3.201911776 ( -0.273368927 x - 0.610834162 y + 0.743068675 z )2 = 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 )2 - 12 z2 - ( x - y/2 + 2 z )2

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