( matrix multiplcation in forth :^]
  make matrix A with '[numbers] columns rows create-mat-a'
  then matrix B with '[numbers] columns rows create-mat-b'
  then do 'multiply' to multiply

  matrix A and B are both 1x1 identity matricies by default

  )

create matrix-a 100 cells allot
create matrix-b 100 cells allot
create matrix-r 100 cells allot

variable a-cols 0 a-cols !
variable a-rows 0 a-rows !
variable b-cols 0 b-cols !
variable b-rows 0 b-rows !
variable r-cols 0 r-cols !
variable r-rows 0 r-rows !

: mat-a   ( col row -- addr ) a-cols @ * + cells matrix-a + ;
: mat-b   ( col row -- addr ) b-cols @ * + cells matrix-b + ;
: mat-r   ( col row -- addr ) r-cols @ * + cells matrix-r + ;

: def-mat-a   ( rows cols -- )
	a-cols ! a-rows ! ;

: def-mat-b   ( cols rows -- )
	b-cols ! b-rows ! ;

: def-mat-r   ( rows col -- )
	r-cols ! r-rows ! ;

: put-mat-a
	cr a-rows @ 0 do ." |" a-cols @ 0 do i j mat-a @ 4 .r loop ." |" cr loop ;
: put-mat-b
	cr b-rows @ 0 do ." |" b-cols @ 0 do i j mat-b @ 4 .r loop ." |" cr loop ;
: put-mat-r
	cr r-rows @ 0 do ." |" r-cols @ 0 do i j mat-r @ 4 .r loop ." |" cr loop ;

: fill-mat-a   ( excepts a-cols a-rows * numbers on stack to fill matrix-a
               eg. for a 3x3 , 1 2 3 4 5 6 fill-mat-a on the stack fills matrix-a like
                   [ 1 2 3 ]
        matrix-a = [ 4 5 6 ]
                   [ 7 8 9 ] )
	0 -1 a-rows @ + do
		0 -1 a-cols @ + do
			i j mat-a !
		-1 +loop
	-1 +loop ;

: fill-mat-b   ( same as fill-mat-a but for matrix-b )
	0 -1 b-rows @ + do
		0 -1 b-cols @ + do
			i j mat-b !
		-1 +loop
	-1 +loop ;

: create-mat-a    ( numbers... cols rows-- ) def-mat-a fill-mat-a ;
: create-mat-b    ( numbers... cols rows-- ) def-mat-b fill-mat-b ;

: can-multiply    ( -- result) a-cols @ b-rows @ = ;
: def-mat-r       ( -- ) a-rows @ r-rows ! b-cols @ r-cols ! ;

: r-get-result    ( col row -- result )
	0 >r \ accumulator
	0 >r \ counter
	begin
		over r@ mat-b @ over r@ swap mat-a @ *
		\ put counter accumulator onto stack
		r> r> swap >r + r> swap >r 1 + dup >r b-rows @ =
	until
	r> r> ( put result on stack ) ;

: r-multiply    ( -- )
	r-rows @ 0 do
		r-cols @ 0 do
			i j r-get-result i j mat-r !
		loop
	loop ;

: multiply
	can-multiply if
		def-mat-R
		r-multiply
		put-mat-a
		." multiplied by  "
		put-mat-b
		." = "
		put-mat-r
	else
		cr ." Can't multiply" cr
	then ;

1 0 0 1 2 2 create-mat-a
1 0 0 1 2 2 create-mat-b