Determinants
Ching-Tang Tseng
Hamilton, New Zealand
16 May 2018
In matrix mathematical operation, many calculations are with
respected to evaluate determinant. Determinant of a matrix A can be represented
by det A or det |A|.
Evaluate determinant by computer program is not a simple
procedure, especially for a high rank matrix. There is a Sarrus formula could
be used in evaluating the determinant of a 3x3 matrix. There are no intuitive
formula could be used for a matrix if its rank is higher than 3. No intuitive
formula factor makes this kind of program difficult to design.
According to a mathematical expression, the determinant of
matrix A is expressed as following:
det |A| = ∑ A(1,k) * SIGN * det
(cofactor(k)) = ∑ A(1,k) * (-1)**(k+1) * det (cofactor(k)) –-------- (1)
Range of k is from 1 to rank n.
To design a program based on this equality has an advantage,
determinant evaluating without division operation. But also has disadvantage,
the programming language must owned recursive feature, owing to both sides of
equality listed above has det (...) operation.
FORTH language has recursive feature but inconvenient to
pass too complicated parameters, for example, something like to pass a
sub-matrix det (cofactor(k)).
Fortunately, ancient able man had ever contributed other
evaluating algorithm, for example, a procedure to evaluate a determinant by the
leading diagonal method, using largest pivots. In this method a reduction is
made matrix to upper triangular form by successive row operations. Determinant
then is the product of all elements on the diagonal line. This non-recursive
algorithm must using division operation, results in inevitable inaccuracy.
Programming code must to do some improvements.
On many websites, in many books, we are able to find many
real codes to evaluate determinants. Even though the designer of those codes did
not attach their resource reference, I am sure they are absolutely not the
algorithm inventor. Here, I would like to introduce an oldest resource from my
personal collected books.
P. N. Corlett and J. D. Tinsley, “Practical Programming”
Cambridge University Press, Alden & Mowbray Ltd. Oxford. Great Britain,
1968, P.185
According to this resource, I rewrote the original ALGOL
code into my ABC FORTH system, all side comments included as well. This code
was a resource 50 years ago. I am citing it now. It can be applied to many
places again via ABC FORTH system even after 50 years.
\ A procedure to evaluate a determinant by the leading
diagonal method, using largest pivots. 2008-08-23
\ Lina64 version 20180508
\ ******************************
\ (1).Data structure declared heap
6 integers N I J K P Q
2 reals a d
20 20 (MATRIX) IA
20 20 MATRIX AA
\ ******************************
\ (2).Data input heap
\ det = 160
: INPUT-IDATA
[[ N = 4 ]]
[[ IA ( 1 1 ) = 1 ]] [[ IA ( 1 2 ) =
2 ]] [[ IA ( 1 3 ) = 3 ]] [[ IA ( 1 4 ) = 4 ]]
[[ IA ( 2 1 ) = 2 ]] [[ IA ( 2 2 ) =
3 ]] [[ IA ( 2 3 ) = 4 ]] [[ IA ( 2 4 ) = 1 ]]
[[ IA ( 3 1 ) = 3 ]] [[ IA ( 3 2 ) =
4 ]] [[ IA ( 3 3 ) = 1 ]] [[ IA ( 3 4 ) = 2 ]]
[[ IA ( 4 1 ) = 4 ]] [[ IA ( 4 2 ) =
1 ]] [[ IA ( 4 3 ) = 2 ]] [[ IA ( 4 4 ) = 3 ]]
;
: INPUT-DATA
INPUT-IDATA
basic
10 for i = 1 to N
20 for j = 1 to N
30 let { AA ( i j ) = i>r ( IA ( i j ) ) }
40 next j
50 next i
60 end ;
\ *******************************
\ (3).Auxiliary instructions heap
: SHOW-AA
BASIC
10 run cr
20 FOR I = 1 TO N
30 FOR J = 1 TO N
40 let { a = AA ( I
J ) }
50 run a fs. 3 spaces
60 NEXT J
70 run cr
80 NEXT I
90 END ;
\ ********************************
\ (4).Main code heap
1.0 e -12 fconstant feps
: (DET)
BASIC
\ 10 RUN INPUT-DATA
20 LET { D = f1.0e0 }
30 FOR K = 1 TO ( N - 1 )
40 LET P = K
50 LET Q = K
60 LET { A = ABS ( AA ( K K ) ) }
70 FOR I = K TO N
80 FOR J = K TO N
90 IF { ABS ( AA ( I
J ) ) > A } THEN 110
100 GOTO 140
110 LET P = I
120 LET Q = J
130 LET { A = ABS ( AA ( I J ) ) }
140 NEXT J
150 NEXT I
\ End of search for largest element.
160 IF P <> K THEN 180
170 GOTO 230
180 FOR J = K TO N
190 LET { A = AA ( P J ) }
200 LET { AA ( P J ) = AA ( K J ) }
210 LET { AA ( K J ) = NEGATE ( A ) }
220 NEXT J
\ End of interchange of rows P and K.
230 IF Q <> K THEN 250
240 GOTO 300
250 FOR I = K TO N
260 LET { A = AA ( I Q ) }
270 LET { AA ( I Q ) = AA ( I K ) }
280 LET { AA ( I K ) = NEGATE ( A ) }
290 NEXT I
\ End of interchange of columns Q and K.
\ Largest element is now the pivot.
300 LET { D = D * AA ( K K ) }
310 IF { ABS ( D ) < feps } THEN 430
320 FOR I = K + 1 TO N
330 FOR J = K + 1 TO N
340 LET { A = AA ( K J ) * AA ( I K ) / AA ( K K ) }
350 LET { AA ( I J ) = AA ( I J ) - A }
360 IF { ABS ( AA ( I J ) ) < ABS ( A ) * feps } THEN 380
370 GOTO 390
380 LET { AA ( I J ) = f0.0e0 }
390 NEXT J
400 NEXT I
\ End of reduction to upper triangular form.
410 NEXT K
420 LET { D = D * AA ( N N ) }
430 IF { ABS ( D )
> feps } THEN 450
440 LET { D = f0.0e0
}
450 LET { D = D }
\ Determinant is now in D
\ 460 PRINT { D }
470 END ;
\ *********************************
\ (5).Manipulation heap
: round ( f1 -- f2 )
over abs nlog 12 > if fround then
;
: DET
BASIC
10 RUN INPUT-DATA
20 run cr ." matrix " n . ." X " n .
." : " cr show-aa cr
30 IF N = 1 THEN 50
40 GOTO 70
50 LET { D = AA ( 1 1 ) }
60 GOTO 120
70 IF N = 2 THEN 90
80 GOTO 110
90 LET { D = AA ( 1 1 ) * AA ( 2 2 ) - AA ( 1 2 ) * AA ( 2 1
) }
100 GOTO 120
110 RUN (DET)
120 PRINT { " determinant = " ; D }
130 let { d = d round }
140 run cr ." determinant = " d f. cr
150 END ;
det
\ ********** the end ***************
Code executed outcomes is showing as following:
matrix 4 X 4 :
1.0 E 0 2.0 E
0 3.0 E 0 4.0 E 0
2.0 E 0 3.0 E
0 4.0 E 0 1.0 E 0
3.0 E 0 4.0 E
0 1.0 E 0 2.0 E 0
4.0 E 0 1.0 E
0 2.0 E 0 3.0 E 0
determinant = 1.5999999999 E 2
determinant = 160.0
OK
(4).Main code heap was a typical model program. This code
was my personal collection 10 years ago(2008). Code in each line had been
treated as close to the original ALGOL code as possible. Other heaps were
prepared in recent days to fit to my new ABC FORTH, which was added on Lina64.
Obviously, word “(det)” had been
consolidated to able to exist forever. It is easy to read and to understand.
You are welcome to share this code or to rewrite it into your program language.
I posted many achievements on this website from my personal
record. I am not pursued to promote my ABC FORTH now, and my floating point number
represented format is thoroughly violated the traditional habits. Those points
are not worth to concern, make sure useful of all my design should be
concerned. The code in this article proved this insisted intention.
Integrate code with the answer was listed above for
reference. The method of data input could be redesigned separately. Data input
could be designed into an interactive input style or created by a separated
functional data generator, or input from a separated pure text data file. The
simplest demonstration is as simple as in (2).Data input heap.
At the beginning in (5).Manipulation heap, there was a “round”
instruction designed by “fround” primitive instruction. “fround” is rounding off last one digit of a floating point number,
no matter how many digits of mantissa has. “round” is rounding off last one digit so long as the
digits of mantissa has bigger than 12. this treatment was arranged to solve
accumulated inaccuracy owing to continuous division operation. In my floating
point system, exponent was always to be put on the top of data stack, mantissa
was always to be put on the next secondary cell.
The real functional meaning of FORTH is based on using it ,
is not based on rebuilding it. ABC FORTH could only to be added on a kernel
system, that was not mean I built a new FORTH, I was still using FORTH to
improve its usage, to make FORTH more convenient in computing aspect only.
FORTH original structure had never been changed after added ABC FORTH. Look
inside the design of “round” you can
understand this meaning.
My personal testing
remarks: Instruction “FROUND”
approved.
沒有留言:
張貼留言