snippet: view plain - save this
1 #! /usr/bin/env python
2 #
3 """Various algorithms for computing the sign of a permutation.
4
5 Ported from http://perplexus.info/show.php?pid=4895&op=sol
6 """
7 __docformat__ = 'reStructuredText'
8
9
10 ## perm_sign by John Burkardt is an example that can be found online
11 ## (see it among the collection at
12 ## http://orion.math.iastate.edu/burkardt/f_src/subset/subset.f90 or
13 ## see http://www.scs.fsu.edu/~burkardt/math2071/perm_sign.m). A
14 ## similar method, that, however, avoids explicit searching, is
15 ## illustrated by the following function::
16
17 def jbsign(p):
18 """Compute sign of permutation `p` by counting the number of
19 interchanges required to change the given permutation into the
20 identity one.
21
22 Examples::
23 >>> jbsign([0,1,2])
24 1
25 >>> jbsign([0,2,1])
26 -1
27 """
28 n = len(p)
29 s = +1
30 # get elements back in their home positions
31 for j in xrange(0,n):
32 q=p[j]
33 if q !=j:
34 p[j],p[q] = p[q],q # interchange p[j] and p[p[j]]
35 s = -s # and account for the interchange
36 # note that q is now in its home position
37 # whether or not an interchange was required
38 return s
39
40
41 ## A somewhat more sophisticated method follows the orbits of elements
42 ## as they are permuted through cycles::
43
44 def sign_by_orbits(p):
45 """Compute sign of permutation `p` by following orbit cycles.
46
47 Examples::
48 >>> sign_by_orbits([0,1,2])
49 1
50 >>> sign_by_orbits([0,2,1])
51 -1
52 """
53 n = len(p)
54 s = +1
55 # follow orbit cycles and tote up signs
56 for j in xrange(0,n):
57 if (p[j] >= 0):
58 # new cycle, follow it, marking each place as "visited"
59 q = j
60 t = +1
61 while (q >= 0):
62 r = p[q]
63 if (r >= 0):
64 p[q] = -1
65 t = -t
66 q = r
67 # factor in contribution of each cycle
68 s *= t
69 return s
70
71
72 ## Lang's Algebra, First Ed. 1965, p. 53, gives an equivalent
73 ## definition of the sign as simply `(-1)^m` where m=n-number of
74 ## orbits. Thus the above can be modified to just count the number of
75 ## orbits and do a little arithmetic at the end, as the following
76 ## function illustrates::
77
78 def sign_lang_formula(p):
79 """Compute sign of permutation `p` as `(-1)^{length of permutation
80 - number of orbits}`.
81
82 Examples::
83 >>> sign_lang_formula([0,1,2])
84 1
85 >>> sign_lang_formula([0,2,1])
86 -1
87 """
88 n = len(p)
89 c = 0
90 # count the orbit cycles
91 for j in xrange(0,n):
92 if (p[j] >= 0):
93 # new cycle, follow it, marking each place
94 q = j
95 while (q >= 0):
96 r = p[q]
97 if (r >= 0):
98 p[q] = -1
99 q = r
100 # increment cycle count
101 c += 1
102 if (n-c) % 2 == 0:
103 return +1
104 else:
105 return -1
106
107
108 ## Another method counts inversions and leads to the following short
109 ## program::
110
111 def sign_counting_inversions(p):
112 """Compute sign of permutation `p` by counting the number of
113 interchanges required to change the given permutation into the
114 identity one.
115
116 Examples::
117 >>> sign_counting_inversions([0,1,2])
118 1
119 >>> sign_counting_inversions([0,2,1])
120 -1
121 """
122 n = len(p)
123 v = 0
124 # count inversions
125 for i in xrange(1,n):
126 for j in xrange(i+1, n):
127 if (p[i] > p[j]):
128 v += 1
129 if (v % 2 == 0):
130 return +1
131 else:
132 return -1
133
134
135 ## And one more: Those with a sufficient background in algebra will
136 ## realize that a simple algorithm exists based on the Laplace
137 ## expansion of the determinant of a permutation matrix. However,
138 ## while rather similar to the algorithm above that counts inversions,
139 ## it turned out to be a bit more complicated to program::
140
141 def lpsign(p):
142 """Compute sign of permutation `p` via Laplace expansion of the
143 determinant of the permutation matrix.
144
145 Examples::
146 >>> lpsign([0,1,2])
147 1
148 >>> lpsign([0,2,1])
149 -1
150 """
151 # Laplace expansion of determinant of permutation
152 # matrix with a 1 in in row p[j] of column j
153 n = len(p)
154 s = 1
155 for j in xrange(1, n):
156 q=p[j]
157 # apply checkerboard sign
158 if (q % 2 == 0):
159 s = -s
160 for k in xrange(j+1, n):
161 r = p[k]
162 # adjust numbers of remaining rows beyond row p[j]
163 # to account for having deleted row p[j]
164 if (r > q):
165 p[k] = r-1
166 return s
167
168
169 def sign_recursive(p):
170 """Recursively compute the sign of permutation `p`.
171
172 A permutation is represented as a linear list: `p` maps
173 `i` to `p[i]`.
174
175 Examples:
176 >>> sign_recursive([0,1,2])
177 1
178 >>> sign_recursive([0,2,1])
179 -1
180 """
181 n = len(p)
182 if 1 >= n:
183 return 1
184 # find highest-numbered element
185 k = p.index(n-1)
186 # remove highest-numbered element for recursion
187 del p[k]
188 # recursively compute
189 if 0 == ((n+k) % 2):
190 s = -1
191 else:
192 s = 1
193 return s * sign_recursive(p)
194
195
196 ## So, based on simplicity, counting inversions appears to have the
197 ## edge over all the other algorithms discussed above. However, based
198 ## on speed, efficiently returning elements to their homes, and
199 ## counting orbits both seem to be quite fast.
200 ##
201 ## Indeed, running 100000 iterations of each function to compute the
202 ## sign of a 5-element permutation, gives::
203 ##
204 ## | $ python2.4 profile.py ./permsign.py
205 ## | ncalls tottime percall cumtime percall function
206 ## | [...]
207 ## | 100002 1.950 0.000 2.480 0.000 jbsign
208 ## | 100002 2.550 0.000 3.140 0.000 sign_counting_inversions
209 ## | 100002 2.810 0.000 3.350 0.000 sign_lang_formula
210 ## | 100002 2.820 0.000 3.510 0.000 sign_by_orbits
211 ## | 100002 3.470 0.000 4.190 0.000 lpsign
212 ## | 100002 11.450 0.000 16.890 0.000 sign_recursive
213 ## | [...]
214 ##
215
216
217
218 ## main: run tests
219
220 if "__main__" == __name__:
221 import doctest
222 doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
0 comments