![在 LaTeX 中执行矩阵运算(加法、乘法、转置等)?](https://linux22.com/image/267830/%E5%9C%A8%20LaTeX%20%E4%B8%AD%E6%89%A7%E8%A1%8C%E7%9F%A9%E9%98%B5%E8%BF%90%E7%AE%97%EF%BC%88%E5%8A%A0%E6%B3%95%E3%80%81%E4%B9%98%E6%B3%95%E3%80%81%E8%BD%AC%E7%BD%AE%E7%AD%89%EF%BC%89%EF%BC%9F.png)
我知道该calc
包可以在 LaTeX 中执行中缀表示法算术......但我想要更多!
我想在 LaTeX 中执行(不一定是用中缀表示法)线性代数运算,例如标量乘法、矩阵加法和乘积,然后在array
或 的某个矩阵环境中打印结果amsmath
。
为什么我要直接在 LaTeX 中执行此操作?为什么我不简单地使用一些线性代数软件,例如 Matlab、Mathematica 等?
好吧,假设我想通过许多数值示例向读者介绍一个详细的线性代数计算。当然,我可以先手动执行所有步骤,然后在输入文件中硬编码每个步骤的结果。但是,这种方法
- 容易出错(LaTeX 的
array
排版不太方便用户使用), - 缺乏可维护性(我可能决定更改示例中的数据,这意味着我必须修改后续的所有内容)。
因此我的问题是:有没有办法在 LaTeX 中轻松地执行线性代数运算?
理想情况下,我希望
模仿 Matlab 定义矩阵的语法(使用逗号作为列分隔符,使用分号作为行尾字符),对其执行操作,提取子矩阵等。语法可能类似于以下内容:
\let\A{1,2,3;4,5,6} \let\b{1;0;0} \let\c\matrixprod{\A,\b} \let\d\submatrix{\c}{(2,1)}
array
有一个在或环境中排版“矩阵对象”的命令matrix
,例如\typesetmatrix[bmatrix]{A}
能够对任意维度(尽管相对较小)的矩阵执行运算(编辑:不仅仅是包装中的 2x2 和 3x3
calculator
)。
目前是否可以通过某些软件包实现这一点?如果不行,我正在考虑撸起袖子实施一些措施,但这可能相当困难,而且我想避免重新发明轮子 :)
编辑关于其他有用的操作:
- diag(提取方阵的对角线)
- 痕迹
- 决定因素
- 规范
- 条件编号
- 逆
甚至更高级的矩阵运算也很棒,但可能很难实现:
- 特征值和特征向量,
- QR、最小二乘法等。
- 奇异值分解,
- 其他常见的矩阵分解。
答案1
以下是一些用于操作任意大小矩阵的代码。目前,它可以执行加法、减法和乘法(以及获取单个条目和转置矩阵等)。条目是l3fp
支持的浮点数(精度为 16 位)。
% 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}
编辑:添加\matget
,提取矩阵中的单个条目。
答案2
calculator
包可能会有帮助。
答案3
矩阵运算(乘法、逆运算、行列式)可以精确计算或使用浮点数(任意精度)
这个 2013 年 11 月的答案于 2017 年 3 月进行了编辑,因为早在 2014 年底,从 xint 中删除宏就导致代码无法使用,并且自 1.1(2014/10/28)以来 xintfrac 不会自动加载 xinttools。
\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}