# Matrix arithmetic

From Rosetta Code
Matrix arithmetic
You are encouraged to solve this task according to the task description, using any language you may know.

For a given matrix, return the determinant and the permanent of the matrix.

The determinant is given by

$\det(A) = \sum_\sigma\sgn(\sigma)\prod_{i=1}^n M_{i,\sigma_i}$

while the permanent is given by

$\operatorname{perm}(A)=\sum_\sigma\prod_{i=1}^n M_{i,\sigma_i}$

In both cases the sum is over the permutations σ of the permutations of 1, 2, ..., n. (A permutation's sign is 1 if there are an even number of inversions and -1 otherwise; see parity of a permutation.)

More efficient algorithms for the determinant are known: LU decomposition, see for example wp:LU decomposition#Computing the determinant. Efficient methods for calculating the permanent are not known.

Cf.

## Contents

C99 code. By no means efficient or reliable. If you need it for serious work, go find a serious library.

`#include <stdio.h>#include <stdlib.h>#include <string.h> double det_in(double **in, int n, int perm){	if (n == 1) return in[0][0]; 	double sum = 0, *m[--n];	for (int i = 0; i < n; i++)		m[i] = in[i + 1] + 1; 	for (int i = 0, sgn = 1; i <= n; i++) {		sum += sgn * (in[i][0] * det_in(m, n, perm));		if (i == n) break; 		m[i] = in[i] + 1;		if (!perm) sgn = -sgn;	}	return sum;} /* wrapper function */double det(double *in, int n, int perm){	double *m[n];	for (int i = 0; i < n; i++)		m[i] = in + (n * i); 	return det_in(m, n, perm);} int main(void){	double x[] = {	0, 1, 2, 3, 4,			5, 6, 7, 8, 9,			10, 11, 12, 13, 14,			15, 16, 17, 18, 19,			20, 21, 22, 23, 24 }; 	printf("det:  %14.12g\n", det(x, 5, 0));	printf("perm: %14.12g\n", det(x, 5, 1)); 	return 0;}`

A method to calculate determinant that might actually be usable:

`#include <stdio.h>#include <stdlib.h>#include <tgmath.h> void showmat(const char *s, double **m, int n){	printf("%s:\n", s);	for (int i = 0; i < n; i++) {		for (int j = 0; j < n; j++)			printf("%12.4f", m[i][j]);			putchar('\n');	}} int trianglize(double **m, int n){	int sign = 1;	for (int i = 0; i < n; i++) {		int max = 0; 		for (int row = i; row < n; row++)			if (fabs(m[row][i]) > fabs(m[max][i]))				max = row; 		if (max) {			sign = -sign;			double *tmp = m[i];			m[i] = m[max], m[max] = tmp;		} 		if (!m[i][i]) return 0; 		for (int row = i + 1; row < n; row++) {			double r = m[row][i] / m[i][i];			if (!r)	continue; 			for (int col = i; col < n; col ++)				m[row][col] -= m[i][col] * r;		}	}	return sign;} double det(double *in, int n){	double *m[n];	m[0] = in; 	for (int i = 1; i < n; i++)		m[i] = m[i - 1] + n; 	showmat("Matrix", m, n); 	int sign = trianglize(m, n);	if (!sign)		return 0; 	showmat("Upper triangle", m, n); 	double p = 1;	for (int i = 0; i < n; i++)		p *= m[i][i];	return p * sign;} #define N 18int main(void){	double x[N * N];	srand(0);	for (int i = 0; i < N * N; i++)		x[i] = rand() % N; 	printf("det: %19f\n", det(x, N));	return 0;}`

This requires the modules from the Permutations and Permutations by swapping Tasks.

Translation of: Python
`import std.algorithm, std.range, std.traits, permutations2,       permutations_by_swapping1; auto prod(Range)(/*in*/ Range r) /*pure nothrow*/ {    return reduce!q{a * b}(cast(ForeachType!Range)1, r);} T permanent(T)(in T[][] a) /*pure nothrow*/in {    foreach (const row; a)        assert(row.length == a[0].length);} body {    auto r = iota(cast()a.length);    T tot = 0;    foreach (sigma; permutations(r.array()))        tot += r.map!(i => a[i][sigma[i]])().prod();    return tot;} T determinant(T)(in T[][] a) /*pure nothrow*/in {    foreach (const row; a)        assert(row.length == a[0].length);} body {    immutable n = a.length;    auto r = iota(n);    T tot = 0;    //foreach (sigma, sign; spermutations(n)) {    foreach (sigma_sign; spermutations(n)) {        const sigma = sigma_sign[0];        immutable sign = sigma_sign[1];        tot += sign * r.map!(i => a[i][sigma[i]])().prod();    }    return tot;} void main() {    import std.stdio;     foreach (const a; [[[1, 2],                        [3, 4]],                        [[1, 2, 3, 4],                        [4, 5, 6, 7],                        [7, 8, 9, 10],                        [10, 11, 12, 13]],                        [[ 0,  1,  2,  3,  4],                        [ 5,  6,  7,  8,  9],                        [10, 11, 12, 13, 14],                        [15, 16, 17, 18, 19],                        [20, 21, 22, 23, 24]]]) {        writefln("[%([%(%2s, %)],\n %)]]", a);        writefln("Permanent: %s, determinant: %s\n",                 permanent(a), determinant(a));    }}`
Output:
```[[ 1,  2],
 [ 3,  4]]
Permanent: 10, determinant: -2

[[ 1,  2,  3,  4],
 [ 4,  5,  6,  7],
 [ 7,  8,  9, 10],
 [10, 11, 12, 13]]
Permanent: 29556, determinant: 0

[[ 0,  1,  2,  3,  4],
 [ 5,  6,  7,  8,  9],
 [10, 11, 12, 13, 14],
 [15, 16, 17, 18, 19],
 [20, 21, 22, 23, 24]]
Permanent: 6778800, determinant: 0```

J has a conjunction for defining verbs which can act as determinant. This conjunction is symbolized as a space followed by a dot.

For people who do not care to sort out what "recursive expansion by minors" means: `-/ .* y` gives the determinant of y and `+/ .* y` gives the permanent of y.

For example, given the matrix:

`   i. 5 5 0  1  2  3  4 5  6  7  8  910 11 12 13 1415 16 17 18 1920 21 22 23 24`

Its determinant is 0. When we use IEEE floating point, we only get an approximation of this result:

`   -/ .* i. 5 5_1.30277e_44`

If we use exact (rational) arithmetic, we get a precise result:

`   -/ .* i. 5 5x0`

The permanent does not have this problem in this example (the matrix contains no negative values and permanent does not use subtraction):

`   +/ .* i. 5 56778800`

As an aside, note that for specific verbs (like -/ .*) J uses an algorithm which is more efficient than the brute force approach.

```П4	ИПE	П2	КИП0	ИП0	П1	С/П	ИП4	/	КП2
L1	06	ИПE	П3	ИП0	П1	Сx	КП2	L1	17
ИП0	ИП2	+	П1	П2	ИП3	-	x#0	34	С/П
ПП	80	БП	21	КИП0	ИП4	С/П	КИП2	-	*
П4	ИП0	П3	x#0	35	Вx	С/П	КИП2	-	<->
/	КП1	L3	45	ИП1	ИП0	+	П3	ИПE	П1
П2	КИП1	/-/	ПП	80	ИП3	+	П3	ИП1	-
x=0	61	ИП0	П1	КИП3	КП2	L1	74	БП	12
ИП0	<->	^	КИП3	*	КИП1	+	КП2	->	L0
82	->	П0	В/О
```

This program calculates the determinant of the matrix of order <= 5. Prior to startup, РE entered 13, entered the order of the matrix Р0, and the elements are introduced with the launch of the program after one of them, the last on the screen will be determinant. Permanent is calculated in this way.

Determinant is a built in function Det

`Permanent[m_List] :=    With[{v = Array[x, Length[m]]},      Coefficient[Times @@ (m.v), Times @@ v]  ]`

`a: matrix([2, 9, 4], [7, 5, 3], [6, 1, 8])$ determinant(a);-360 permanent(a);900`

The determinant is built in:

`matdet(M)`

and the permanent can be defined as

`matperm(M)=my(n=#M,t);sum(i=1,n!,t=numtoperm(n,i);prod(j=1,n,M[j,t[j]]))`

Using the module file spermutations.py from Permutations by swapping. The algorithm for the determinant is a more literal translation of the expression in the task description and the Wikipedia reference.

`from itertools import permutationsfrom operator import mulfrom math import fsumfrom spermutations import spermutations def prod(lst):    return reduce(mul, lst, 1) def perm(a):    n = len(a)    r = range(n)    s = permutations(r)    return fsum(prod(a[i][sigma[i]] for i in r) for sigma in s) def det(a):    n = len(a)    r = range(n)    s = spermutations(n)    return fsum(sign * prod(a[i][sigma[i]] for i in r)                for sigma, sign in s) if __name__ == '__main__':    from pprint import pprint as pp     for a in (             [             [1, 2],              [3, 4]],              [             [1, 2, 3, 4],             [4, 5, 6, 7],             [7, 8, 9, 10],             [10, 11, 12, 13]],                     [             [ 0,  1,  2,  3,  4],             [ 5,  6,  7,  8,  9],             [10, 11, 12, 13, 14],             [15, 16, 17, 18, 19],             [20, 21, 22, 23, 24]],        ):        print('')        pp(a)        print('Perm: %s Det: %s' % (perm(a), det(a)))`
Sample output
```[[1, 2], [3, 4]]
Perm: 10 Det: -2

[[1, 2, 3, 4], [4, 5, 6, 7], [7, 8, 9, 10], [10, 11, 12, 13]]
Perm: 29556 Det: 0

[[0, 1, 2, 3, 4],
 [5, 6, 7, 8, 9],
 [10, 11, 12, 13, 14],
 [15, 16, 17, 18, 19],
 [20, 21, 22, 23, 24]]
Perm: 6778800 Det: 0```

The second matrix above is that used in the Tcl example. The third matrix is from the J language example. Note that the determinant seems to be 'exact' using this method of calculation without needing to resort to other than Pythons default numbers.

` #lang racket(require math)(define determinant matrix-determinant) (define (permanent M)  (define n (matrix-num-rows M))  (for/sum ([σ (in-permutations (range n))])    (for/product ([i n] [σi σ])      (matrix-ref M i σi)))) `

The determinant is provided by the linear algebra package in Tcllib. The permanent (being somewhat less common) requires definition, but is easily described:

Library: Tcllib (Package: math::linearalgebra)
Library: Tcllib (Package: struct::list)
`package require math::linearalgebrapackage require struct::list proc permanent {matrix} {    for {set plist {};set i 0} {$i<[llength $matrix]} {incr i} {	lappend plist $i    }    foreach p [::struct::list permutations $plist] {	foreach i $plist j $p {	    lappend prod [lindex $matrix $i $j]	}	lappend sum [::tcl::mathop::* {*}$prod[set prod {}]]    }    return [::tcl::mathop::+ {*}$sum]}`

Demonstrating with a sample matrix:

`set mat {    {1 2 3 4}    {4 5 6 7}    {7 8 9 10}    {10 11 12 13}}puts [::math::linearalgebra::det $mat]puts [permanent $mat]`
Output:
```1.1315223609263888e-29
29556
```

Tweet
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 2, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "github"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "english_confidence": 0.8654286861419678, "perplexity_flag": "middle"}