Threaded Mode | Linear Mode

Werner · 10-12-2023, 11:18 AM

New, shorter version, with comments:
(I also prefix 2-byte instructions with '_' to easily determine the byte count)

I : descriptor
0 : r row index
1 : c col index
2 : i counter
3 : P pivots
RESULT : matrix

207 bytes

  001 LBL B        Factor and Solve
  002 STO I             X Y RESULT
  003 CLX          In:  X Y R
  004 +            Out: R = inv(X).Y
  005 RCL DIM I    PA=LU factorisation
  006 E3
  008 /
  009 STO 2
  010 CLX
  011 STO 3
  012 LBL 1        Outer loop i=1..n
  013 RCL DIM I
  014 RCL 2
  015 INT
  016 -
  017_STOx 3       P := P*(n-i);
  018_ISG 2
  019 RCL 2
  020 STO 0
  021 STO 1
  022 0            max |Ari|, r=i..n
  023 ENTER        X   Y   Z   T
+>024 LBL 2        G   -   p   -
| 025 RCL 0
| 026 R^
| 027 RCL (i)
| 028 ABS
| 029 R^           G |Ari| p   r
| 030 TEST 8       X<Y?
| 031 Rv
| 032 ISG 0
+-033 GTO 2
  034 RCL 2
  035 R^
  036 STO 0
  037 -
  038_STO- 3       Save pivot
  039 RCL DIM I    n  n
  040 STO 1
  041 -            0
  042 TEST 6       X#Y?
  043 GSB 8        Aic :=: Apc; c=n..1
  044 RCL 2
  045 STO 0
  046 STO 1
  047 RCL (i)
  048 ENTER
  049 ISG 0
  050 TEST 6       X#Y? - AFF
  051 GTO 0        end
  052 X=0?
  053 GTO 1
+>054 LBL 3        Ari := Ari/Aii; r=i+1..n
| 055_STO/ (i)
| 056 ISG 0
+-057 GTO 3
  058 RCL 2
  059 STO 0
  060 ISG 0
+>061 LBL 4        Rank-1 update
| 062 RCL 2        for r=i+1 to n
| 063 STO 1          t := Ari;
| 064uRCL (i)        for c=i+1 to n
| 065 ABS              Arc := Arc - t*Aic;
| 066 GSB 9          next c;
| 067 ISG 0        next r;
+-068 GTO 4
  069 GTO 1

  070 LBL C        Solve only
  071 STO I
  072 CLX
  073 +            Set RESULT
  074 LBL 0        Solve LUX = PB
  075 RCL RESULT
  076 X<> I        Work on B
  077 STO RESULT
  078 RCL DIM I
  079 CLX
  080 E3
  082 /
  083 STO 2        i=0.00n
  084 RCL 3
+>085 LBL 5        Z=PB
| 086 ENTER        for i=1 to n
| 087_ISG 2          pi := P/(n-i)!;
| 088 RCL DIM I      pi := INT(pi);
| 089 STO 1          P := P - pi*(n-i)!;
| 090 CLX            r := i + pi;
| 091 RCL 2          Bic :=: Brc; c=m..1
| 092 INT          next i;
| 093 STO 0
| 094 -
| 095 x!
| 096 /
| 097 LASTX
| 098 X<>Y
| 099 INT          pi
| 100_STO+ 0
| 101 x
| 102 -            Pnew Pold
| 103 TEST 6       X#Y?
| 104 GSB 8
| 105 TEST 1       X>0?
+-106 GTO 5
  107 RCL 2
  108 FRAC
  109 1
  110 +
  111 STO 2
+>112 LBL 6        LY=Z
| 113 RCL 2        for i=1 to n-1
| 114 STO 0          for r=i+1 to n
| 115 ISG 0            t := Ari;
| 116 GSB 7            Brc := Brc - t*Bic; c=1..m
| 117_ISG 2          next r;
+-118 GTO 6        next i;
  119 RCL DIM I
  120 X<>Y
  121 STO 2
+>122_LBL .6       UX=Y
| 123 RCL 2        for i=n to 1
| 124 STO 0          Bic := Bic/Aii; c=1..m
| 125 GSB 0          for r=1 to i-1
| 126_LBL .7  <-+      t := Ari;
| 127_STO/ (i)  |      Brc := Brc - t*Bic; c=1..m
| 128 ISG 1     |    next r;
| 129_GTO .7  --+  next i;
| 130 RCL 2
| 131 1
| 132 -
| 133 E3
| 135 /
| 136 STO 0        r=1..i-1
| 137 ISG 0
| 138 GSB 7
| 139_DSE 2
+-140_GTO .6
  141 RCL I
  142 STO RESULT
  143 RTN
                   Subroutines
  144 LBL 7        Rank-1 modification
  145 GSB 0
  146 ABS
  147 GSB 9
  148 ISG 0
  149 GTO 7
  150 RTN

  151 LBL 0        recall Rri
  152 RCL DIM I
  153 E3
  155 /
  156 STO 1
  157 ISG 1        c=1.00m
  158 RCL RESULT
  159 X<> I
  160 RCL 0
  161 RCL 2
  162_RCL g (i)
  163 X<>Y
  164 STO I
  165 X<>Y
  166 RTN

  167 LBL 8        swap subroutine
  168 RCL 2        Aic :=: Arc; DSE c
  169 RCL 1             X Y Z T 0 1 2
  170_RCL g (i)    In:  X       r c i
  171 X<> (i)      Out: X       r   i
  172 RCL 2
  173 RCL 1
  174_STO g (i)
  175 X<>Y
  176 DSE 1
  177 GTO 8
  178 RTN

  179 LBL 9        axpy subroutine
  180 RCL 2        Arc := Arc - t*Aic; ISG c
  181 RCL 1            L X Y Z T 0 1 2
  182_RCL g (i)    In: t         r c i
  183 LASTX
  184 x
  185_STO- (i)
  186 ISG 1
  187 GTO 9
  188 RTN

Cheers, Werner