/usr/share/axiom-20170501/src/algebra/PERMAN.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 | )abbrev package PERMAN Permanent
++ Authors: Johannes Grabmeier, Oswald Gschnitzer
++ Date Created: 7 August 1989
++ Date Last Updated: 23 August 1990
++ References:
++ Henryk Minc: Evaluation of Permanents,
++ Proc. of the Edinburgh Math. Soc.(1979), 22/1 pp 27-32.
++ Nijenhuis and Wilf : Combinatorical Algorithms, Academic
++ Press, New York 1978.
++ S.G.Williamson, Combinatorics for Computer Science,
++ Computer Science Press, 1985.
++ Description:
++ Permanent implements the functions permanent, the
++ permanent for square matrices.
Permanent(n,R) : SIG == CODE where
n : PositiveInteger
R : Ring with commutative("*")
I ==> Integer
L ==> List
V ==> Vector
SM ==> SquareMatrix(n,R)
VECTPKG1 ==> VectorPackage1(I)
NNI ==> NonNegativeInteger
PI ==> PositiveInteger
GRAY ==> GrayCode
SIG ==> with
permanent : SM -> R
++ permanent(x) computes the permanent of a square matrix x.
++ The permanent is equivalent to
++ the \spadfun{determinant} except that coefficients have
++ no change of sign. This function
++ is much more difficult to compute than the
++ determinant. The formula used is by H.J. Ryser,
++ improved by [Nijenhuis and Wilf, Ch. 19].
++ Note that permanent(x) choose one of three algorithms, depending
++ on the underlying ring R and on n, the number of rows (and
++ columns) of x:\br
++ if 2 has an inverse in R we can use the algorithm of
++ [Nijenhuis and Wilf, ch.19,p.158]; if 2 has no inverse,
++ some modifications are necessary:\br
++ if n > 6 and R is an integral domain with characteristic
++ different from 2 (the algorithm works if and only 2 is not a
++ zero-divisor of R and characteristic()$R ^= 2,
++ but how to check that for any given R ?),
++ the local function permanent2 is called;\br
++ else, the local function permanent3 is called
++ (works for all commutative rings R).
CODE ==> add
-- local functions:
permanent2: SM -> R
permanent3: SM -> R
x : SM
a,b : R
i,j,k,l : I
permanent3(x) ==
-- This algorithm is based upon the principle of inclusion-
-- exclusion. A Gray-code is used to generate the subsets of
-- 1,... ,n. This reduces the number of additions needed in
-- every step.
sgn : R := 1
k : R
a := 0$R
vv : V V I := firstSubsetGray(n)$GRAY
-- For the meaning of the elements of vv, see GRAY.
w : V R := new(n,0$R)
j := 1 -- Will be the number of the element changed in subset
while j ^= (n+1) repeat -- we sum over all subsets of (1,...,n)
sgn := -sgn
b := sgn
if vv.1.j = 1 then k := -1
else k := 1 -- was that element deleted(k=-1) or added(k=1)?
for i in 1..(n::I) repeat
w.i := w.i +$R k *$R x(i,j)
b := b *$R w.i
a := a +$R b
vv := nextSubsetGray(vv,n)$GRAY
j := vv.2.1
if odd?(n) then a := -a
a
permanent(x) ==
-- If 2 has an inverse in R, we can spare half of the calcu-
-- lation needed in "permanent3": This is the algorithm of
-- [Nijenhuis and Wilf, ch.19,p.158]
n = 1 => x(1,1)
two : R := (2:I) :: R
half : Union(R,"failed") := recip(two)
if (half case "failed") then
if n < 7 then return permanent3(x)
else return permanent2(x)
sgn : R := 1
a := 0$R
w : V R := new(n,0$R)
-- w.i will be at first x.i and later lambda.i in
-- [Nijenhuis and Wilf, p.158, (24a) resp.(26)].
rowi : V R := new(n,0$R)
for i in 1..n repeat
rowi := row(x,i) :: V R
b := 0$R
for j in 1..n repeat
b := b + rowi.j
w.i := rowi(n) - (half*b)$R
vv : V V I := firstSubsetGray((n-1): PI)$GRAY
-- For the meaning of the elements of vv, see GRAY.
n :: I
b := 1
for i in 1..n repeat
b := b * w.i
a := a+b
j := 1 -- Will be the number of the element changed in subset
while j ^= n repeat -- we sum over all subsets of (1,...,n-1)
sgn := -sgn
b := sgn
if vv.1.j = 1 then k := -1
else k := 1 -- was that element deleted(k=-1) or added(k=1)?
for i in 1..n repeat
w.i := w.i +$R k *$R x(i,j)
b := b *$R w.i
a := a +$R b
vv := nextSubsetGray(vv,(n-1) : PI)$GRAY
j := vv.2.1
if not odd?(n) then a := -a
two * a
permanent2(x) ==
c : R := 0
sgn : R := 1
if (not (R has IntegralDomain))
-- or (characteristic()$R = (2:NNI))
-- compiler refuses to compile the line above !!
or (sgn + sgn = c)
then return permanent3(x)
-- This is a slight modification of permanent which is
-- necessary if 2 is not zero or a zero-divisor in R, but has
-- no inverse in R.
n = 1 => x(1,1)
two : R := (2:I) :: R
a := 0$R
w : V R := new(n,0$R)
-- w.i will be at first x.i and later lambda.i in
-- [Nijenhuis and Wilf, p.158, (24a) resp.(26)].
rowi : V R := new(n,0$R)
for i in 1..n repeat
rowi := row(x,i) :: V R
b := 0$R
for j in 1..n repeat
b := b + rowi.j
w.i := (two*(rowi(n)))$R - b
vv : V V I := firstSubsetGray((n-1): PI)$GRAY
n :: I
b := 1
for i in 1..n repeat
b := b *$R w.i
a := a +$R b
j := 1 -- Will be the number of the element changed in subset
while j ^= n repeat -- we sum over all subsets of (1,...,n-1)
sgn := -sgn
b := sgn
if vv.1.j = 1 then k := -1
else k := 1 -- was that element deleted(k=-1) or added(k=1)?
c := k * two
for i in 1..n repeat
w.i := w.i +$R c *$R x(i,j)
b := b *$R w.i
a := a +$R b
vv := nextSubsetGray(vv,(n-1) : PI)$GRAY
j := vv.2.1
if not odd?(n) then a := -a
b := two ** ((n-1):NNI)
(a exquo b) :: R
|