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

[Valentin Albillo]

 
Hi, all,

As it happens, I've been here for 9 years to the day since I joined back in 2015, and to commemorate that momentous event (plus today it's San Valentin's Day, which also helps,) here you are, a brand new SRC #014 dealing with the accurate computation of NxN real matrix determinants and in particular of those matrices whose elements are all integer, where the determinant's value should mandatorily be integer too but frequently isn't when computing it using our beloved HP calculators' DET functionality.

The Problem

HP vintage models such as the HP-15C, the HP-71B and the RPL models compute the determinant of a square matrix via its LU decomposition, which is a fast, efficient way to do it but involves pivoting and frequently incurs in rounding errors which at times can severely degrade the accuracy of the result, even to the point of rendering it meaningless, with no correct digits whatsoever, and even if the matrix is non-singular. At best, the determinant of all-integer matrices will very frequently be output as a non-integer value.

For example,using MATRIX 9 (HP-15C) or DET (HP-71B) we typically get results like these:

     | 0  3  4 |   Det (15C): 1.999999998
M1 = | 3  1  4 | , Det (71B): 1.99999999995
     | 1  1  2 |   Det Exact: 2

     | 29  18   9 |   Det (15C): -262.0000005
M2 = | 32 -28 -22 | , Det (71B): -262.000000023
     | 18  25  15 |   Det Exact: -262

     | -19  41  22   7 |   Det (15C): -2383.999891
M3 = |   5  19 -14   0 | , Det (71B): -2383.99999934
     | -36  16   9  26 |   Det Exact: -2384
     |  42 -38  14 -38 |

and worst cases include my own Albillo matrices featured in the PDF article HP Article VA016 - Mean Matrices, such as for instance:

       | 58  71  67  36  35  19  60 |
       | 50  71  71  56  45  20  52 |
       | 64  40  84  50  51  43  69 |   Det (15C): 1.080204421
AM#1 = | 31  28  41  54  31  18  33 | , Det (71B): 0.97095056196
       | 45  23  46  38  50  43  50 |   Det Exact: 1
       | 41  10  28  17  33  41  46 |
       | 66  72  71  38  40  27  69 |

       | 13  72  57  94  90  92  35 |
       | 40  93  90  99  01  95  66 |
       | 48  91  71  48  93  32  67 |   Det (15C): -2.956799886
AM#7 = | 07  93  29  02  24  24  07 | , Det (71B): 0.0699243217409 
       | 41  84  44  40  82  27  49 |   Det Exact: 1
       | 03  72  06  33  97  34  04 |
       | 43  82  66  43  83  29  61 |

The Sleuthing

As stated above, the internal LU decomposition and subsequent processing to compute the determinant from it will usually incur in rounding errors, which can accumulate to the point that for difficult (but non-singular) matrices the result will be meaningless or having significantly degraded accuracy, frequently outputting integer determinants as non-integers, even for 2x2 matrices, let alone larger ones.

To attempt solving this annoying problem, some vintage RPL models had a flag setting that would check if all elements were integer and if so it would round the computed result to an integer value. However, this wasn't foolproof at all and at times this ad-hoc rounding would go awry and result in a ±1 error, usually much bigger than just leaving the non-integer value alone. Worse, the user wasn't informed that this rounding had taken place so in the end the cure was worse than the disease.

What to do ? Well, a possible solution would be to use a different algorithm, in particular one which doesn't involve internal arithmetic operations with real numbers, most likely to appear as the intermediate result of non-exact divisions (e.g. 1/3, 293/177, ...) and in 2005 (19 years ago as of 2024 !) I wrote a program to implement this idea, MDET, which you can find featured in this article:

HP Program VA711 - HP-71B Exact Integer Determinants and Permanents

"MDET uses a recursive general expansion by minors procedure and works for any dimensions from 2x2 upwards, though it’s only reasonably efficient for low-order N because computation time grows like N*N!. It produces exact integer results for integer matrices (provided there are no intermediate overflows) [...]"

MDET uses just multiplications and additions/subtractions but no divisions in sight so it certainly delivers the goods and will compute exactly all determnants above. However, its recursive nature means that using it with matrices larger than, say, 7x7 or 10x10 (depending on the calc's speed) would really take too long because the execution time increases super-exponentially.

Thus, again, what to do ? As usual, the idea is correct but the chosen algorithm isn't optimal, we need to use a much faster one, in particular faster than O(N*N!), and this is my implementation for the HP-15C of such an algorithm (the HP-71B version is quite straightforward so I'm not including it here.)

The HP-15C Implementation

This small 50-step (56 bytes = 8 registers) program will accurately compute in polynomial time (not super-exponential time like MDET) the determinant of an arbitrary NxN real matrix A (not necessarily all-integer) for 2 ≤ N ≤ 8 (depending on available memory.)

It takes no inputs (but the user must have previously dimensioned and populated the real square matrix A) and it outputs the computed determinant to the display.

Program listing

   ►LBL A          001- 42,21,11     STO 0          027-    44  0
    RCL DIM A      002- 45,23,11     STO 1          028-    44  1
    STO I          003-    44 25     CLX            029-    43 35
    STO 2          004-    44  2     STO E          030-    44 15
    DSE 2          005- 42, 5, 2     RCL B          031-    45 12
    RCL MATRIX A   006- 45,16,11     CHS            032-       16
    STO MATRIX B   007- 44,16,12     DSE 0          033- 42, 5, 0
   ►LBL 2          008- 42,21, 2    ►LBL 3          034- 42,21, 3
    RCL MATRIX B   009- 45,16,12     RCL 0          035-    45  0
    STO MATRIX E   010- 44,16,15     STO 1          036-    44  1
    RCL I          011-    45 25     X<>Y           037-       34
    STO 0          012-    44  0     STO E          038-    44 15
   ►LBL 0          013- 42,21, 0     RCL- B         039- 45,30,12
    STO 1          014-    44  1     DSE 0          040- 42, 5, 0
    DSE 1          015- 42, 5, 1     GTO 3 ►        041-    22  3
    CLX            016-    43 35     RCL MATRIX E   042- 45,16,15
   ►LBL 1          017- 42,21, 1     RCL MATRIX A   043- 45,16,11
    STO E          018-    44 15     RESULT B       044- 42,26,12
    DSE 1          019- 42, 5, 1      x             045-       20
    GTO 1 ►        020-    22  1     CHS            046-       16
    DSE 0          021- 42, 5, 0     DSE 2          047- 42, 5, 2
     1             022-        1     GTO 2 ►        048-    22  2
    RCL 0          023-    45  0     MATRIX 1       049- 42,16, 1
    TEST 6         024- 43,30, 6     RCL B          050-    45 12
    GTO 0 ►        025-    22  0    
    RCL I          026-    45 25

Notes:

• This stand-alone program works for NxN matrices (2x2 ≤ N ≤ 8x8, depending on the HP-15 physical/virtual model's available memory,) and runs in polynomial time.
   
• It uses matrices A (the input) and auxiliary matrices B and E, all of them NxN, as well as registers R₀-R₂, R_I and labels A, 0-3, but no subroutines or flags.
   
• The input matrix A is not affected by the computation and so remains unaltered and available for further use without having to re-input it or restore its elements back.
   
• The program uses the end of program memory to end execution. To call it instead from another program as a subroutine, add a RTN instruction at its very end (RTN   051- 43 32). It will then return to the calling program with the determinant's value in the X stack register.
   
• The required available memory to compute an NxN determinant is 3*N²+9 registers (program itself included,) so a vintage physical HP-15C/64 can do matrices up to 4x4, the CE/192 can do up to 7x7 and the DM15/M1B can do up to 8x8.

Worked examples

Although all examples below deal with all-integer matrices, the program of course works for arbitrary matrices with non-integer elements and produces their determinants with improved accuracy as well. Assume we're using an HP-15C CE in 192-register mode throughout.


Example 1. Accurately compute the determinant of the following 4x4 all-integer matrix and compare the result with the determinant computed using the buil-in instruction (MATRIX 9):

     | -19  41  22   7 |
M3 = |   5  19 -14   0 |
     | -36  16   9  26 |
     |  42 -38  14 -38 |

- Initialize: allocate memory for register R₂ and clear all matrices to 0x0:

     2, DIM (i), MATRIX 0            (MEM: 02 183 08-0)

- Dimension and populate the input matrix:

     4, ENTER, DIM A, USER, MATRIX 1,

     -19, STO A,  41, STO A,  22, STO A,   7, STO A, 
       5, STO A,  19, STO A, -14, STO A,   0, STO A, 
     -36, STO A,  16, STO A,   9, STO A,  26, STO A, 
      42, STO A, -38, STO A,  14, STO A, -38, STO A

- Compute its determinant:

     FIX 6, GSB A  ->  -2,384.000000

- Compare the result with the one produced by the built-in function (MATRIX 9):
   (no need to re-input the matrix as it's left unaltered by the program)

     RCL MATRIX A, MATRIX 9  ->  -2383.999891

As you can see, the built-in instruction returns a non-integer value with degraded accuracy in the last three places, while the program returns the exact integer value.


Example 2. Ditto for my 7x7 matric AM#1:

       | 50  71  71  56  45  20  52 |
       | 64  40  84  50  51  43  69 |
AM#1 = | 31  28  41  54  31  18  33 |
       | 45  23  46  38  50  43  50 |
       | 41  10  28  17  33  41  46 |
       | 66  72  71  38  40  27  69 |

- Initialize: allocate memory for register R₂ and clear all matrices to 0x0:

     2, DIM (i), MATRIX 0            (MEM: 02 183 08-0)

- Dimension and populate the input matrix:

     7, ENTER, DIM A, USER, MATRIX 1,

     50, STO A, ...,  69, STO A

- Compute its determinant:

     FIX 9, GSB A  ->  1.000000000

- Compare the result with the one produced by the built-in function (MATRIX 9):

     RCL MATRIX A, MATRIX 9  ->  1.080204421

This time the built-in instruction returns a non-integer value with very severely degraded accuracy in the last eight places, while the program again returns the exact integer value.


Example 3. Last but not least, ditto for my 7x7 matric AM#7:

        | 13 72 57 94 90 92 35 |
        | 40 93 90 99 01 95 66 |
        | 48 91 71 48 93 32 67 |
AM #7 = | 07 93 29 02 24 24 07 |
        | 41 84 44 40 82 27 49 |
        | 03 72 06 33 97 34 04 |
        | 43 82 66 43 83 29 61 |

- Initialize: allocate memory for register R₂ and clear all matrices to 0x0:

     2, DIM (i), MATRIX 0            (MEM: 02 183 08-0)

- Dimension and populate the input matrix:

     7, ENTER, DIM A, USER, MATRIX 1,

     13, STO A, ...,  61, STO A

- Compute its determinant:

     FIX 9, GSB A  ->  1.000000000

- Compare the result with the one produced by the built-in function (MATRIX 9):

     RCL MATRIX A, MATRIX 9  ->  -2.956799886

And once more the built-in instruction returns a non-integer value with accuracy so degraded that it loses all 10 digits (and even the sign is wrong !) while the program returns the exact integer value once more.


That'll be all, Over and Out.

V.