Exact results with the HP-71B Message #53 Posted by Valentin Albillo on 28 Apr 2005, 5:03 a.m., in response to message #1 by Valentin Albillo
Hi, all:
Though using DET to compute the Albillo's matrices determinants does result in quite poor results (due to their extreme ill-conditioning and inherent loss of precision caused by the many division operations the internal triangularization algorithm entails), this doesn't mean we can't obtain exact results in a 71B.
In fact, we just need another approach that doesn't involve divisions at any step, and one such approach is to use the "determinant expansion by minors" technique. This algorithm is of considerable theoretical interest, and can be recursively used to compute any determinant using just multiplications and additions, never a division in sight. As such it is highly inefficient for matrices above 4x4, say, but for just the purpose of showing how you can in theory compute those difficult determinants exactly, it will do.
Thanks to the powerful HP-71B's BASIC implementation, a simple version of the expansion by minors can be coded as a recursive subprogram in just 4 lines (with arbitrary line numbers):
100 SUB XDET(A(,),D) @ DIM N,E,I,J,K @ N=UBND(A,1)
110 IF N=2 THEN D=A(1,1)*A(2,2)-A(1,2)*A(2,1) @ END ELSE DIM B(N-1,N-1) @ D=0
120 FOR K=1 TO N @ FOR I=2 TO N @ C=1 @ FOR J=1 TO N @ IF J#K THEN B(I-1,C)=A(I,J) @ C=C+1
130 NEXT J @ NEXT I @ CALL XDET(B,E) @ D=D-(-1)^K*A(1,K)*E @ NEXT K
You must previously dimension the matrix (with OPTION BASE 1) and a real variable to return the value of the determinant, then populate the matrix and call the subprogram, like this:
Albillo Matrix 1:
From the keyboard do:
>OPTION BASE 1 @ DIM A(7,7) @ MAT INPUT A [ENTER]
A(1,1)? 58,71,67,36,35,19,60,50,71,71,56,45,20,52,64,40,84,50,51,43,69 [ENTER]
A(4,1)? 31,28,41,54,31,18,33,45,23,46,38,50,43,50,41,10,28,17,33,41,46 [ENTER]
A(7,1)? 66,72,71,38,40,27,69 [ENTER]
>DET(A) [ENTER]
.97095056196 (inexact result with DET)
>CALL XDET(A,D) @ DISP D [ENTER]
1 (exact result with XDET)
Albillo Matrix 2:
>MAT INPUT A [ENTER]
A(1,1)? 86,69,53,33,18,10,87,70,65,63,55,20,25,73,44,53,63,67,45,44,49 [ENTER]
A(4,1)? 32,27,50,84,36,12,33,28,41,39,64,51,45,33,11,20,13,17,19,33,15 [ENTER]
A(7,1)? 88,73,55,35,21,18,90 [ENTER]
>DET(A) [ENTER]
1.04720740631 (inexact result with DET)
>CALL XDET(A,D) @ DISP D [ENTER]
1 (exact result with XDET)
Albillo Matrix 3:
>MAT INPUT A [ENTER]
A(1,1)? 60,53,61,65,50,37,64,49,78,67,67,24,10,50,31,52,59,60,32,19,33 [ENTER]
A(4,1)? 18,51,50,82,35,27,21,10,18,30,36,38,12,11,6,18,14,33,13,33,10 [ENTER]
A(7,1)? 61,57,63,72,51,45,66 [ENTER]
>DET(A) [ENTER]
.914527613074 (inexact result with DET)
>CALL XDET(A,D) @ DISP D [ENTER]
1 (exact result with XDET)
Please keep in mind that this version of XDET is just a proof-of-concept implementation, to show that it can be done, and as such has no error handling (the matrix must be at least 2x2, for instance), and is extremely inefficient, as it has to recursively call itself (N-1)!/2 times to compute an NxN determinant.
This means that for each of our three examples it calls itself (7-1)!/2 = 360 times, and it will take some 2 min. in Emu71 running on a 2.4 Ghz PC, and one or two hours on a real 71B. It also uses up large amounts of RAM as it deepens the recursion level.
However, it can be easily optimized in a number of ways. For instance, a preliminary search to locate rows/columns containing 0's, then expanding by minors along these rows/columns would mean avoiding a whole sub-branch of recursion for each 0 located. If no 0's are present, they can be created by suitably adding and subtracting rows.
Also, when N is 2, XDET refrains from recursing again to "compute" the required 1x1 determinants, but computes instead the 2x2 determinant directly (halving the number of recursive calls needed). This could be done for N=3 instead, which would result in greater time savings.
Best regards from V.
|