Determinant and permanent
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
while the permanent is given by
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.
C
C99 code. By no means efficient or reliable. If you need it for serious work, go find a serious library. <lang C>#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; }</lang> A method to calculate determinant that might actually be usable: <lang c>#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 18
int 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; }</lang>
D
This requires the modules from the Permutations and Permutations by swapping Tasks.
<lang d>import std.algorithm, std.range, std.traits, permutations2,
permutations_by_swapping1;
auto prod(Range)(Range r) nothrow @safe @nogc {
return reduce!q{a * b}(ForeachType!Range(1), r);
}
T permanent(T)(in T[][] a) nothrow @safe in {
assert(a.all!(row => row.length == a[0].length));
} body {
auto r = a.length.iota;
T tot = 0;
foreach (const sigma; r.array.permutations)
tot += r.map!(i => a[i][sigma[i]]).prod;
return tot;
}
T determinant(T)(in T[][] a) nothrow in {
assert(a.all!(row => row.length == a[0].length));
} body {
immutable n = a.length;
auto r = n.iota;
T tot = 0;
//foreach (sigma, sign; n.spermutations) {
foreach (const sigma_sign; n.spermutations) {
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 (/*immutable*/ 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",
a.permanent, a.determinant);
}
}</lang>
- 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
Fortran
Please find the compilation and example run at the start of the comments in the f90 source. Thank you.
<lang FORTRAN> !-*- mode: compilation; default-directory: "/tmp/" -*- !Compilation started at Sat May 18 23:25:42 ! !a=./F && make $a && $a < unixdict.txt !f95 -Wall -ffree-form F.F -o F ! j example, determinant: 7.00000000 ! j example, permanent: 5.00000000 ! maxima, determinant: -360.000000 ! maxima, permanent: 900.000000 ! !Compilation finished at Sat May 18 23:25:43
! NB. example computed by J ! NB. fixed seed random matrix ! _2+3 3?.@$5 ! 2 _1 1 !_1 _2 1 !_1 _1 _1 ! ! (-/ .*)_2+3 3?.@$5 NB. determinant !7 ! (+/ .*)_2+3 3?.@$5 NB. permanent !5
!maxima example !a: matrix([2, 9, 4], [7, 5, 3], [6, 1, 8])$ !determinant(a); !-360 ! !permanent(a); !900
! compute permanent or determinant
program f
implicit none real, dimension(3,3) :: j, m data j/ 2,-1, 1,-1,-2, 1,-1,-1,-1/ data m/2, 9, 4, 7, 5, 3, 6, 1, 8/ write(6,*) 'j example, determinant: ',det(j,3,-1) write(6,*) 'j example, permanent: ',det(j,3,1) write(6,*) 'maxima, determinant: ',det(m,3,-1) write(6,*) 'maxima, permanent: ',det(m,3,1)
contains
recursive function det(a,n,permanent) result(accumulation)
! setting permanent to 1 computes the permanent.
! setting permanent to -1 computes the determinant.
real, dimension(n,n), intent(in) :: a
integer, intent(in) :: n, permanent
real, dimension(n-1, n-1) :: b
real :: accumulation
integer :: i, sgn
if (n .eq. 1) then
accumulation = a(1,1)
else
accumulation = 0
sgn = 1
do i=1, n
b(:, :(i-1)) = a(2:, :i-1)
b(:, i:) = a(2:, i+1:)
accumulation = accumulation + sgn * a(1, i) * det(b, n-1, permanent)
sgn = sgn * permanent
enddo
endif
end function det
end program f </lang>
FunL
From the task description: <lang funl>def sgn( p ) = product( (if s(0) < s(1) xor i(0) < i(1) then -1 else 1) | (s, i) <- p.combinations(2).zip( (0:p.length()).combinations(2) ) )
def perm( m ) = sum( product(m(i, sigma(i)) | i <- 0:m.length()) | sigma <- (0:m.length()).permutations() )
def det( m ) = sum( sgn(sigma)*product(m(i, sigma(i)) | i <- 0:m.length()) | sigma <- (0:m.length()).permutations() )</lang>
Laplace expansion (recursive): <lang funl>def perm( m )
| m.length() == 1 and m(0).length() == 1 = m(0, 0) | otherwise = sum( m(i, 0)*perm(m(0:i, 1:m.length()) + m(i+1:m.length(), 1:m.length())) | i <- 0:m.length() )
def det( m )
| m.length() == 1 and m(0).length() == 1 = m(0, 0) | otherwise = sum( (-1)^i*m(i, 0)*det(m(0:i, 1:m.length()) + m(i+1:m.length(), 1:m.length())) | i <- 0:m.length() )</lang>
Test using the first set of definitions (from task description): <lang funl>matrices = [
( (1, 2), (3, 4)), ( (-2, 2, -3), (-1, 1, 3), ( 2, 0, -1)), ( ( 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)) ]
for m <- matrices
println( m, 'perm: ' + perm(m), 'det: ' + det(m) )</lang>
- Output:
((1, 2), (3, 4)), perm: 10, det: -2 ((-2, 2, -3), (-1, 1, 3), (2, 0, -1)), perm: 10, det: 18 ((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
Go
Importing the permute packge from the Permutations by swapping task. <lang go>package main
import (
"fmt" "permute"
)
func determinant(m [][]float64) (d float64) {
p := make([]int, len(m))
for i := range p {
p[i] = i
}
it := permute.Iter(p)
for s := it(); s != 0; s = it() {
pr := 1.
for i, σ := range p {
pr *= m[i][σ]
}
d += float64(s) * pr
}
return
}
func permanent(m [][]float64) (d float64) {
p := make([]int, len(m))
for i := range p {
p[i] = i
}
it := permute.Iter(p)
for s := it(); s != 0; s = it() {
pr := 1.
for i, σ := range p {
pr *= m[i][σ]
}
d += pr
}
return
}
var m2 = [][]float64{
{1, 2},
{3, 4}}
var m3 = [][]float64{
{2, 9, 4},
{7, 5, 3},
{6, 1, 8}}
func main() {
fmt.Println(determinant(m2), permanent(m2)) fmt.Println(determinant(m3), permanent(m3))
}</lang>
- Output:
-2 10 -360 900
J
J has a conjunction for defining verbs which can act as determinant. This conjunction is symbolized as a space followed by a dot.
For example, given the matrix:
<lang J> i. 5 5
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</lang>
Its determinant is 0. When we use IEEE floating point, we only get an approximation of this result:
<lang J> -/ .* i. 5 5 _1.30277e_44</lang>
If we use exact (rational) arithmetic, we get a precise result:
<lang J> -/ .* i. 5 5x 0</lang>
The permanent does not have this problem in this example (the matrix contains no negative values and permanent does not use subtraction):
<lang J> +/ .* i. 5 5 6778800</lang>
As an aside, note that for specific verbs (like -/ .*) J uses an algorithm which is more efficient than the brute force approach.
jq
Recursive definitions
<lang jq># Eliminate row i and row j def except(i;j):
reduce del(.[i])[] as $row ([]; . + [$row | del(.[j]) ] );
def det:
def parity(i): if i % 2 == 0 then 1 else -1 end;
if length == 1 and (.[0] | length) == 1 then .[0][0]
else . as $m
| reduce range(0; length) as $i
(0; . + parity($i) * $m[0][$i] * ( $m | except(0;$i) | det) )
end ;
def perm:
if length == 1 and (.[0] | length) == 1 then .[0][0]
else . as $m
| reduce range(0; length) as $i
(0; . + $m[0][$i] * ( $m | except(0;$i) | perm) )
end ;</lang>
Examples <lang jq>def matrices:
[ [1, 2], [3, 4]],
[ [-2, 2, -3], [-1, 1, 3], [ 2, 0, -1]],
[ [ 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]]
"Determinants: ", (matrices | det), "Permanents: ", (matrices | perm)</lang>
- Output:
<lang sh>$ jq -n -r -f Matrix_arithmetic.jq Determinants: -2 18 0 0 Permanents: 10 10 29556 6778800</lang>
Determinant via LU Decomposition
The following uses the jq infrastructure at LU decomposition to achieve an efficient implementation of det/0: <lang jq># Requires lup/0 def det:
def product_diagonal: . as $m | reduce range(0;length) as $i (1; . * $m[$i][$i]); def tidy: if . == -0 then 0 else . end; lup | (.[0]|product_diagonal) as $l | if $l == 0 then 0 else $l * (.[1]|product_diagonal) | tidy end ;
</lang> Examples
Using matrices/0 as defined above: <lang jq>matrices | det</lang>
- Output:
$ /usr/local/bin/jq -M -n -f LU.rc 2 -18 0 0
Julia
The determinant of a matrix A can be computed by the built-in function
<lang julia>det(A)</lang>
The following function computes the permanent of a matrix A from the definition: <lang julia>function perm(A)
m, n = size(A)
if m != n; throw(ArgumentError("permanent is for square matrices only")); end
sum(σ -> prod(i -> A[i,σ[i]], 1:n), permutations(1:n))
end</lang>
Example output: <lang julia>julia> A = [2 9 4; 7 5 3; 6 1 8] julia> det(A), perm(A) (-360.0,900)</lang>
МК-61/52
П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.
Maple
<lang Maple>M:=<<2|9|4>,<7|5|3>,<6|1|8>>:
with(LinearAlgebra):
Determinant(M); Permanent(M);</lang> Output:
-360
900
Mathematica
Determinant is a built in function Det <lang Mathematica>Permanent[m_List] :=
With[{v = Array[x, Length[m]]},
Coefficient[Times @@ (m.v), Times @@ v]
]</lang>
Maxima
<lang maxima>a: matrix([2, 9, 4], [7, 5, 3], [6, 1, 8])$
determinant(a); -360
permanent(a); 900</lang>
Nim
Using the permutationsswap module from Permutations by swapping: <lang nim>import sequtils, permutationsswap
type Matrix[M,N: static[int]] = array[M, array[N, float]]
proc det[M,N](a: Matrix[M,N]): float =
let n = toSeq 0..a.high for sigma, sign in n.permutations: var x = sign.float for i in n: x *= a[i][sigma[i]] result += x
proc perm[M,N](a: Matrix[M,N]): float =
let n = toSeq 0..a.high for sigma, sign in n.permutations: var x = 1.0 for i in n: x *= a[i][sigma[i]] result += x
const
a = [ [1.0, 2.0]
, [3.0, 4.0]
]
b = [ [ 1.0, 2, 3, 4]
, [ 4.0, 5, 6, 7]
, [ 7.0, 8, 9, 10]
, [10.0, 11, 12, 13]
]
c = [ [ 0.0, 1, 2, 3, 4]
, [ 5.0, 6, 7, 8, 9]
, [10.0, 11, 12, 13, 14]
, [15.0, 16, 17, 18, 19]
, [20.0, 21, 22, 23, 24]
]
echo "perm: ", a.perm, " det: ", a.det echo "perm: ", b.perm, " det: ", b.det echo "perm: ", c.perm, " det: ", c.det</lang> Output:
perm: 10.0 det: -2.0 perm: 29556.0 det: 0.0 perm: 6778800.0 det: 0.0
PARI/GP
The determinant is built in: <lang parigp>matdet(M)</lang> and the permanent can be defined as <lang parigp>matperm(M)=my(n=#M,t);sum(i=1,n!,t=numtoperm(n,i);prod(j=1,n,M[j,t[j]]))</lang>
Perl 6
Uses the permutations generator from the Permutations by swapping task. This implementation is naive and brute-force (slow) but exact.
<lang perl6>sub insert( $x, @xs) { [@xs[0..$_-1], $x, @xs[$_..*]] for 0..@xs } sub order ($sg, @xs) { $sg > 0 ?? @xs !! @xs.reverse }
multi σ_permutations ([]) { [] => 1 }
multi σ_permutations ([$x, *@xs]) {
σ_permutations(@xs).map({ order($_.value, insert($x, $_.key)) }) Z=> (1,-1) xx *
}
sub m_arith ( @a, $op ) {
note "Not a square matrix" and return
if [||] map { @a.elems cmp @a[$_].elems }, ^@a;
[+] map {
my $permutation = .key;
my $term = $op eq 'perm' ?? 1 !! .value;
for $permutation.kv -> $i, $j { $term *= @a[$i][$j] };
$term
}, σ_permutations [^@a];
}
- Testing ###########
my @tests = (
[
[ 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 ]
]
);
sub dump (@matrix) {
say $_».fmt: "%3s" for @matrix, ;
}
for @tests -> @matrix {
say 'Matrix:'; @matrix.&dump; say "Determinant:\t", @matrix.&m_arith: <det>; say "Permanent: \t", @matrix.&m_arith: <perm>; say '-' x 25;
} </lang>
Output
Matrix: 1 2 3 4 Determinant: -2 Permanent: 10 ------------------------- Matrix: 1 2 3 4 4 5 6 7 7 8 9 10 10 11 12 13 Determinant: 0 Permanent: 29556 ------------------------- Matrix: 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 Determinant: 0 Permanent: 6778800 -------------------------
PowerShell
<lang PowerShell> function det-perm ($array) {
if($array) {
$size = $array.Count
function prod($A) {
$prod = 1
if($A) { $A | foreach{$prod *= $_} }
$prod
}
function generate($sign, $n, $A) {
if($n -eq 1) {
$i = 0
$prod = prod @($A | foreach{$array[$i++][$_]})
[pscustomobject]@{det = $sign*$prod; perm = $prod}
}
else{
for($i = 0; $i -lt ($n - 1); $i += 1) {
generate $sign ($n - 1) $A
if($n % 2 -eq 0){
$i1, $i2 = $i, ($n-1)
$A[$i1], $A[$i2] = $A[$i2], $A[$i1]
}
else{
$i1, $i2 = 0, ($n-1)
$A[$i1], $A[$i2] = $A[$i2], $A[$i1]
}
$sign *= -1
}
generate $sign ($n - 1) $A
}
}
$det = $perm = 0
generate 1 $size @(0..($size-1)) | foreach{
$det += $_.det
$perm += $_.perm
}
[pscustomobject]@{det = "$det"; perm = "$perm"}
} else {Write-Error "empty array"}
} det-perm 5 det-perm @(@(1,0,0),@(0,1,0),@(0,0,1)) det-perm @(@(0,0,1),@(0,1,0),@(1,0,0)) det-perm @(@(4,3),@(2,5)) det-perm @(@(2,5),@(4,3)) det-perm @(@(4,4),@(2,2)) </lang> Output:
det perm --- ---- 5 5 1 1 -1 1 14 26 -14 26 0 16
Python
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.
<lang python>from itertools import permutations from operator import mul from math import fsum from 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)))</lang>
- 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.
Racket
<lang racket>
- 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))))
</lang>
REXX
<lang rexx>/* REXX ***************************************************************
- Test the two functions determinant and permanent
- using the matrix specifications shown for other languages
- 21.05.2013 Walter Pachl
- /
Call test ' 1 2',
' 3 4',2
Call test ' 1 2 3 4',
' 4 5 6 7',
' 7 8 9 10',
'10 11 12 13',4
Call test ' 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',5
Exit
test: /**********************************************************************
- Show the given matrix and compute and show determinant and permanent
- /
Parse Arg as,n asc=as Do i=1 To n
ol= Do j=1 To n Parse Var asc a.i.j asc ol=ol right(a.i.j,3) End Say ol End
Say 'determinant='right(determinant(as),7) Say ' permanent='right(permanent(as),7) Say copies('-',50) Return</lang>
<lang rexx>/* REXX ***************************************************************
- determinant.rex
- compute the determinant of the given square matrix
- Input: as: the representation of the matrix as vector (n**2 elements)
- 21.05.2013 Walter Pachl
- /
Parse Arg as
n=sqrt(words(as))
Do i=1 To n
Do j=1 To n
Parse Var as a.i.j as
End
End
Select
When n=2 Then det=a.1.1*a.2.2-a.1.2*a.2.1
When n=3 Then det= a.1.1*a.2.2*a.3.3,
+a.1.2*a.2.3*a.3.1,
+a.1.3*a.2.1*a.3.2,
-a.1.3*a.2.2*a.3.1,
-a.1.2*a.2.1*a.3.3,
-a.1.1*a.2.3*a.3.2
Otherwise Do
det=0
Do k=1 To n
det=det+((-1)**(k+1))*a.1.k*determinant(subm(k))
End
End
End
Return det
subm: Procedure Expose a. n /**********************************************************************
- compute the submatrix resulting when row 1 and column k are removed
- Input: a.*.*, k
- Output: bs the representation of the submatrix as vector
- /
Parse Arg k
bs=
do i=2 To n
Do j=1 To n
If j=k Then Iterate
bs=bs a.i.j
End
End
Return bs
sqrt: Procedure /**********************************************************************
- compute and return the (integer) square root of the given argument
- terminate the program if the argument is not a square
- /
Parse Arg nn Do n=1 By 1 while n*n<nn End If n*n=nn Then Return n Else Do Say 'invalid number of elements:' nn 'is not a square.' Exit End</lang>
<lang rexx>/* REXX ***************************************************************
- permanent.rex
- compute the permanent of a matrix
- I found an algorithm here:
- http://www.codeproject.com/Articles/21282/Compute-Permanent-of-a-Matrix-with-Ryser-s-Algorit
- see there for the original author.
- translated it to REXX (hopefully correctly) to REXX
- and believe that I can "publish" it here, on rosettacode
- when I look at the copyright rules shown there:
- http://www.codeproject.com/info/cpol10.aspx
- 20.05.2013 Walter Pachl
- /
Call init arg(1) /* initialize the matrix (n and a.* */ sum=0 rowsumprod=0 rowsum=0 chi.=0 c=2**n Do k=1 To c-1 /* loop all 2^n submatrices of A */
rowsumprod = 1
chis=dec2binarr(k,n) /* characteristic vector */
Do ci=0 By 1 While chis<>
Parse Var chis chi.ci chis
End
Do m=0 To n-1 /* loop columns of submatrix #k */
rowsum = 0
Do p=0 To n-1 /* loop rows and compute rowsum */
mnp=m*n+p
rowsum=rowsum+chi.p*A.mnp
End
rowsumprod=rowsumprod*rowsum /* update product of rowsums */
/* (optional -- use for sparse matrices) */
/* if (rowsumprod == 0) break; */
End
sum=sum+((-1)**(n-chi.n))*rowsumprod
End
Return sum /**********************************************************************
- Notes
- 1.The submatrices are chosen by use of a characteristic vector chi
- (only the columns are considered, where chi[p] == 1).
- To retrieve the t from Ryser's formula, we need to save the number
- n-t, as is done in chi[n]. Then we get t = n - chi[n].
- 2.The matrix parameter A is expected to be a one-dimensional integer
- array -- should the matrix be encoded row-wise or column-wise?
- -- It doesn't matter. The permanent is invariant under
- row-switching and column-switching, and it is Screenshot
- - per_inv.gif .
- 3.Further enhancements: If any rowsum equals zero,
- the entire rowsumprod becomes zero, and thus the m-loop can be broken.
- Since if-statements are relatively expensive compared to integer
- operations, this might save time only for sparse matrices
- (where most entries are zeros).
- 4.If anyone finds a polynomial algorithm for permanents,
- he will get rich and famous (at least in the computer science world).
- /
/**********************************************************************
- At first, we need to transform a decimal to a binary array
- with an additional element
- (the last one) saving the number of ones in the array:
- /
dec2binarr: Procedure
Parse Arg n,dim ol='n='n 'dim='dim res.=0 pos=dim-1 Do While n>0 res.pos=n//2 res.dim=res.dim+res.pos n=n%2 pos=pos-1 End res_s= Do i=0 To dim res_s=res_s res.i End Return res_s
init: Procedure Expose a. n /**********************************************************************
- a.* (starting with index 0) contains all array elements
- n is the dimension of the square matrix
- /
Parse Arg as n=sqrt(words(as)) a.=0 Do ai=0 By 1 While as<>
Parse Var as a.ai as End
Return
sqrt: Procedure /**********************************************************************
- compute and return the (integer) square root of the given argument
- terminate the program if the argument is not a square
- /
Parse Arg nn Do n=1 By 1 while n*n<nn End If n*n=nn Then Return n Else Do Say 'invalid number of elements:' nn 'is not a square.' Exit End</lang>
Output:
1 2 3 4 determinant= -2 permanent= 10 -------------------------------------------------- 1 2 3 4 4 5 6 7 7 8 9 10 10 11 12 13 determinant= 0 permanent= 29556 -------------------------------------------------- 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 determinant= 0 permanent=6778800 --------------------------------------------------
Ruby
Matrix in the standard library provides a method for the determinant, but not for the permanent. <lang ruby>require 'matrix'
class Matrix
# Add "permanent" method to Matrix class
def permanent
r = (0...row_count).to_a # [0,1] (first example), [0,1,2,3] (second example)
r.permutation.inject(0) do |sum, sigma|
sum += sigma.zip(r).inject(1){|prod, (row, col)| prod *= self[row, col] }
end
end
end
m1 = Matrix[[1,2],[3,4]] # testcases from Python version
m2 = Matrix[[1, 2, 3, 4], [4, 5, 6, 7], [7, 8, 9, 10], [10, 11, 12, 13]]
m3 = Matrix[[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]]
[m1, m2, m3].each do |m|
puts "determinant:\t #{m.determinant}", "permanent:\t #{m.permanent}"
puts
end</lang>
- Output:
determinant: -2 permanent: 10 determinant: 0 permanent: 29556 determinant: 0 permanent: 6778800
Tcl
The determinant is provided by the linear algebra package in Tcllib. The permanent (being somewhat less common) requires definition, but is easily described:
<lang tcl>package require math::linearalgebra package 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]
}</lang> Demonstrating with a sample matrix: <lang tcl>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]</lang>
- Output:
1.1315223609263888e-29 29556