Perform matrix operations (addition, product, transpose, etc.) in LaTeX?
Here is some code to manipulate matrices of any size. Currently, it can perform additions, subtractions, and multiplication (as well as fetching individual entries, and transposing a matrix, for instance). Entries are floating points that l3fp supports (16 digits of precision).
% Programming-level functions: \fpm_new:N, \fpm_set:Nn, \fpm_gset:Nn,
% \fpm_add:NNN, \fpm_sub:NNN, \fp_neg:NN, \fp_transpose:NN, \fp_mul:NNN.
%
% Expandable programming-level functions: \fpm_lines:N, \fpm_columns:N,
% \fpm_get:Nnn.
%
% Document-level functions: \matnew, \matset, \matgset, \matadd,
% \matsub, \matmul, \mattypeset.
%
\RequirePackage{expl3}
{
\ExplSyntaxOn
%
% Programming-level code, for adding, multiplying, matrices. A matrix
% of size |MxN| is stored as a token list of the form
%
% \s__fpm { M } { N } { {a11} ... {a1N} } ... { {aM1} ... {aMN} } ;
%
% where |\s__fpm| is a marker used to recognize matrices, |M| and |N|
% are non-negative integers, and |aij| are floating point numbers.
%
% (1) Variables.
%
\cs_new_eq:NN \s__fpm \scan_stop: % A marker.
\tl_const:Nn \c_empty_fpm { \s__fpm { 0 } { 0 } ; }
\cs_new_eq:NN \l__fpm_tmpa_fpm \c_empty_fpm
\seq_new:N \l__fpm_lines_seq
\int_new:N \l__fpm_lines_A_int
\int_new:N \l__fpm_lines_B_int
\int_new:N \l__fpm_columns_A_int
\int_new:N \l__fpm_columns_B_int
\tl_new:N \l__fpm_matrix_A_tl
\tl_new:N \l__fpm_matrix_B_tl
\tl_new:N \l__fpm_matrix_C_tl
\seq_new:N \l__fpm_matrix_A_seq
\seq_new:N \l__fpm_matrix_B_seq
\seq_new:N \l__fpm_one_line_A_seq
\seq_new:N \l__fpm_one_line_B_seq
\tl_new:N \l__fpm_one_line_A_tl
\int_new:N \l__fpm_tmpa_int
%
% (3) Storing matrices.
%
\cs_new_protected:Npn \fpm_new:N #1
{ \cs_new_eq:NN #1 \c_empty_fpm }
\cs_new_protected_nopar:Npn \fpm_set:Nn
{ \__fpm_set:NNn \tl_set:Nx }
\cs_new_protected_nopar:Npn \fpm_gset:Nn
{ \__fpm_set:NNn \tl_gset:Nx }
\cs_new_protected:Npn \__fpm_set:NNn #1#2#3
{
\seq_set_split:Nnn \l__fpm_lines_seq { ; } {#3}
\seq_set_filter:NNn \l__fpm_lines_seq \l__fpm_lines_seq
{ ! \tl_if_empty_p:n {##1} }
%
% Now all lines are non-empty.
%
\tl_clear:N \l__fpm_matrix_A_tl
\int_zero:N \l__fpm_lines_A_int
\int_zero:N \l__fpm_columns_A_int
\seq_map_inline:Nn \l__fpm_lines_seq
{
\int_incr:N \l__fpm_lines_A_int
\seq_set_from_clist:Nn \l__fpm_one_line_A_seq {##1}
\int_set:Nn \l__fpm_tmpa_int { \seq_count:N \l__fpm_one_line_A_seq }
\int_compare:nNnT \l__fpm_columns_A_int = \c_zero
{ \int_set_eq:NN \l__fpm_columns_A_int \l__fpm_tmpa_int }
\int_compare:nNnF \l__fpm_tmpa_int = \l__fpm_columns_A_int
{ \seq_map_break:n { \msg_error:nn { fpm } { invalid-size } } }
\tl_put_right:Nx \l__fpm_matrix_A_tl
{ { \seq_map_function:NN \l__fpm_one_line_A_seq \__fpm_set_aux:n } }
}
#1 #2
{
\s__fpm
{ \int_use:N \l__fpm_lines_A_int }
{ \int_use:N \l__fpm_columns_A_int }
\l__fpm_matrix_A_tl
;
}
}
\cs_new:Npn \__fpm_set_aux:n #1 { { \fp_to_tl:n {#1} } }
%
% (4) Extracting the size of a matrix, and its contents.
% |#1| is the matrix, |#2|, |#3| integer variables receiving the
% number of lines and of columns, and |#4| a token list receiving the
% contents of the matrix.
%
\cs_new_protected:Npn \__fpm_get_parts:NNNN #1#2#3#4
{ \exp_after:wN \__fpm_get_parts:NnnwNNN #1 #2 #3 #4 }
\cs_new_protected:Npn \__fpm_get_parts:NnnwNNN \s__fpm #1#2#3 ; #4#5#6
{
\int_set:Nn #4 {#1}
\int_set:Nn #5 {#2}
\tl_set:Nn #6 {#3}
}
%
% (5) Some expandable functions: getting one entry, getting the size.
%
\cs_new:Npn \fpm_lines:N #1
{ \exp_after:wN \__fpm_lines:NnnwN #1 \use_i:nn }
\cs_new:Npn \fpm_columns:N #1
{ \exp_after:wN \__fpm_lines:NnnwN #1 \use_ii:nn }
\cs_new:Npn \__fpm_lines:NnnwN \s__fpm #1#2#3 ; #4 { #4 {#1} {#2} }
\cs_new:Npn \fpm_get:Nnn #1#2#3
{ \exp_after:wN \__fpm_get:Nnnwnn #1 #2 #3 }
\cs_new:Npn \__fpm_get:Nnnwnn \s__fpm #1#2#3 ; #4#5
{ \exp_args:Nf \tl_item:nn { \tl_item:nn {#3} {#4} } {#5} }
%
% (6) Summing matrices
%
\cs_new_protected_nopar:Npn \fpm_add:NNN { \__fpm_add:NNNN + }
\cs_new_protected_nopar:Npn \fpm_sub:NNN { \__fpm_add:NNNN - }
\cs_new_protected:Npn \__fpm_add:NNNN #1#2#3#4
{
\tl_set:Nn \l__fpm_sign_tl {#1}
\__fpm_get_parts:NNNN #3
\l__fpm_lines_A_int \l__fpm_columns_A_int \l__fpm_matrix_A_tl
\__fpm_get_parts:NNNN #4
\l__fpm_lines_B_int \l__fpm_columns_B_int \l__fpm_matrix_B_tl
\int_compare:nNnTF \l__fpm_lines_A_int = \l__fpm_lines_B_int
{
\int_compare:nNnTF \l__fpm_columns_A_int = \l__fpm_columns_B_int
{ \__fpm_add:N #2 }
{ \msg_error:nn { fpm } { invalid-size } }
}
{ \msg_error:nn { fpm } { invalid-size } }
}
\cs_new_protected:Npn \__fpm_add:N #1
{
\seq_set_split:NnV \l__fpm_matrix_A_seq { } \l__fpm_matrix_A_tl
\seq_set_split:NnV \l__fpm_matrix_B_seq { } \l__fpm_matrix_B_tl
\tl_clear:N \l__fpm_matrix_C_tl
\seq_mapthread_function:NNN
\l__fpm_matrix_A_seq
\l__fpm_matrix_B_seq
\__fpm_add_lines:nn
\tl_set:Nx #1
{
\s__fpm
{ \int_use:N \l__fpm_lines_A_int }
{ \int_use:N \l__fpm_columns_A_int }
\l__fpm_matrix_C_tl
;
}
}
\cs_new_protected:Npn \__fpm_add_lines:nn #1#2
{
\seq_set_split:Nnn \l__fpm_one_line_A_seq { } {#1}
\seq_set_split:Nnn \l__fpm_one_line_B_seq { } {#2}
\tl_put_right:Nx \l__fpm_matrix_C_tl
{
{
\seq_mapthread_function:NNN
\l__fpm_one_line_A_seq
\l__fpm_one_line_B_seq
\__fpm_add_entries:nn
}
}
}
\cs_new:Npn \__fpm_add_entries:nn #1#2
{ { \fp_to_tl:n { #1 \l__fpm_sign_tl #2 } } }
%
% (7) Negating all entries.
%
\cs_new_protected:Npn \fpm_neg:NN #1#2
{ \tl_set:Nx #1 { \exp_after:wN \__fpm_neg:Nnnw #2 } }
\cs_new:Npn \__fpm_neg:Nnnw \s__fpm #1#2#3 ;
{ \s__fpm {#1} {#2} \tl_map_function:nN {#3} \__fpm_neg_aux:n ; }
\cs_new:Npn \__fpm_neg_aux:n #1
{ { \tl_map_function:nN {#1} \__fpm_neg_auxii:n } }
\cs_new:Npn \__fpm_neg_auxii:n #1
{ { \fp_to_tl:n { - #1 } } }
%
% (8) Transposing a matrix.
%
\cs_new_protected:Npn \fpm_transpose:NN #1#2
{
\__fpm_get_parts:NNNN #2
\l__fpm_lines_A_int \l__fpm_columns_A_int \l__fpm_matrix_A_tl
\seq_set_split:NnV \l__fpm_matrix_A_seq { } \l__fpm_matrix_A_tl
\tl_clear:N \l__fpm_matrix_B_tl
\prg_replicate:nn { \l__fpm_columns_A_int }
{
\tl_put_right:Nx \l__fpm_matrix_B_tl
{ { \seq_map_function:NN \l__fpm_matrix_A_seq \__fpm_wrap_head:n } }
\seq_set_map:NNn \l__fpm_matrix_A_seq \l__fpm_matrix_A_seq
{ \tl_tail:n {##1} }
}
\tl_set:Nx #1
{
\s__fpm
{ \int_use:N \l__fpm_columns_A_int }
{ \int_use:N \l__fpm_lines_A_int }
\l__fpm_matrix_B_tl
;
}
}
\cs_new:Npn \__fpm_wrap_head:n #1 { { \tl_head:n {#1} } }
%
% (9) Multiplying matrices.
%
\cs_new_protected:Npn \fpm_mul:NNN #1#2#3
{
\int_compare:nNnTF { \fpm_columns:N #2 } = { \fpm_lines:N #3 }
{
\fpm_transpose:NN \l__fpm_tmpa_fpm #3
\__fpm_get_parts:NNNN #2
\l__fpm_lines_A_int \l__fpm_columns_A_int \l__fpm_matrix_A_tl
\__fpm_get_parts:NNNN #3
\l__fpm_lines_B_int \l__fpm_columns_B_int \l__fpm_matrix_B_tl
\tl_set:Nx #1
{
\s__fpm
{ \int_use:N \l__fpm_lines_A_int }
{ \int_use:N \l__fpm_columns_B_int }
\tl_map_function:NN \l__fpm_matrix_A_tl \__fpm_mul_line:n
;
}
}
{ \msg_error:nn { fpm } { invalid-size } }
}
\cs_new:Npn \__fpm_mul_line:n #1
{ { \exp_after:wN \__fpm_mul_line:Nnnwn \l__fpm_tmpa_fpm {#1} } }
\cs_new:Npn \__fpm_mul_line:Nnnwn \s__fpm #1#2#3 ; #4
{ \__fpm_mul_line:nn {#4} #3 \q_recursion_tail \q_recursion_stop }
\cs_new:Npn \__fpm_mul_line:nn #1#2
{
\quark_if_recursion_tail_stop:n {#2}
{
\fp_to_tl:n
{
\__fpm_mul_one:nwn #1 \use_none_delimit_by_q_stop:w
\q_mark #2 \q_nil \q_stop
0
}
}
\__fpm_mul_line:nn {#1}
}
\cs_new:Npn \__fpm_mul_one:nwn #1#2 \q_mark #3
{ #1 * #3 + \__fpm_mul_one:nwn #2 \q_mark }
%
%
% Messages.
%
\msg_new:nnn { fpm } { invalid-size }
{ Sizes~of~matrices~or~lines~don't~match. }
}
\RequirePackage{amsmath, siunitx}
{
\ExplSyntaxOn
%
% Turning matrices into arrays for display.
%
\cs_new_protected:Npn \fpm_to_array:N #1
{
\begin{pmatrix}
\exp_after:wN \__fpm_to_array:Nnnw #1
\end{pmatrix}
}
\cs_new_eq:NN \__fpm_newline: ? % Dummy def.
\cs_new_protected:Npn \__fpm_to_array:Nnnw \s__fpm #1#2#3 ;
{
\cs_gset_nopar:Npn \__fpm_newline:
{ \cs_gset_nopar:Npn \__fpm_newline: { \\ } }
\tl_map_inline:nn {#3}
{
\__fpm_newline:
\seq_set_split:Nnn \l__fpm_one_line_A_seq { } {##1}
\seq_set_map:NNn \l__fpm_one_line_A_seq \l__fpm_one_line_A_seq
{ \__fpm_to_array_entry:n {####1} }
\seq_use:Nnnn \l__fpm_one_line_A_seq { & } { & } { & }
}
}
\cs_new_protected:Npn \__fpm_to_array_entry:n #1
{
\str_case:nnn {#1}
{
{ nan } { \text{nan} }
{ inf } { \infty }
{ -inf } { -\infty }
}
{ \num{#1} }
}
}
\RequirePackage{xparse}
\ExplSyntaxOn
%
% Document-level functions.
%
\NewDocumentCommand { \matnew } { m } { \fpm_new:N #1 }
\NewDocumentCommand { \matset } { mm } { \fpm_set:Nn #1 {#2} }
\NewDocumentCommand { \matgset } { mm } { \fpm_gset:Nn #1 {#2} }
\NewDocumentCommand { \matadd } { mmm } { \fpm_add:NNN #1 #2 #3 }
\NewDocumentCommand { \matsub } { mmm } { \fpm_sub:NNN #1 #2 #3 }
\NewDocumentCommand { \matneg } { mm } { \fpm_neg:NN #1 #2 }
\NewDocumentCommand { \mattranspose } { mm } { \fpm_transpose:NN #1 #2 }
\NewDocumentCommand { \matmul } { mmm } { \fpm_mul:NNN #1 #2 #3 }
\NewDocumentCommand { \mattypeset } { m }
{ \fpm_to_array:N #1 }
\DeclareExpandableDocumentCommand { \matget } { mmm }
{ \fp_to_tl:n { \fpm_get:Nnn #1 {#2} {#3} } }
\ExplSyntaxOff
\documentclass{article}
\begin{document}
\matnew \X
\matnew \Y
\matnew \Z
\matset \X { 1 , 2 + 3 ; 4 , 3.4e22 }
\matset \Y { 3 , 4 ; -5 , 6 }
\begin{align}
\matadd \Z \X \Y
\mattypeset \Z & = \mattypeset \X + \mattypeset \Y \\
\matsub \Z \X \Y
\mattypeset \Z & = \mattypeset \X - \mattypeset \Y \\
\matmul \Z \X \Y
\mattypeset \Z & = \mattypeset \X \times \mattypeset \Y \\
\matmul \Z \Y \X
\mattypeset \Z & = \mattypeset \Y \times \mattypeset \X
\end{align}
\(X[1,2] = \matget\X{1}{2}\).
\end{document}
Edit: added \matget, which extracts an individual entry in the matrix.
calculator package might help.
Matrix operations (muliplications, inverses, determinants) either exactly or with floats (arbitrary precision)
this answer from November 2013 was edited in March 2017 because already in late 2014, removal of a macro from xint had rendered code not usable, and since 1.1 (2014/10/28) xintfrac does not load xinttools automatically.
\documentclass{article}
\usepackage[paperheight=100cm,vscale=0.9]{geometry}
\usepackage{xintfrac}
\RequirePackage{array}
% november 8-11, 2013
\catcode`_ 11
% update 2017/03/23, because some macros stopped being defined by later
% versions of xint...
% A. xintfrac stopped loading xinttools at 1.1 (2014/10/28)
\usepackage{xinttools}
% B. \XINTinFloatSum got removed at 1.1a (2014/11/07)
% \lverb|1.09a: quick write-up, for use by \xintfloatexpr, will need to be
% thought through again. Renamed (and slightly modified) in 1.09h. Should be
% extended for optional precision. Should be rewritten for optimization. |
\def\XINTinFloatSum {\romannumeral0\XINTinfloatsum }%
\def\XINTinfloatsum #1{\expandafter\XINT_floatsum_a\romannumeral-`0#1\relax }%
\def\XINT_floatsum_a #1{\expandafter\XINT_floatsum_b
\romannumeral0\XINTinfloat[\XINTdigits]{#1}\Z }%
\def\XINT_floatsum_b #1\Z #2%
{\expandafter\XINT_floatsum_c\romannumeral-`0#2\Z {#1}\Z}%
\def\XINT_floatsum_c #1%
{\xint_gob_til_relax #1\XINT_floatsum_e\relax\XINT_floatsum_d #1}%
\def\XINT_floatsum_d #1\Z
{\expandafter\XINT_floatsum_b\romannumeral0\XINTinfloatadd {#1}}%
\def\XINT_floatsum_e #1\Z #2\Z { #2}%
% C. \XINT_Abs got removed at 1.2i (2016/12/13)
\def\XINT_Abs #1{\romannumeral0\XINT_abs #1}%
% end of update
\makeatletter
\let\MAT_xintfloatsum\XINTinFloatSum
\newcount\MAT_cnta
\newcount\MAT_cntb
\newcount\MAT_cntc
\newcount\MAT_cntd
\newcount\MAT_cnte
\def\MATset {\def\MAT_xintin {\xintRaw}\MATset_ }%
\def\MATfloatset {\def\MAT_xintin {\XINTinFloat [\XINTdigits]}\MATset_ }%
\def\MATset_ #1#2{%
\def\MATset_name{#1}%
\edef\MAT_tmpa {#2}%
\MAT_cnta \xint_c_ % sets \MAT_cnta to zero
\expandafter\MATset_a
\romannumeral0\expandafter\xintzapspaces\expandafter{\MAT_tmpa};!;%
}%
\def\MATset_a {\futurelet\XINT_token\MATset_b }%
\def\MATset_b #1;{\def\MAT_tmpa{#1}%
\ifx\XINT_token;\expandafter\MATset_w
\else
\ifx\XINT_token!%
\expandafter\expandafter\expandafter\MATset_x
\else
\expandafter\expandafter\expandafter\MATset_c
\fi\fi }%
\def\MATset_w !;{\MATset_x }%
\def\MATset_x {\expandafter\def
\csname MAT@\expandafter\string\MATset_name {I}\expandafter\endcsname
\expandafter {\the\MAT_cnta }%
\expandafter\def
\csname MAT@\expandafter\string\MATset_name {J}\expandafter\endcsname
\expandafter {\the\MAT_cntb }%
\expandafter\edef \MATset_name [##1]%
{\noexpand\csname MAT@\expandafter\string\MATset_name
\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
%
\def\MAT_in #1,#2,{\xint_bye #2\xint_gobble_iv\xint_bye
{\the\numexpr #1}{\the\numexpr #2}\xint_gobble_iii
{\xintZapSpaces{#1}}}%
%
\def\MATset_c {\advance\MAT_cnta \xint_c_i % row count ++
\MAT_cntb \xint_c_ % column count intially zero
\expandafter\MATset_d\romannumeral0\expandafter
\xintzapspaces\expandafter {\MAT_tmpa},!,}%
\def\MATset_d {\futurelet\XINT_token\MATset_e }%
\def\MATset_e #1,{\ifx\XINT_token!\expandafter\MATset_a
\else
\advance\MAT_cntb \xint_c_i
\expandafter\def
\csname MAT@\expandafter\string\MATset_name
{\the\MAT_cnta}{\the\MAT_cntb}\expandafter\endcsname
\expandafter{\romannumeral-`0\MAT_xintin{\xintZapSpacesB{#1}}}%
\expandafter\MATset_d\fi
}%
% \MATdef
\def\MATdef {\def\MAT_xintin {\xintRaw}%
\MATdef_ }%
\def\MATfloatdef {\def\MAT_xintin {\XINTinFloat [\XINTdigits]}%
\MATdef_ }%
% #3 should be a replacement text with #1 and #2 for horizontal and vertical
% indices, which can be expanded to its final result inside an \edef, and this
% result must be parsable by the xint macros.
% WARNING! version of NOV 10 defined only square matrices, this one of NOV 11
% defines *rectangular matrices and has one more argument*
\def\MATdef_ #1#2#3#4{%
\MAT_cnta #2\relax
\MAT_cntb #3\relax
\def\MAT_tmpa ##1##2{#4}%
\MAT_cntc \xint_c_i % =1
\xintloop
{\expandafter\def\expandafter\MAT_tmpc\expandafter
{\expandafter{\the\MAT_cntc}}%
\MAT_cntd \xint_c_i %=1
\xintloop
\expandafter\def\expandafter\MAT_tmpd\expandafter
{\expandafter{\the\MAT_cntd}}%
\edef\MAT_tmpb {\expandafter\expandafter\expandafter\MAT_tmpa
\expandafter\MAT_tmpc\MAT_tmpd}%
\expandafter\def
\csname MAT@\string#1\MAT_tmpc\MAT_tmpd\expandafter\endcsname
\expandafter {\romannumeral-`0\MAT_xintin
{\expandafter\xintZapSpacesB\expandafter{\MAT_tmpb}}}%
\ifnum\MAT_cntd<\MAT_cntb
\advance\MAT_cntd \xint_c_i
\repeat
\ifnum\MAT_cntc<\MAT_cnta
\advance\MAT_cntc \xint_c_i
}\repeat
\expandafter\def
\csname MAT@\string#1{I}\expandafter\endcsname\expandafter {\the\MAT_cnta}%
\expandafter\def
\csname MAT@\string#1{J}\expandafter\endcsname\expandafter {\the\MAT_cntb}%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% \MATsetentry
\def\MATsetentry {\def\MAT_xintin {\xintRaw}%
\MATsetentry_ }%
\def\MATfloatsetentry {\def\MAT_xintin
{\XINTinFloat [\XINTdigits]}%
\MATsetentry_ }%
\def\MATsetentry_ #1#2#3{%
\edef\MAT_tmpa {#3}%
\expandafter\def
\csname MAT@\string#1\MAT_in #2,\xint_bye,\expandafter\endcsname\expandafter
{\romannumeral-`0\MAT_xintin
{\expandafter\xintZapSpaces\expandafter{\MAT_tmpa}}}%
}%
% NOTA BENE
% use of \xintFor is for ease of coding. In an official package, I would use
% special loops for optimal efficiency (the \xintFor is a general tool which has
% safeguards against situations which do not arise here, like groups suddenly
% closing)
% 10 november:
% Current version has already replaced use of \xintFor by \xintloop in a number
% of places, notably for the computation of inverses and determinants.
% but I leave \xintFor in a number of macros.
% Improvements from using less \edef's in various places
\def\MATrelax #1{%
\toks2 \expandafter {\romannumeral-`0\xintSeq {1}{#1[I]}}%
\toks4 \expandafter {\romannumeral-`0\xintSeq {1}{#1[J]}}%
\xintFor* ##1 in {\the\toks2 }
\do{\xintFor* ##2 in {\the\toks4 }
\do{\expandafter\let\csname MAT@\string#1{##1}{##2}\endcsname\relax }}%
\expandafter\let\csname MAT@\string#1{I}\endcsname \relax
\expandafter\let\csname MAT@\string#1{J}\endcsname \relax
\let #1\relax
}%
\def\MATlet #1#2{%
\toks2 \expandafter {\romannumeral-`0\xintSeq {1}{#2[I]}}%
\toks4 \expandafter {\romannumeral-`0\xintSeq {1}{#2[J]}}%
\xintFor* ##1 in {\the\toks2 }
\do{\xintFor* ##2 in {\the\toks4 }
\do{\expandafter\let
\csname MAT@\string#1{##1}{##2}\expandafter\endcsname
\csname MAT@\string#2{##1}{##2}\endcsname
}}%
\expandafter\edef\csname MAT@\string#1{I}\endcsname {#2[I]}%
\expandafter\edef\csname MAT@\string#1{J}\endcsname {#2[J]}%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% \MATapply
% argument #1 is \macro or \macro {arg1}..{argn} where \macro is a macro with
% n+1 arguments.
\def\MATapply #1#2{%
\toks2 \expandafter {\romannumeral-`0\xintSeq {1}{#2[I]}}%
\toks4 \expandafter {\romannumeral-`0\xintSeq {1}{#2[J]}}%
\xintFor* ##1 in {\the\toks2 }
\do{\xintFor* ##2 in {\the\toks4 }
\do
{\toks@ {#1}%
\expandafter\edef
\csname MAT@\string#2{##1}{##2}\expandafter\expandafter\expandafter
\endcsname\expandafter\expandafter\expandafter
{\expandafter\the\expandafter\toks@\expandafter
{\romannumeral-`0\csname MAT@\string#2{##1}{##2}\endcsname }}%
}%
}%
}%
% TRANSPOSE
% Code rewritten to illustrate how one can proceed with \xintloop and counts
% rather than \xintFor.
\def\MATtranspose #1#2{%
\MAT_cnta #2[I]\relax
\MAT_cntb #2[J]\relax
\MAT_cntd \xint_c_i
\xintloop {%
\toks0 \expandafter{\the\MAT_cntd}%
\MAT_cnte \xint_c_i
\xintloop
\toks2 \expandafter{\the\MAT_cnte}%
\expandafter\let
\csname MAT@_tmp{\the\toks2}{\the\toks0}\expandafter\endcsname
\csname MAT@\string#2{\the\toks0}{\the\toks2}\endcsname
\ifnum \MAT_cnte < \MAT_cntb \advance\MAT_cnte \xint_c_i
\repeat
\ifnum \MAT_cntd < \MAT_cnta \advance\MAT_cntd \xint_c_i
}\repeat
\MAT_cntd \xint_c_i
\xintloop {%
\toks0 \expandafter{\the\MAT_cntd}%
\MAT_cnte \xint_c_i
\xintloop
\toks2 \expandafter{\the\MAT_cnte}%
\expandafter\let
\csname MAT@\string#1{\the\toks0}{\the\toks2}\expandafter\endcsname
\csname MAT@_tmp{\the\toks0}{\the\toks2}\endcsname
\ifnum \MAT_cnte < \MAT_cnta \advance\MAT_cnte \xint_c_i
\repeat
\ifnum \MAT_cntd < \MAT_cntb \advance\MAT_cntd \xint_c_i
}\repeat
\expandafter\def\csname MAT@\string#1{I}\expandafter\endcsname
\expandafter {\the\MAT_cntb }%
\expandafter\def\csname MAT@\string#1{J}\expandafter\endcsname
\expandafter {\the\MAT_cnta }%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% SCALAR MULTIPLICATION
\def\MATsmul {\def\MAT_xintin {\xintRaw}%
\def\MAT_MUL {\xintMul}%
\MATsmul_ }%
\def\MATfloatsmul {\def\MAT_xintin {\XINTinFloat [\XINTdigits]}%
\def\MAT_MUL {\XINTinFloatMul}%
\MATsmul_ }%
\def\MATsmul_ #1#2#3{%
\edef\MAT_tmpa {#2}%
\expandafter\def\expandafter\MAT_tmpa\expandafter
{\romannumeral-`0\MAT_xintin
{\expandafter\xintZapSpaces\expandafter{\MAT_tmpa}}}%
\toks0 \expandafter {\romannumeral-`0\xintSeq {1}{#3[I]}}%
\toks2 \expandafter {\romannumeral-`0\xintSeq {1}{#3[J]}}%
\xintFor* ##1 in {\the\toks0 }
\do{\xintFor* ##2 in {\the\toks2 }
\do{\expandafter
\def\csname MAT@\string#1{##1}{##2}\expandafter\endcsname
\expandafter{\romannumeral-`0\MAT_MUL\MAT_tmpa {#3[##1,##2]}}%
}%
}%
\expandafter\edef\csname MAT@\string#1{I}\endcsname {#3[I]}%
\expandafter\edef\csname MAT@\string#1{J}\endcsname {#3[J]}%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% ADDITION
\def\MATadd {\def\MAT_ADD ##1##2{\xintIrr {\xintAdd {##1}{##2}}[0]}%
\MATadd_ }%
\def\MATfloatadd {\def\MAT_ADD {\XINTinFloatAdd }\MATadd_ }%
\def\MATadd_ #1#2#3{%
\edef\MAT_tmpa {\xintSeq {1}{#2[I]}}%
\edef\MAT_tmpb {\xintSeq {1}{#2[J]}}%
\xintFor* ##1 in \MAT_tmpa
\do{\xintFor* ##2 in \MAT_tmpb
\do{\expandafter\def\csname MAT@_tmp{##1}{##2}\expandafter\endcsname
\expandafter{\romannumeral-`0\MAT_ADD {#2[##1,##2]}{#3[##1,##2]}}%
}%
}%
\xintFor* ##1 in \MAT_tmpa
\do{\xintFor* ##2 in \MAT_tmpb
\do{\expandafter\let
\csname MAT@\string#1{##1}{##2}\expandafter\endcsname
\csname MAT@_tmp{##1}{##2}\endcsname
}%
}%
\expandafter\edef\csname MAT@\string#1{I}\endcsname {#2[I]}%
\expandafter\edef\csname MAT@\string#1{J}\endcsname {#2[J]}%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% MULTIPLICATION
\def\MATmul {\def\MAT_MUL {\xintMul }%
\def\MAT_SUM ##1{\xintIrr {\xintSum {##1}}[0]}%
\MATmul_ }%
\def\MATfloatmul {\def\MAT_MUL {\XINTinFloatMul}%
\def\MAT_SUM {\MAT_xintfloatsum}%
\MATmul_ }%
\def\MATmul_ #1#2#3{%
\edef\MAT_tmpa {\xintSeq {1}{#2[I]}}%
\edef\MAT_tmpb {\xintSeq {1}{#3[J]}}%
\edef\MAT_tmpc {\xintSeq {1}{#2[J]}}%
\xintFor* ##1 in \MAT_tmpa
\do{\xintFor* ##2 in \MAT_tmpb
\do{%
\def\MAT_tmpd ####1{\MAT_MUL {#2[##1,####1]}{#3[####1,##2]}}%
\expandafter
\def\csname MAT@_tmp{##1}{##2}\expandafter\endcsname
\expandafter
{\romannumeral-`0\MAT_SUM{\xintApply\MAT_tmpd\MAT_tmpc}}%
}%
}%
\xintFor* ##1 in \MAT_tmpa
\do{\xintFor* ##2 in \MAT_tmpb
\do{\expandafter\let
\csname MAT@\string#1{##1}{##2}\expandafter\endcsname
\csname MAT@_tmp{##1}{##2}\endcsname
}%
}%
\expandafter\edef\csname MAT@\string#1{I}\endcsname {#2[I]}%
\expandafter\edef\csname MAT@\string#1{J}\endcsname {#3[J]}%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% IDENTITY MATRIX
\def\MATid {\def\MAT_tmpf{/1}\MAT_id }%
\def\MATfloatid {\def\MAT_tmpf{}\MAT_id }%
\def\MAT_id #1#2{%
\MAT_cntc #2\relax
\MAT_cnta \xint_c_i % 1
\xintloop
{\expandafter\def\expandafter\MAT_tmpa \expandafter{\the\MAT_cnta}%
\MAT_cntb \xint_c_i % 1
\xintloop
\expandafter\edef
\csname MAT@\string#1{\MAT_tmpa}{\the\MAT_cntb}\endcsname
{\ifnum\MAT_cntb=\MAT_cnta 1\else 0\fi \MAT_tmpf[0]}%
\ifnum\MAT_cntb<\MAT_cntc
\advance\MAT_cntb \xint_c_i
\repeat
\ifnum\MAT_cnta<\MAT_cntc
\advance\MAT_cnta \xint_c_i
}\repeat
\expandafter\def\csname MAT@\string#1{I}\expandafter\endcsname
\expandafter {\the\MAT_cntc}%
\expandafter\def\csname MAT@\string#1{J}\expandafter\endcsname
\expandafter {\the\MAT_cntc}%
\edef #1[##1]%
{\noexpand\csname
MAT@\string#1\noexpand\MAT_in ##1,\noexpand\xint_bye,\endcsname }%
}%
% INVERSES AND DETERMINANTS
\def\MATinv {\def\MATinvordet_Ia{\MATinv_Ia}%
\def\MATinvordet_II{\MATinv_II}%
\MATinvordet }
\def\MATdet {\def\MAT_det {1/1[0]}% initial value
\def\MATinvordet_Ia{\MATdet_Ia}%
\def\MATinvordet_II{\edef\MAT_det{\xintIrr{\MAT_det}[0]}%
\MATdet_end}%
\MATinvordet }
\def\MATfloatinv {\def\MATinvordet_Ia{\MATinv_Ia}%
\def\MATinvordet_II{\MATinv_II}%
\MATfloatinvordet }
\def\MATfloatdet {\def\MAT_det {1[0]}% initial value
\def\MATinvordet_Ia{\MATdet_Ia}%
\def\MATinvordet_II{\edef\MAT_det{\xintFloat{\MAT_det}}\MATdet_end}%
\MATfloatinvordet }
\def\MATandinverse #1{\def\MAT_name {#1}\MATinv_II }%
\def\MATinvordet #1#2{%
\def\MAT_ZERO {0/1[0]}%
\def\MAT_DIV ##1##2{\xintIrr{\xintDiv {##1}{##2}}}%
\def\MAT_SUB ##1##2{\xintIrr{\xintSub {##1}{##2}}}%
\def\MAT_MUL {\xintMul }%
\MATid \MAT_invN {#2[I]}%
\MATlet\MAT_invM #2%
\def\MAT_name {#1}%
% \MAT_cntc is the size of the matrix. Will NOT be changed in subroutines.
\MAT_cntc #2[I]\relax
\MAT_cnta \xint_c_i
\MATinvordet_I
}%
\def\MATfloatinvordet #1#2{%
\def\MAT_ZERO {0.e0}%
\def\MAT_DIV {\XINTinFloatDiv }%
\def\MAT_SUB {\XINTinFloatSub }%
\def\MAT_MUL {\XINTinFloatMul }%
\MATfloatid \MAT_invN {#2[I]}%
\MATlet\MAT_invM #2%
\def\MAT_name {#1}%
% \MAT_cntc is the size of the matrix. Will NOT be changed in subroutines.
\MAT_cntc #2[I]\relax
\MAT_cnta \xint_c_i
\MATinvordet_I
}%
\def\MATinvordet_I {\ifnum\MAT_cnta>\MAT_cntc
\expandafter\MATinvordet_II
\else\expandafter\MATinvordet_Ia
\fi }%
\def\MATinv_II {\ifnum\MAT_cnta=\xint_c_i
\expandafter\MATinv_end
\else\expandafter\MATinv_IIa
\fi }%
\def\MATinv_end {\expandafter\MATlet\MAT_name\MAT_invN }%
\def\MATdet_end {\expandafter\let\MAT_name\MAT_det }
\catcode`! 11
\def\MATinv_Ia {%
\MAT_cntb \MAT_cnta\relax
\xintloop
\xintifZero {\MAT_invM [\MAT_cntb,\MAT_cnta]}
{\advance\MAT_cntb \xint_c_i
\ifnum\MAT_cntb>\MAT_cntc \MATinv_!\MATinvordet_I\fi
\iftrue}
{\iffalse}%
\repeat
\MATinv_Ipivot
\ifnum\MAT_cntb>\MAT_cnta \MATinv_exc\fi
\advance\MAT_cnta \xint_c_i
\MATinvordet_I
}%
\def\MATdet_Ia {%
\MAT_cntb \MAT_cnta\relax
\xintloop
\xintifZero {\MAT_invM [\MAT_cntb,\MAT_cnta]}
{\advance\MAT_cntb \xint_c_i
\ifnum\MAT_cntb>\MAT_cntc \MATdet_!\MATinvordet_I\fi
\iftrue}
{\iffalse}%
\repeat
\MATinv_Ipivot
\ifodd\numexpr\MAT_cntb-\MAT_cnta\relax
\edef\MAT_det{\xintOpp {\MAT_det}}%
\fi
\edef\MAT_det {\MAT_MUL {\MAT_pivot}{\MAT_det}}%
\ifnum\MAT_cntb>\MAT_cnta \MATinv_exc\fi
\advance\MAT_cnta \xint_c_i
\MATinvordet_I
}%
\def\MATinv_! #1\fi{\fi
\xintbreakloopanddo
{NOT INVERTIBLE \on@line\typeout{NOT INVERTIBLE \on@line}%
\MATinv_end \def\MAT_tmpa ##1#1{}\MAT_tmpa }%
}%
\def\MATdet_! #1\fi{\fi
\xintbreakloopanddo
{\edef\MAT_det{\MAT_ZERO}%
\MATdet_end \def\MAT_tmpa ##1#1{}\MAT_tmpa }%
}%
\catcode`! 12
\def\MATinv_IIa {%
\advance\MAT_cnta -\xint_c_i
\MATinv_IIpivot
\MATinv_II
}%
\def\MATinv_exc {%
% we optimize as we only need to do in M the indices > \MAT_cnta
% and in N the indices at most \MAT_cntb
\toks0 \expandafter{\the\MAT_cnta}%
\toks2 \expandafter{\the\MAT_cntb}%
% first we do in matrix M, column indices > "a"
\MAT_cntd \MAT_cnta
\xintloop
\ifnum \MAT_cntd<\MAT_cntc
\advance \MAT_cntd \xint_c_i
\toks4 \expandafter{\the\MAT_cntd}%
\expandafter\def\expandafter\MAT_tmpd\expandafter
{\csname MAT@\string\MAT_invM{\the\toks0}{\the\toks4}\endcsname }%
\expandafter\def\expandafter\MAT_tmpe\expandafter
{\csname MAT@\string\MAT_invM{\the\toks2}{\the\toks4}\endcsname }%
\expandafter\let\expandafter\MAT_tmpc\MAT_tmpd
\expandafter\expandafter\expandafter\let\expandafter\MAT_tmpd\MAT_tmpe
\expandafter\let\MAT_tmpe\MAT_tmpc
\repeat
% Then we do in matrix N, column indices <= "b"
\MAT_cntd \xint_c_i % 1
\xintloop
\toks4 \expandafter{\the\MAT_cntd}%
\expandafter\def\expandafter\MAT_tmpd\expandafter
{\csname MAT@\string\MAT_invN{\the\toks0}{\the\toks4}\endcsname }%
\expandafter\def\expandafter\MAT_tmpe\expandafter
{\csname MAT@\string\MAT_invN{\the\toks2}{\the\toks4}\endcsname }%
\expandafter\let\expandafter\MAT_tmpc\MAT_tmpd
\expandafter\expandafter\expandafter\let\expandafter\MAT_tmpd\MAT_tmpe
\expandafter\let\MAT_tmpe\MAT_tmpc
\ifnum \MAT_cntd<\MAT_cntb
\advance\MAT_cntd \xint_c_i
\repeat
}%
\def\MATinv_Ipivot {%
% does pivot simplification on both matrices M and N
% pivot is from matrice M at location (cntb,cnta)
\expandafter\def\expandafter\MAT_tmpa\expandafter {\the\MAT_cnta}%
\expandafter\def\expandafter\MAT_tmpb\expandafter {\the\MAT_cntb}%
\expandafter\let\expandafter\MAT_pivot
\csname MAT@\string\MAT_invM{\the\MAT_cntb}{\the\MAT_cnta}\endcsname
\MAT_cntd \MAT_cnta
\xintloop
\ifnum\MAT_cntd<\MAT_cntc
\advance\MAT_cntd\xint_c_i
% divide in M all entries to the right of the pivot by pivot
\expandafter\def\expandafter\MAT_tmpd\expandafter {\the\MAT_cntd}%
\expandafter
\edef\csname MAT@\string\MAT_invM{\MAT_tmpb}{\MAT_tmpd}\endcsname
{\MAT_DIV{\csname MAT@\string\MAT_invM{\MAT_tmpb}{\MAT_tmpd}\endcsname }
{\MAT_pivot}}%
\repeat
\MAT_cntd \xint_c_i
\xintloop
% divide in N all elements on the "b" row with column indices at most
% equal to "b" by the pivot value
\edef\MAT_tmpd {\the\MAT_cntd}
\expandafter
\edef\csname MAT@\string\MAT_invN{\MAT_tmpb}{\MAT_tmpd}\endcsname
{\MAT_DIV{\csname MAT@\string\MAT_invN{\MAT_tmpb}{\MAT_tmpd}\endcsname }%
{\MAT_pivot}%
}%
\ifnum\MAT_cntd<\MAT_cntb
\advance\MAT_cntd \xint_c_i
\repeat
% we now will simplify the next rows, in both matrices M and N
% Again we don't have to do all entries: >a in M and <= b in N
\MAT_cntd \MAT_cntb % will be increased by 1, row index
\xintloop
{% will not create a group!
\ifnum\MAT_cntd<\MAT_cntc
\advance\MAT_cntd \xint_c_i % we start with the "b+1" row
% We are working with row \cntd
\edef\MAT_tmpd {\the\MAT_cntd}%
% we need the (\cntd, \cnta) entry
\edef\MAT_tmpf
{\csname MAT@\string\MAT_invM{\MAT_tmpd}{\MAT_tmpa}\endcsname }%
% We now multiply by tmpf the cntb row and subtract it from the cntd row
% this sets to zero the (cntd,cnta) entry:
% in matrix M, only need to look at columns to the right
\MAT_cnte\MAT_cnta % necessarily cnta< size of M, as cnta<= cntb<cntd
\advance\MAT_cnte \xint_c_i
\xintloop
\edef\MAT_tmpe {\the\MAT_cnte}%
\expandafter
\edef\csname MAT@\string\MAT_invM{\MAT_tmpd}{\MAT_tmpe}\endcsname
{\MAT_SUB{\csname MAT@\string\MAT_invM{\MAT_tmpd}{\MAT_tmpe}\endcsname }
{\MAT_MUL \MAT_tmpf
{\csname MAT@\string\MAT_invM{\MAT_tmpb}{\MAT_tmpe}\endcsname }}%
}%
\ifnum\MAT_cnte<\MAT_cntc
\advance\MAT_cnte \xint_c_i
\repeat% end of subloop for matrix M, row "d", columns "e>=a"
% we now do the row "d" in matrix N, columns "e<=b"
\MAT_cnte \xint_c_i
\xintloop
\edef\MAT_tmpe {\the\MAT_cnte}%
\expandafter
\edef\csname MAT@\string\MAT_invN{\MAT_tmpd}{\MAT_tmpe}\endcsname
{\MAT_SUB{\csname MAT@\string\MAT_invN{\MAT_tmpd}{\MAT_tmpe}\endcsname }
{\MAT_MUL \MAT_tmpf
{\csname MAT@\string\MAT_invN{\MAT_tmpb}{\MAT_tmpe}\endcsname }}%
}%
\ifnum\MAT_cnte<\MAT_cntb
\advance\MAT_cnte \xint_c_i
\repeat% end of subloop for matrix N, row "d"
}\repeat
}%
\def\MATinv_IIpivot {%
% does pivot simplification on matrices M and N
% M is now upper triangular with 1's on the diagonal
% pivot = 1 is in the \MAT_cnta row. We simplify rows above.
% There is no need to keep track of the computations for M itself
% Only need to read M data and modify rows of N accordingly
\expandafter\def\expandafter\MAT_tmpa\expandafter {\the\MAT_cnta}%
\MAT_cntb \MAT_cnta
\xintloop
{% will not create a group!
\ifnum\MAT_cntb>\xint_c_i
\advance\MAT_cntb -\xint_c_i
\expandafter\def\expandafter\MAT_tmpb\expandafter {\the\MAT_cntb}%
\expandafter\let\expandafter\MAT_tmpf
\csname MAT@\string\MAT_invM{\MAT_tmpb}{\MAT_tmpa}\endcsname
\MAT_cntd\xint_c_i
\xintloop
\expandafter\def\expandafter\MAT_tmpd\expandafter {\the\MAT_cntd}%
\expandafter
\edef\csname MAT@\string\MAT_invN{\MAT_tmpb}{\MAT_tmpd}\endcsname
{\MAT_SUB
{\csname MAT@\string\MAT_invN{\MAT_tmpb}{\MAT_tmpd}\endcsname }
{\MAT_MUL \MAT_tmpf
{\csname MAT@\string\MAT_invN{\MAT_tmpa}{\MAT_tmpd}\endcsname }}%
}%
\ifnum\MAT_cntd<\MAT_cntc
\advance\MAT_cntd \xint_c_i
\repeat
}\repeat
}%
% DISPLAYING MACROS
\makeatother
\def\MATraw {\MATrawwith {\MATrawone}}%
\def\MATrawone {\xintPRaw}%
\def\MATrawwith #1#2{%
\xintListWithSep {; }%
{\xintApply { \MAT_raw_row {#1}#2}{\xintSeq {1}{#2[I]}}}%
}%
\def\MAT_raw_row #1#2#3{%
\xintListWithSep {, }%
{\xintApply { \MAT_raw_one {#1}#2{#3}}{\xintSeq {1}{#2[J]}}}%
}%
\def\MAT_raw_one #1#2#3#4{#1{#2[#3,#4]}}%
%% MATH MODE DISPLAYING
\newcommand\MATdisplay [1][1.25]{\MATdisplaywith [#1]{\MATdisplayone}}
\def\MATdisplayone {\xintSignedFrac}
\newcolumntype\MATdisplaycoltype {c}
\newcolumntype\MATdisplaypreamble [1]{@{}*{#1[J]}\MATdisplaycoltype@{}}
\newcommand\MATdisplaywith [3][1.25]
{\left(\def\arraystretch{#1}%
\begin{array}{\MATdisplaypreamble {#3}}
\xintListWithSep {\\}
{\xintApply { \MAT_display_row {#2}#3}{\xintSeq {1}{#3[I]}}}
\end{array}\right)%
}%
\def\MAT_display_row #1#2#3{%
\xintListWithSep {&}
{\xintApply{ \MAT_display_one {#1}#2{#3}}{\xintSeq {1}{#2[J]}}}%
}%
\def\MAT_display_one #1#2#3#4{#1{#2[#3,#4]}}%
\def\MATminus {\expandafter\MAT_minus_a\romannumeral-`0}%
\def\MAT_minus_a {\futurelet\XINT_token\MAT_minus_b }%
\def\MAT_minus_b {\ifx\XINT_token-\else\phantom{-}\fi }%
\usepackage {siunitx}
\usepackage {numprint}
\newcommand{\MATfloatdisplay}[1][\XINTdigits]
{\MATfloatdisplaywith [#1]{\MATfloatone}}%
\def\MATfloatone #1{\expandafter\MAT_flone\romannumeral-`0#1!}%
\def\MAT_flone #1.#2e#3!{%
\xintSgnFork{\xintiiSgn{\XINT_Abs #3}}%
{}{#1.#2}{#1.#2\times 10^{#3}}}%
\newcommand{\MATfloatdisplaywith}[3][\XINTdigits]
{\left(\edef\MAT_tmpa{#1}%
\begin{array}{\MATdisplaypreamble{#3}}
\xintListWithSep {\\}
{\xintApply { \MAT_fldisplay_row {#2}#3}{\xintSeq {1}{#3[I]}}}%
\end{array}\right)}%
\def\MAT_fldisplay_row #1#2#3{%
\xintListWithSep {&}
{\xintApply{ \MAT_fldisplay_one {#1}#2{#3}}{\xintSeq {1}{#2[J]}}}}%
\def\MAT_fldisplay_one #1#2#3#4{#1{\xintFloat [\MAT_tmpa]{#2[#3,#4]}}}%
\catcode`_ 8
\begin{document}
see the images.
\end{document}
