Quad Forms

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: /

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 ) = 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" )

• 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

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 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 !

hp41programs

Quadratic Forms for the HP-41