(HP71B) integer determinant - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (HP71B) integer determinant (/thread-21302.html) (HP71B) integer determinant - Albert Chan - 02-16-2024 10:36 PM Inspired from [VA] SRC #014 - HP-15C & clones: Accurate NxN Determinants And, tested with examples from thread. HP71B Code translated from @pycoder's Python code see https://en.wikipedia.org/wiki/Bareiss_algorithm 10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N? ";N 20 DIM M(N,N) @ MAT INPUT M @ S=1 @ P=1 30 FOR I=1 TO N-1 @ IF M(I,I) THEN 100 40 FOR J=I+1 TO N @ IF M(J,I) THEN 60 50 NEXT J @ S=0 @ GOTO 140 ! no pivot 60 FOR K=1 TO N @ VARSWAP M(I,K),M(J,K) @ NEXT K @ S=-S 100 P0=P @ P=M(I,I) 110 FOR J=I+1 TO N @ FOR K=I+1 TO N 120 X=P*M(J,K)-M(J,I)*M(I,K) @ M(J,K)=IROUND(X/P0) 130 NEXT K @ NEXT J @ NEXT I 140 DISP "DET =";S*M(N,N) >run N? 3 M(1,1)? 0,3,4 M(2,1)? 3,1,4 M(3,1)? 1,1,2 DET = 2 >run N? 3 M(1,1)? 29,18,9 M(2,1)? 32,-28,-22 M(3,1)? 18,25,15 DET =-262 >run N? 4 M(1,1)? -19,41,22,7 M(2,1)? 5,19,-14,0 M(3,1)? -36,16,9,26 M(4,1)? 42,-38,14,-38 DET =-2384 >run N? 7 M(1,1)? 58,71,67,36,35,19,60 M(2,1)? 50,71,71,56,45,20,52 M(3,1)? 64,40,84,50,51,43,69 M(4,1)? 31,28,41,54,31,18,33 M(5,1)? 45,23,46,38,50,43,50 M(6,1)? 41,10,28,17,33,41,46 M(7,1)? 66,72,71,38,40,27,69 DET = 1 >run N? 7 M(1,1)? 13,72,57,94,90,92,35 M(2,1)? 40,93,90,99,01,95,66 M(3,1)? 48,91,71,48,93,32,67 M(4,1)? 07,93,29,02,24,24,07 M(5,1)? 41,84,44,40,82,27,49 M(6,1)? 03,72,06,33,97,34,04 M(7,1)? 43,82,66,43,83,29,61 DET = 1 RE: (HP71B) integer determinant - Albert Chan - 02-16-2024 11:16 PM OP last example, step by step. (pivot taken from top left corner) Because of delayed divison, all intermediate matrix cells are integers. $$\begin{vmatrix} 13 & 72 & 57 & 94 & 90 & 92 & 35 \\ 40 & 93 & 90 & 99 & 1 & 95 & 66 \\ 48 & 91 & 71 & 48 & 93 & 32 & 67 \\ 7 & 93 & 29 & 2 & 24 & 24 & 7 \\ 41 & 84 & 44 & 40 & 82 & 27 & 49 \\ 3 & 72 & 6 & 33 & 97 & 34 & 4 \\ 43 & 82 & 66 & 43 & 83 & 29 & 61 \end{vmatrix}$$ $$= \begin{vmatrix} -1671 & -1110 & -2473 & -3587 & -2445 & -542 \\ -2273 & -1813 & -3888 & -3111 & -4000 & -809 \\ 705 & -22 & -632 & -318 & -332 & -154 \\ -1860 & -1765 & -3334 & -2624 & -3421 & -798 \\ 720 & -93 & 147 & 991 & 166 & -53 \\ -2030 & -1593 & -3483 & -2791 & -3579 & -712 \end{vmatrix} ÷ (13)^5$$ $$= \begin{vmatrix} 38961 & 67363 & -227290 & 86655 & 9221 \\ 63024 & 215349 & 235401 & 175269 & 49188 \\ 68055 & 74718 & -175932 & 89907 & 25026 \\ 73431 & 118071 & 71283 & 114078 & 36831 \\ 31431 & 61531 & -201373 & 78243 & 6884 \end{vmatrix} ÷ (-1671)^4$$ $$= \begin{vmatrix} -2480387 & -14061151 & -818259 & -799084 \\ 1001377 & -5154838 & 1432938 & -207961 \\ 207282 & -11650143 & 1148157 & -453540 \\ -167578 & 419953 & -194358 & 12937 \end{vmatrix} ÷ (38961)^3$$ $$= \begin{vmatrix} 689574353 & -70194683 & 33777575 \\ 816495643 & -68742161 & 33125188 \\ -87215049 & 8854004 & -4260611 \end{vmatrix} ÷ (-2480387)^2$$ $$= \begin{vmatrix} -3995675528 & 1909767603 \\ 6667765 & -3186916 \end{vmatrix} ÷ (689574353)^1$$ $$= \begin{vmatrix} 1 \\ \end{vmatrix} ÷ (-3995675528)^0 = 1$$ Ref: Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination, by Erwin H. Bareiss RE: (HP71B) integer determinant - Albert Chan - 02-17-2024 05:03 PM It is hard to see why delayed division (by 1 step) preserve integer. Perhaps we can start with smaller 3×3 matrix Cas> m := [[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]] Say we pick a11 as pivot, and we clear pivot column. Cas> m[2] := a11*m[2] - a21*m[1] Cas> m[3] := a11*m[3] - a31*m[1] $$\left(\begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ 0 & a_{11}\; a_{22}- a_{12}\; a_{21} & a_{11}\; a_{23} - a_{13}\; a_{21} \\ 0 & a_{11}\; a_{32}- a_{12}\; a_{31} & a_{11}\; a_{33} - a_{13}\; a_{31} \end{array}\right)$$ This new matrix still have integer entries, with determinant scaled up by pivot^2 If we divide by pivot right the way, entries may have fractional part. If we remove pivot row/col, new matrix have determinant scaled up by pivot^1 Fractional part issue remains, which is easily shown setting pivot = 0 Cas> m := subMat(m, [2,2], [3,3]) Cas> m(a11=0) $$\left(\begin{array}{ccc} - a_{12}\; a_{21} & - a_{13}\; a_{21} \\ - a_{12}\; a_{31} & - a_{13}\; a_{31} \end{array}\right)$$ However, if we delay division by 1 step, cells are divisible by previous pivot. (-a12*a21) * (-a13*a31) - (-a12*a31) * (-a13*a21) = 0 RE: (HP71B) integer determinant - robve - 02-18-2024 03:05 PM Thanks for sharing your observations, which very nicely clarify the algorithmic steps. The one step delay to divide is "math-magic". I wonder if line 30 30 FOR I=1 TO N @ IF M(I,I) THEN 100 should actually run to N-1 and not N, because iteration N is not useful as it skips all zero-trip loops running from I+1 to N? - Rob RE: (HP71B) integer determinant - Albert Chan - 02-18-2024 03:33 PM Hi, robve You are right! Code corrected. All it needed is to reduce from N×N to 1×1 matrix, loops = N-1 RE: (HP71B) integer determinant - J-F Garnier - 02-21-2024 08:23 AM (02-16-2024 10:36 PM)Albert Chan Wrote:  HP71B Code translated from @pycoder's Python code ... 120 X=P*M(J,K)-M(J,I)*M(I,K) @ M(J,K)=IROUND(X/P0) ... Hi Albert, You had to artificially round the X/P0 term to an integer, because the terms P*M(J,K) and M(J,I)*M(I,K) exceed 1e12 at some points (up to 1e15). This seems to be an advantage of the Bird's algorithm used by Valentin in his SRC #014 to have all matrix terms less than 1e10 (for Valentin's examples of course) so the algorithm runs fine on a 10-digit machine such as the 15c. To be more precise, intermediate values during inner products of the matrix multiplication may be up to 1e12, but this is correctly handled by the internal extended accuracy code used during this operation (13 digits on the 15c), and there is no loss of data. Here is my HP-71B implementation of the Bird's algorithm (a translation/adaptation of Valentin's code): 10 ! SRC14B 20 OPTION BASE 1 30 INPUT "N? ";N 40 DIM A(N,N) 50 MAT INPUT A 100 CALL BDET(A,D) 110 PRINT "DET=";D 120 END 130 ! 430 SUB BDET(A(,),D) ! Bird's algorithm 440   N=UBND(A,1) @ DIM B(N,N),E(N,N) 460   MAT B=A 470   FOR K=1 TO N-1 480     MAT E=B 490     FOR I=2 TO N @ FOR J=1 TO I-1 @ E(I,J)=0 @ NEXT J @ NEXT I 520     E(N,N)=0 @ X=0-B(N,N) 530     FOR I=N-1 TO 1 STEP -1 @ E(I,I)=X @ X=X-B(I,I) @ NEXT I 580     MAT E=-E @ MAT B=E*A 590  !  PRINT "B=" @ MAT PRINT B; ! for debug 600   NEXT K 610   D=B(1,1) 620 END SUB The examples I posted in Valentin's thread are failing with the 71B implementation of the Bareiss' algorithm, but are correct processed with this Bird's implementation. J-F Ref: birds-algorithm-for-computing-determinants RE: (HP71B) integer determinant - Albert Chan - 02-21-2024 01:31 PM (02-21-2024 08:23 AM)J-F Garnier Wrote:  The examples I posted in Valentin's thread are failing with the 71B implementation of the Bareiss' algorithm, but are correct processed with this Bird's implementation. Thanks for HP71B Bird's implementation! It work! For fair comparison, I replace dot product using built-in DOT, for extra internal precisions. First, row/column order is swapped (order does not matter, just more efficient the other way) < 120 X=P*M(J,K)-M(J,I)*M(I,K) @ M(J,K)=IROUND(X/P0) > 120 X=P*M(K,J)-M(K,I)*M(I,J) @ M(K,J)=IROUND(X/P0) --> X = [M(I,I), -M(I,J)] • [M(K,J), M(K,I)] = P • Q 10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N? ";N 20 DIM M(N,N),P(2),Q(2) @ MAT INPUT M @ S=1 @ P(1)=1 25 DISP "BUILT-IN DET =";DET(M) 30 FOR I=1 TO N-1 @ IF M(I,I) THEN 100 40 FOR J=I+1 TO N @ IF M(J,I) THEN 60 50 NEXT J @ S=0 @ GOTO 140 ! no pivot 60 FOR K=1 TO N @ VARSWAP M(I,K),M(J,K) @ NEXT K @ S=-S 100 P0=P(1) @ P(1)=M(I,I) 110 FOR J=I+1 TO N @ P(2)=-M(I,J) 115 FOR K=I+1 TO N @ Q(1)=M(K,J) @ Q(2)=M(K,I) 120 X=DOT(P,Q) @ M(K,J)=IROUND(X/P0) 130 NEXT K @ NEXT J @ NEXT I 140 DISP "BAREISS DET =";S*M(N,N) Both of your examples now calculated correctly. >run N? 7 M(1,1)? 274,213,400,322,341,308,446 M(2,1)? 202,210,383,295,360,295,450 M(3,1)? 154,175,332,175,322,361,413 M(4,1)? 176,147,265,272,328,277,351 M(5,1)? 111,131,249,182,324,296,340 M(6,1)? 155,179,329,218,381,399,444 M(7,1)? 147,127,229,174,245,258,297 BUILT-IN DET =-11.4070209786 BAREISS DET = 1 >run N? 7 M(1,1)? 477,389,402,515,358,409,289 M(2,1)? 302,273,282,322,280,283,205 M(3,1)? 278,231,339,319,343,254,214 M(4,1)? 432,360,406,502,391,359,319 M(5,1)? 475,316,509,649,543,393,288 M(6,1)? 299,304,351,369,459,346,221 M(7,1)? 526,561,442,441,371,491,445 BUILT-IN DET =-54.0791415377 BAREISS DET = 1 OP version actually not too bad, getting DET of 3 for 1st, -1 for 2nd We could calculate DET(M MOD 10) MOD 10, to correct for this. (this is a nice check anyway)