[VA] SRC #014 - HP-15C & clones: Accurate NxN Determinants

[Valentin Albillo]

 
Hi, all,

Well, today it's February, 29^th so 15 days have elapsed since I posted my OP and it's high time to add some final comments.

First of all, thanks to all of you for your interest and most especially to those of you who contributed with code and various improvements, namely J-F Garnier and Werner. Frankly, I expected to get more posts since the subject matter is both interesting, instructive and useful, and the HP-15C CE is still very much fashionable but alas, at a meager 7 posts and 1,000 views in 15 days it wasn't to be. Now I'll comment on some of your messages and various matters:

Gene Wrote:Bravo, Valentin!

Thanks a lot for your kind appreciation, Gene.

J-F Garnier Wrote:Another great piece of 15c code from Valentin, and nice to come accross the famous AM matrices again !

Thank you very much for your appreciation and kind words re my code and difficult matrices.

J-F Garnier Wrote:And here we come to the readability of RPN code... I had to first translate the code to a high level language (actually the HP-71B BASIC) to be able to understand it.

You didn't post your 15C/RPN to 71B/BASIC program in this thread so I'm posting my original straightforward HP-71B BASIC version below.

J-F Garnier Wrote:The algorithm (as far as I understand it) is really interesting, not only it provides an integer result for integer input matrices (within certain limits of course), but it doesn't use any division at all, and so doesn't have to choose a non-zero pivot. This simplifies the code a lot.

Indeed. Bird's algorithm doesn't use any divisions at all and that's a big plus, completely eliminating having to care about pivots, having to work with floating point values and having to deal with rounding errors.

J-F Garnier Wrote:The drawback is that it uses 2 auxiliary matrices, which is a real limitation on small machines (such as the original 15C) [...]

Using two extra matrices is unavoidable, as Bird's method essentially relies on matrix multiplication and this usually requires three different matrices on the HP-15C: A=BxC with distinct matrices A, B, C.

J-F Garnier Wrote:This SRC is a great illustration of a recent algorithm, but it works on the 15c (or the 71B) only if the input matrix elements are small enough to avoid any integer overflow in the intermediate calculations.

Yes, the size of intermediate results is a problem for many algorithms, including the naive one which in this case would entail adding 5040 terms (each of them a product of 7 elements,) which are sure to be more than 10-12 digits long if the elements are 3-digit long or more, so Bird's algorithm does no worse than others in this regard and usually it does better.

J-F Garnier Wrote:This SRC achieved its goal, as far as I'm concerned: learn new (even recent) methods for computing determinants, such as the Bareiss and Bird algorithms, and understand (a little bit) how they are working as well as their limitations on our machines.

That's the idea !

J-F Garnier Wrote:Thanks, Valentin, for this SRC !

You're welcome, J-F, again thanks to you for your continued appreciation.

Werner Wrote:The Bird algorithm (thanks, J-F) implemented for the 42S.
I made two changes: - negate the input matrix A once, - shrink the matrix B (remove a row) at every step, as the last row is always filled with zeroes. O, and stack-only. Of course.

Excellent !

Werner Wrote:.. and, my 15C version. - works for n=1, too, - twice as fast, - uses 2n fewer matrix elements, - uses R0 and R1 only (not I), - 45 bytes instead of 57

Most excellent ! The idea of removing a zero row at every step is really clever. Congratulations !

---

Now for my original BASIC program for the HP-71B. It is this 4-line, 191-byte suprogram called BIDET (BIrd's DETerminant, pun most definitely intended):

100  SUB BIDET(A(,),D) @ N=UBND(A,1) @ DIM B(N,N),E(N,N) @ MAT B=A @ FOR K=1 TO N-1 @ MAT E=B
110  FOR I=2 TO N @ FOR J=1 TO I-1 @ E(I,J)=0 @ NEXT J @ NEXT I
120  FOR I=1 TO N @ S=0 @ FOR J=I+1 TO N @ S=S-B(J,J) @ NEXT J @ E(I,I)=S @ NEXT I
130  MAT B=E*A @ MAT B=-B @ NEXT K @ D=B(1,1)

Let's try it with J-F's second 7x7 matrix, namely:

| 477 389 402 515 358 409 289 |
| 302 273 282 322 280 283 205 |
| 278 231 339 319 343 254 214 |
| 432 360 406 502 391 359 319 |
| 475 316 509 649 543 393 288 |
| 299 304 351 369 459 346 221 |
| 526 561 442 441 371 491 445 |

Directly from the command line, first execute the following initialization:

      >DESTROY ALL @ OPTION BASE 1 @ DIM A(7,7),D @ MAT INPUT A

        (... enter all 49 elements ...)

and then let's compute the determinant while also gauging the exactness of the results:

      >CFLAG INX @ CALL BIDET(A,D) @ D,FLAG(INX)

           1             0

so the determinant is computed as 1 and the INeXact flag is 0 (false) so the value is exact.

Now let's compare with the assembly-language DET function:

      >CFLAG INX @ DET(A),FLAG(INX)

          -53.8832026193  1

and this time the INeXact flag is 1 (true) so the computed value is truly inexact. A whole lot !

As for my HP-15C program returning 7269230 instead of 1 (J-F dixit), it can possibly be made to return the exact value by combining it with partitioned matrix techniques, as decribed in my article.

We could try partitions into 4 blocks (from 1x1... to 6x6...) but perhaps the combo 4x4, 4x3, 3x4, 3x3 would possibly avoid intermediate overflows and result in the exact value being produced. This is an area worthy of research for those interested.

Last but not least, the same way that J-F translated my 15C/RPN program to 71B/BASIC, the reverse translation is easy to perform manually line by line and results in translating my 71/BASIC subprogram into my exact 15C/RPN original program, like this:

MAT B=A @ FOR K=1 TO N-1 @ MAT E=B

 ►LBL A
  RCL DIM A      N
  STO I          RI=N
  STO 2
  DSE 2          K=N-1
  RCL MATRIX A
  STO MATRIX B   MAT B=A
 ►LBL 2          FOR K loop
  RCL MATRIX B
  STO MATRIX E   MAT E=B

FOR I=2 TO N @ FOR J=1 TO I-1 @ E(I,J)=0 @ NEXT J @ NEXT I

  RCL I          N
  STO 0          I=N
 ►LBL 0          FOR I loop
  STO 1
  DSE 1          J=I-1
  CLX            0
 ►LBL 1          FOR J loop
  STO E          E(I,J)=0
  DSE 1
  GTO 1►         NEXT J
  DSE 0
   1
  RCL 0
  TEST 6 (X#Y?)
  GTO 0►         NEXT I

S=-B(N,N) @ E(N,N)=0 @ FOR I=N-1 TO 1 STEP -1 @ E(I,I)=S @ S=S-B(I,I) @ NEXT I

  RCL I          N
  STO 0          I=N
  STO 1          J=N
  CLX            0
  STO E          E(N,N)=0
  RCL B          B(N,N)
  CHS            S=-B(N,N)
  DSE 0          I=N-1
 ►LBL 3          FOR I loop
  RCL 0          I
  STO 1          J=I
  X<>Y           S
  STO E          E(I,I)=S
  RCL- B         S=S-B(I,I)
  DSE 0
  GTO 3►         NEXT I
 
MAT B=E*A @ MAT B=-B @ NEXT K @ D=B(1,1)
  
  RCL MATRIX E
  RCL MATRIX A
  RESULT B
   x             MAT B=E*A
  CHS            MAT B=-B
  DSE 2
  GTO 2►         NEXT K
  MATRIX 1       I=1, J=1
  RCL B          B(1,1)
  

Easy, right ?

That's all for now, best regards.

V.