Perform matrix operations (addition, product, transpose, etc.) in LaTeX?

Question

I know the calc package can perform infix-notation arithmetic in LaTeX... but I want more!

I'd like to perform (not necessarily infix-notation) linear-algebra operations such as scalar multiplication, matrix addition and product in LaTeX, and then print the result in an array or in some matrix environment from amsmath.

Why would I want to do that in LaTeX directly? Why do I not simply use some linear-algebra software, such as Matlab, Mathematica, etc.?

Well, suppose I want to walk my readers through a detailed linear-algebra calculation with many numerical examples. Of course, I could perform all the steps manually first and then hardcode the result of each step in my input file. However, this approach

is prone to errors (LaTeX' arrays are not very user-friendly to typeset),
lacks maintainability (I may decide to change the data in my example, which means I have to modify everything that follows).

Hence my question: Is there a way of easily performing linear-algebra operations in LaTeX?

Ideally, I would like to

mimick Matlab's syntax for defining matrices (using commas as column separators and semicolons as end-of-row characters), performing operations on them, extracting sub matrices, etc.. The syntax could be something like the following:
```
\let\A{1,2,3;4,5,6}
\let\b{1;0;0}
\let\c\matrixprod{\A,\b}
\let\d\submatrix{\c}{(2,1)}
```
have a command that typesets a "matrix object" in an array or matrix environment, e.g.
```
\typesetmatrix[bmatrix]{A}
```
be able to perform operations on matrices of arbitrary--albeit relatively small--dimensions (edit: not just 2x2 and 3x3 as in the calculator package).

Is that currently possible with some package? If not, I'm considering rolling up my sleeves and implementing something, but this would probably prove quite difficult, and I would like to avoid reinventing the wheel :)

Edit about other operations that would be useful:

diag (extraction of the diagonal of a square matrix)
trace
determinant
norm(s)
condition number(s)
inverse

Even more advanced matrix operations that would be awesome, but probably tough to implement:

eigenvalues & eigenvectors,
QR, least squares etc.
SVD,
other common matrix factorisations.

Give me a week or so. What operations do you need to be able to do? I'm planning to start with additions/subtractions and multiplications only, at first, but presumably it shouldn't be too hard to write code for exponentiation and exp (using the power series). What other operations are interesting? — Bruno Le Floch, May 05 '13 at 17:50
@BrunoLeFloch Wow, you're picking up the gauntlet! Have you been planning on coding such a package for a long time? I'd be interested in having a chat about the approach you're thinking about... Ok, other operations that would be useful, aside from exponentiation and exp: trace and determinant. — jub0bs, May 05 '13 at 18:40
For now I will not implement a "proper" package, because it may be better to add a 'matrix' data type to LaTeX3's l3fp. Currently, this package (l3fp) evaluates floating point expressions (expandably). If I can find a good way of defining data types for l3fp, and if I can convince the other members of the LaTeX3 team, it will be possible to write a package that adds matrix operations, letting us make use of the already existing parser: for instance, \fpa_set:Nn \A { matrix(1,2;4,5) } \fpa_set:Nn \B { 1 + \A + \A * \A / 2 - exp(\A) }. — Bruno Le Floch, May 05 '13 at 18:52
I will of course try to make sure that what I implement short-term can still be supported with the 'better' approach, but I cannot give any guarantee. — Bruno Le Floch, May 05 '13 at 18:54
Native support for some linear-algebra stuff in LaTeX3 would be great! Is progress on l3pf accessible to mere mortals somewhere (GitHub perhaps)? — jub0bs, May 05 '13 at 19:00
The package is most likely available in your distribution, in the l3kernel bundle. It's loaded in documents through the expl3 package ('Experimental LaTeX3'), minimal example: \documentclass{article}\usepackage{expl3}\begin{document}\ExplSyntaxOn\fp_eval:n{sin(123/2)}\ExplSyntaxOff\end{document}. Cutting edge progress can be found in the GitHub repository, but the version on CTAN is reasonably up-to-date, as we have been less active in the last month. — Bruno Le Floch, May 05 '13 at 19:21
Why not use LuaLaTeX, or a preprocessor in some other language (e.g. Python/NumPy)? Of course it is possible to do anything in plain LaTeX, but... — leftaroundabout, May 05 '13 at 19:24
@leftaroundabout LuaLaTeX is an option, and I am hoping to eventually add an option to speed up floating point computations by passing them to Lua, but XeTeX is still used quite a lot, and the LuaTeX engine is still a moving target. — Bruno Le Floch, May 07 '13 at 10:55
If you give LuaLaTeX a try maybe this could be helpful for you: http://tex.stackexchange.com/q/73543/10570. Scott H. did a lot of work to create a LaTeX-Lua(TeX) interface for matrix calculations. — Holle, May 07 '13 at 14:44
What I did was export the matrices from a MATLAB/Octave script in LaTeX format, so each time the script is run, it writes the matrices into .tex files, then \input them. — marczellm, Nov 07 '13 at 19:47
@marczellm If the matrices are defined in a Matlab/Octave script file, that's easy enough, but how do you export computation results in text format, ready to be \input inside a .tex file? — jub0bs, Nov 07 '13 at 19:53
Assume M is an Octave matrix. This expression is the key: strcat("\\begin{bmatrix}\n",strrep(strrep(mat2str(M,4)," ","&"),";","\\\\\n")(2:end-1),"\n\\end{bmatrix}\n"); Then write the result to a file. MATLAB and number format needs some tweaking. — marczellm, Nov 07 '13 at 21:53
@marczellm Thanks! I hadn't even thought of doing that. I'll experiment with it. I'm sure it would be possible to do something in Mathematica also. — jub0bs, Nov 07 '13 at 22:28

score 53 · Accepted Answer · edited Jun 11 '17 at 17:53

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.

enter image description here

This is really nice! Are you going to (1) develop this further and (2) make it into a package and put it on CTAN? :) — Svend Tveskæg, May 06 '13 at 01:08
@SvendTveskæg: Probably not before next year. As I mentioned in comments to the question, this is not quite the correct approach for writing a linear algebra package. I would like to reuse the (expandable) expression parser that I wrote for l3fp: one could then write expressions involving matrices, e.g., \matset\B{(1 + \A/5)^5}, and have them be evaluated on the fly. Adding data types to l3fp is not easy to do without slowing down the core stuff (floating points usd internally for rotations etc.). Also, I have to decide how extensible to make the data-type system. — Bruno Le Floch, May 06 '13 at 13:10
I completey missed your comments to the other answer when I wrote the other comment; sorry. — Svend Tveskæg, May 07 '13 at 03:10
Bruno, I must be missing something (I've never used any LaTeX3 before), but I seem unable to compile the code you posted. — jub0bs, May 07 '13 at 08:40
@Jubobs What error message do you get? My guess is that your version of expl3 (hence of l3fp) is too old. — Bruno Le Floch, May 07 '13 at 10:52
Sorry. No errors after updating latex3 packages and siunitx. Thanks! — jub0bs, May 07 '13 at 11:15
Bruno, that's powerful stuff! I'm seriously considering awarding you a bounty, not only to thank you, but to attract more users to your answer. — jub0bs, May 07 '13 at 11:32
I'm not sure but maybe you are "reinvent" a style file which already exists: http://tex.stackexchange.com/q/73543/10570. Perhaps you can use/merge some functions from Scott's style file. — Holle, May 07 '13 at 15:02
@Holle The main difference is that Scott H.'s style file is only meant to work with LuaTeX, as it defers calculations to Lua, while the code I wrote above (and am planning to rewrite eventually) works will all three engines (pdfTeX, XeTeX, LuaTeX). I am planning to provide an option to defer calculations to the LuaTeX engine, but this will have to be requested explicitly by the user, so that results do not depend on where the file is compiled. — Bruno Le Floch, May 07 '13 at 19:07
@Jubobs Thanks! I think that it would make more sense for the bounty to wait until I change the internals of l3fp a bit to accommodate types, and reimplement the package using that mechanism. For that, give me at least a month, perhaps four. — Bruno Le Floch, May 07 '13 at 19:09
@BrunoLeFloch You deserve it :) I have yet to get familiar with LaTeX3 syntax, but I'm somewhat familiar with IEEE-754. Is there any way I can be of help? — jub0bs, May 07 '13 at 19:17
@Jubobs: The actual code of l3fp is hell to read, because of optimisations. Providing a good set of tests would be very helpful (the result of basic arithmetic is defined exactly by the IEEE standard). Otherwise, there are some missing functions, mostly because I couldn't (or didn't have time to) think of a good approach: exact mod(a,b), exactly rounded sqrt(x), atan2, asin, acos, log(x,base), cosh, sinh, tanh. See 'Disclaimer and roadmap' section in l3fp chapter of source3.pdf. — Bruno Le Floch, May 07 '13 at 19:33
@BrunoLeFloch By "16 digits of precision", I assume you're referring to doubles. Do you intend to make l3fp support all binary floating-point formats defined in IEEE-754-2008? — jub0bs, May 31 '13 at 22:53
@Jubobs: IEEE-754-2008 defines both binary and decimal formats. Because TeX's integer arithmetic (despite using binary internally) outputs decimal digits, it is technically much easier, and most likely quite a bit faster, to manipulate decimal floating point numbers rather than binary. I decided to only support 16 digit floating point numbers, which correspond to IEEE's decimal64 format (the min and max exponents are a bit different, I'll fix that). I might support other format for interchange, not calculations. I might add hooks where someone could implement an arbitrary precision unit. — Bruno Le Floch, Jun 03 '13 at 11:08
@BrunoLeFloch: How would I make a macro to return a new matrix that can be used by \matset? — geometrian, Nov 04 '13 at 23:36
@IanMallett Where does the data originate? \matset expects rows separated by semicolons, and cells in a row are separated by commas (or so it seems from a cursory glance at the code). — Bruno Le Floch, Nov 05 '13 at 00:39
@BrunoLeFloch It's a macro that generates a 3x3 rotation matrix. I'm still somewhat new to LaTeX, and I tried permutations on what seemed obvious, but nothing really worked. — geometrian, Nov 05 '13 at 02:11
@IanMallett It is hard to know what to answer without having more precise code. You ask a question which includes the code you have so far, links to the current question, and asks how to get \matset\A{11,12,13;21,22,23;31,32,33} from your code. — Bruno Le Floch, Nov 06 '13 at 15:43
@BrunoLeFloch http://tex.stackexchange.com/questions/142626/applying-matrix-operation-code — geometrian, Nov 06 '13 at 17:19
@BrunoLeFloch can you add an example of getting a value from a matrix? — geometrian, Nov 07 '13 at 00:41
@IanMallett Use \fpm_get:Nnn within your code or \matget in a document setting (I have just added the latter in my last edit). — Bruno Le Floch, Nov 07 '13 at 17:49
@BrunoLeFloch: Awesome thanks; I had tried to do this myself but failed. Working through some related issues (http://tex.stackexchange.com/questions/142875/tikz-dimension-too-large-matrix-transform) now, but thanks for the matrix code! — geometrian, Nov 08 '13 at 00:51

kiss my armpit · Answer 2 · 2013-05-05T16:39:38.867

25

calculator package might help.

enter image description here

edited May 05 '13 at 16:39

answered May 05 '13 at 15:04

kiss my armpit

36,086

@MatthewLeingang: See the last update at the bottom. – kiss my armpit May 05 '13 at 16:39

score 20 · Answer 3 · 2017-03-23T17:07:24.500

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}

I have edited answer in March 2017 but had to remove even more comments to obey max size cap. The edit was needed because since 11/2014 a (non public) macro used here was not defined by xintfrac anymore, and also xinttools needed explicit loading. — , Mar 23 '17 at 17:02
I founded and use now, your brilliant macro for operations on matrices and I would tell you if you have transformed it in a package, which will be interesting because I think that it's very performant. On another hand, I want to tell you if there is some similar macros for computing eigenvalues and eigenvectors — Faouzi Bellalouna, Jan 16 '19 at 04:06
Other thing. One can work wit usual numbers, and with fractions in matrices. How can I hold with more "literal" expressions generic elements of matrices? something like \sqrt{2}, for example ? — Faouzi Bellalouna, Jan 16 '19 at 04:11
@FaouziBellalouna Thanks for your suggestions. Unfortunately, my replies will be a bit disappointing: 1. I have not made it into a package, 2. I have not written additional macros for eigenvalues and eigenvectors, 3. The macros of this answer are only for numeric values. However, it is possible to combine this answer with my package polexpr to do some work with polynomial entries; but not rational functions. Then one can find numerically all real roots. But there is much extra work to do here... — , Jan 16 '19 at 09:03
@FaouziBellalouna You will perhaps be interested by https://tex.stackexchange.com/a/360116/4686 — , Jan 16 '19 at 09:06
Thank you for your kindly reply and beautiful work. I'll see your package and suggestions and tell you my opinion and/or use if any — Faouzi Bellalouna, Jan 16 '19 at 09:55
Your macro \MATrowreduce, based on the Gauss algorithm, is really very interesting, because it gives the rank of the matrices in the general case, and the eigenvalues too for the square matrix case. Now one macro for computing the eigenvectors will be wellcome — Faouzi Bellalouna, Jan 16 '19 at 15:33
@FaouziBellalouna I don't think LU reduction can help immediately in finding eigenvalues. — , Jan 16 '19 at 22:58
Beware that \MATdet produces the same determinant for both matrices 0,4,2;0,8,7;3,1,9 and 3,1,9;0,8,7;0,4,2 which is incorrect. It seems the issue comes from the lines \ifodd\numexpr\MAT_cntb-\MAT_cnta\relax \edef\MAT_det{\xintOpp {\MAT_det}}% \fi I do not know how to fix it yet. — Chilote, Jun 27 '23 at 09:37

score 18 · Answer 4 · edited Jun 11 '17 at 17:52

enter image description here

This is one more option you can check out. Asymptote supports matrix operations, and here is a brief example to demonstrate what is possible. It includes matrix expressions, transpose and inverse. Usage:

define matrices inside the asy environment along with operations on them;
define TeX names with matrixdata function, e.g.: matrixdata("D^T",transpose(a*(b-a)));, here a TeX name for typesetting a matrix is D^T, and the matrix is a result of the matrix expression transpose(a*(b-a)), where a,b were previously defined.
access matrix data inside a standard matrix environment with \mxData{} , e.g. \mxData{D^T}

Example matr.tex:

\documentclass[10pt,a4paper]{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage[inline]{asymptote}
\usepackage[left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}

\begin{asydef}
typedef real[][] matrix;

string smatrixdata(string texName, matrix a){
  string s="\expandafter\gdef\csname texMatrix["+texName+"]\endcsname{";
  for(int i=0;i<a.length;++i){
    for(int j=0;j<a[0].length;++j){
      s+=((j==0)?"":"&")+format("%#.3f",a[i][j]);
    }
    if(i<a.length-1){
      s+="\\"+'\n';
    }
  }
  s+="}";
  return s;
}

void matrixdata(string texName, matrix a){
  tex(smatrixdata(texName, a));
};

\end{asydef}

\gdef\mxData#1{\ifcsname texMatrix[#1]\endcsname\relax\csname texMatrix[#1]\endcsname\relax\fi}


\begin{document}

\begin{asy}
  matrix a={
  {-8,-1,6,0},
  {10,-4,-5,-5},
  {-2,-5,-8,2},
  {-4,7,9,-3},
  };

  matrix b={
  {3,-3,9,-9},
  {-5,9,6,7},
  {-9,-8,-6,1},
  {7,-4,-9,9},
  }; 

  matrix a_squared;
  a_squared=a*a;

  matrixdata("A",a);
  matrixdata("B",b);
  matrixdata("C",a*b);
  matrixdata("A^2",a_squared);
  matrixdata("D",a*(b-a));

  matrixdata("D^T",transpose(a*(b-a)));

  matrix va={ 
    {10},
    {20}
  };

  matrix vb={ 
    {1,2,3,4,5,6}
  };

  matrixdata("va",va);
  matrixdata("vb",vb);
  matrixdata("va*vb",va*vb);
  matrixdata("A^-1",inverse(a));
  matrixdata("A*A^-1",a*inverse(a));

\end{asy}


\begin{align}
A&=\left[
\begin{matrix}
\mxData{A}
\end{matrix}
\right]
\\
B&=\left[
\begin{matrix}
\mxData{B}
\end{matrix}
\right]
\\
C=A\times B&=\left[
\begin{matrix}
\mxData{C}
\end{matrix}
\right]
\\
A^2&=\left[
\begin{matrix}
\mxData{A^2}
\end{matrix}
\right]
\\
D=A\times (B-A)&=\left[
\begin{matrix}
\mxData{D}
\end{matrix}
\right]
\\
D^T&=\left[
\begin{matrix}
\mxData{D^T}
\end{matrix}
\right]
\\
a&=\left[
\begin{matrix}
\mxData{va}
\end{matrix}
\right]
\\
b&=\left[
\begin{matrix}
\mxData{vb}
\end{matrix}
\right]
\\
a\times b&=\left[
\begin{matrix}
\mxData{va*vb}
\end{matrix}
\right]
\\
A^{-1}&=\left[
\begin{matrix}
\mxData{A^-1}
\end{matrix}
\right]
%
\\
A\times A^{-1}&=\left[
\begin{matrix}
\mxData{A*A^-1}
\end{matrix}
\right]
%
\end{align}

\end{document}

To process it with latexmk, create file latexmkrc:

sub asy {return system("asy '$_[0]'");}
add_cus_dep("asy","eps",0,"asy");
add_cus_dep("asy","pdf",0,"asy");
add_cus_dep("asy","tex",0,"asy");

and run latexmk -pdf matr.tex.

As for the other operations you mentioned, feel free to add them as a functions inside the asydef block. I suppose, the C-implementations of the algorithms could be found somewhere, and since the Asymptote syntax is very similar, a translation should not be difficult.

Can Asymptote only deal with numeric values or also with fractions and symbolic math (like sqrt(2))? — Weißer Kater, Feb 07 '20 at 18:02
@user125730: Directly, no, I don't think so. But you can try to combine python asymptote module with Sympy. — g.kov, Feb 08 '20 at 02:44

Bruno · Answer 5 · 2013-05-07T14:25:46.800

16

You could also use the sagetex package, working with the free software Sage.

Pros:

Maintainability
Full power of Sage: matrices, but also polynomials, plots, etc... and any kind of operations (such as the ones required in the edit!)
Don't reinvent the wheel, build a bike!
Easy export to — or integration into — LaTeX
Easy inclusion of the source code if this is needed
Free software!

Cons:

Needs Sage on your computer, or a server to perform computations
Needs some compilation outside LaTeX

edited May 07 '13 at 14:25

answered May 07 '13 at 11:43

Bruno

1,975

Thanks for your answer. I'm not familiar with Sage. I'd like to keep everything in LaTeX if possible, though. – jub0bs May 07 '13 at 23:19
1

As long as you deal with simple computations, it makes sense to stay within LaTeX. The package sagetex can be a good alternative if you need more complicated stuffs. I think that the things you want (such as the determinant or the condition number) may already need some nontrivial implementations, for which I would clearly prefer using a real computer algebra software. Of course, it's up to you! ;-) – Bruno May 08 '13 at 09:24
Yes, some of the useful operations I list in my question are nontrivial, but one may need them for exposition purposes. I'll look into sagetex though; thanks again. – jub0bs May 08 '13 at 10:22

Perform matrix operations (addition, product, transpose, etc.) in LaTeX?

5 Answers5

Matrix operations (muliplications, inverses, determinants) either exactly or with floats (arbitrary precision)

Linked