( calc_det4)
( 16th March 2016)

( Calculating a 4x4 determinant by simplifying)
4 FLD !

0 VARIABLE a{n,m}      ( array element)
0 VARIABLE n           ( row number)
0 VARIABLE m           ( column number)
4 CONSTANT #rows    4 CONSTANT #columns


( example of array for order 4 - 4 rows and 4 columns)
#rows #columns ARRAY[] det4

: zero_det4 #rows 0 DO  #columns 0 DO   0 J I det4 !   LOOP  LOOP ;

: dsp_det4 #rows 0 DO CR #columns  0 DO  J I det4 @ F. SPACE LOOP  LOOP CR ;

: load_det4 ( a0,a1,a2,a3...)
   0 #rows 1- DO DUP >R I  det4    ( a0,a1,a2,a3...)
          !  R>                    ( a0,a1,a2,row...)
       -1 +LOOP DROP ;

( values to be floating point)
: init_det4
  zero_det4
  1  FINT  2 FINT -1 FINT  1 FINT   0  load_det4
  -2 FINT  1 FINT  3 FINT -2 FINT   1  load_det4 
  1  FINT  -2 FINT 1 FINT  1 FINT   2  load_det4 
  -3 FINT  1 FINT  2 FINT -1 FINT   3  load_det4 ;

( make column zero be zero apart from a{0,0})
( a{3,n}=a{3,n} - a{1,n}*a{3,0}/a{1,0} - eg for bottom row - #3) 
: nrow_mod  ( row number...) n !   ( save row number)
   n @    m @ det4 @    ( a{3,0}...)
   n @ 1- m @ det4 @    ( a{3,0},a{1,0}...)
   F/                 ( ratio...)
   #columns 0 DO DUP
               n @ 1- I  det4 @    ( ratio,ratio,a{1,n}...) 
               F*  FMNF            ( ratio,-a{1,n}*a{3,0}/a{1,0}...)
               n @    I  det4      ( ratio,-a{1,n}*a{3,0}/a{1,0},addr...)
               F+!                 ( ratio,a{3,n} - a{1,n}*a{3,0}/a{1,0}...)
              LOOP DROP ;

: ncol_mod  ( column no...)  m !
   n @   m @     det4 @    ( a{0,2}...)
   n @   m @ 1-  det4 @    ( a{0,2},a{0,1}...)
   F/                    ( ratio...)
   #columns 0 DO DUP
               I  n @ 1- det4 @
               F* FMNF
               I  n @    det4
               F+!
              LOOP DROP ;


: all_rows ( end row,start row..)
   ( 1 #rows 1-) DO I nrow_mod    CR dsp_det4  -1 +LOOP ;

: all_cols
   ( 1 #columns 1-) DO I ncol_mod   CR dsp_det4   -1 +LOOP ;   

( NB m is column number; n is row number)
: calc4x4
   init_det4
   ( all_rows)
    0 m !  1 #rows 1-     all_rows ( CR ." Here1")
    1 m !  2 #rows 1-     all_rows ( CR ." Here2" KEY DROP)
   ( 0 n !  1 #columns 1-  all_cols CR KEY DROP
    1 n !  2 #columns 1-  all_cols CR KEY DROP)
   2 2 det4 @
   3 3 det4 @  F*
   2 3 det4 @  3 2 det4 @ F*  F-
   1 1 det4 @  F*
   0 0 det4 @  F* 
   CR ." Determinant Value = " F.  ;

  
