Rosetta Code had a new
task that had this expression for calculating the determinant of a square matrix (Leibniz's formula):
You can read the
wikipedia article to find out more about sigma and sgn or
this article which took me step-by-step through to understanding it.
Although the Leibniz formula is inefficient for calculation in general I decided to use it, which meant I needed to create permutations of n items,
σ, and the sign of each permutation
Swapping signs
The sign of a permutation is defined as
the parity of the number of swaps between any two items of a permutation needed to generate that permutation from the initial permutation.
So for example, given a starting perm of [0,1,2] the permutation [2,0,1] can be achieved by first swapping the 2 and the 0 to form [2,1,0], then swapping the 1 and 0 to form [2,0,1. That is two swaps in total. The sign of the permutation is given as -1 if an odd number of swaps are used in generating it otherwise as +1 if an even (or no) swaps are needed. (The sign of the permutation will not change even if other swaps, or redundant swaps are used to generate it from the same starting perm).
Pythons itertools.permutations function was of unknown sign, but I picked up on the
Steinhaus–Johnson–Trotter algorithm that generates permutations in such a way that it guarantees that successive permutations only differ by the swapping of two items - this meant that the sign would start as +1 and alternate in sign between successive items.
The Steinhaus–Johnson–Trotter Python program
I decided to create a new Rosetta code task for generating both permutations and their sign with a view for their subsequent use in calculating determinants:
from operator import itemgetter
DEBUG = False # like the built-in __debug__
def spermutations(n):
"""permutations by swapping. Yields: perm, sign"""
sign = 1
p = [[i, 0 if i == 0 else -1] # [num, direction]
for i in range(n)]
if DEBUG: print ' #', p
yield tuple(pp[0] for pp in p), sign
while any(pp[1] for pp in p): # moving
i1, (n1, d1) = max(((i, pp) for i, pp in enumerate(p) if pp[1]),
key=itemgetter(1))
sign *= -1
if d1 == -1:
# Swap down
i2 = i1 - 1
p[i1], p[i2] = p[i2], p[i1]
# If this causes the chosen element to reach the First or last
# position within the permutation, or if the next element in the
# same direction is larger than the chosen element:
if i2 == 0 or p[i2 - 1][0] > n1:
# The direction of the chosen element is set to zero
p[i2][1] = 0
elif d1 == 1:
# Swap up
i2 = i1 + 1
p[i1], p[i2] = p[i2], p[i1]
# If this causes the chosen element to reach the first or Last
# position within the permutation, or if the next element in the
# same direction is larger than the chosen element:
if i2 == n - 1 or p[i2 + 1][0] > n1:
# The direction of the chosen element is set to zero
p[i2][1] = 0
if DEBUG: print ' #', p
yield tuple(pp[0] for pp in p), sign
for i3, pp in enumerate(p):
n3, d3 = pp
if n3 > n1:
pp[1] = 1 if i3 < i2 else -1
if DEBUG: print ' # Set Moving'
if __name__ == '__main__':
from itertools import permutations
for n in (3, 4):
print '\nPermutations and sign of %i items' % n
sp = set()
for i in spermutations(n):
sp.add(i[0])
print('Perm: %r Sign: %2i' % i)
#if DEBUG: raw_input('?')
# Test
p = set(permutations(range(n)))
assert sp == p, 'Two methods of generating permutations do not agree'
Determinants using Leibniz's formula
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 = list(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 = list(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)))
Sign() of the times
I was still thinking about how to generate the sign of Pythons in-built permutation generator, and a little thought made me think that somewhere there was a sort routine that sorted by swapping two elements. If you count the swaps then you can calculate the sign if the number of swaps is even or odd.
def selectionSort(lst):
for i in range(0,len(lst)-1):
mn = min(range(i,len(lst)), key=lst.__getitem__)
lst[i],lst[mn] = lst[mn],lst[i]
return lst
Unfortunately the sorter does unnecessary swaps so I made a slight modification to create a routine that generated a sign of any permutation of the integers 0..n:
def perm_parity(lst):
'''\
Given a permutation of the digits 0..N in order as a list,
returns its parity (or sign): +1 for even parity; -1 for odd.
'''
parity = 1
for i in range(0,len(lst)-1):
if lst[i] != i:
parity *= -1
mn = min(range(i,len(lst)), key=lst.__getitem__)
lst[i],lst[mn] = lst[mn],lst[i]
return parity
if __name__ == '__main__':
from itertools import permutations
for p in permutations(range(3)):
l = list(p)
print "%2i %r" % (perm_parity(l), p)
I published the sign generator on ActiveState as I didn't think it was enough for a Rosetta Code task.
(If you run the code you can see that their is in fact a pattern to the signs generated by successive members Pythons permutations).
Patterns
If you look at a few of the series generated by the Steinhaus–Johnson–Trotter Python program you can see patterns that point to a recursive method of generating the same sequences. I coded that up too.
Serendipity doo dah
A few days later I came across a Stack-exchange question that needed to generate all permutations when there were duplicate items in the initial list. (It goes on to need the permutations in a specific order, but I wasn't thinking of that bit at first).
Of coure I started by generating possible perm then using a set to remove duplicates, but it gets very (very) wasteful, so I started looking around for an algorithm to generate what are sometimes called lexicographically correct permutations. for example the l-perms of [1,1,1,2,2,2] are the twenty terms:
1 1 1 2 2 2
1 1 2 1 2 2
1 2 1 1 2 2
2 1 1 1 2 2
2 1 1 2 1 2
1 2 1 2 1 2
1 1 2 2 1 2
1 2 2 1 1 2
2 1 2 1 1 2
2 2 1 1 1 2
2 2 1 1 2 1
2 1 2 1 2 1
1 2 2 1 2 1
1 1 2 2 2 1
1 2 1 2 2 1
2 1 1 2 2 1
2 1 2 2 1 1
1 2 2 2 1 1
2 2 1 2 1 1
2 2 2 1 1 1
Which is a lot less than 6! = 720 terms of the standard permutation.
I spent an hour playing with my recursive permutation generator until I hit on the idea that if you had to do permutations of 123XX then you could think of each (identical) X as being X0 and X1, and that valid l-perms are the perms when X0 and X1 appear in only one order.
Now the recursive algorithm for perms introduces each item of the initial perm one-by-one. If you restrict each X as it is introduced to never be introduced 'above' and X that is already in the partial perm item that is being expanded, then, in effect, X0 and X1 can never 'cross' and so will stay in the same 'order'.
I made the change and it seems to work for me.
Lexicographic permutation generator in Python
def l_perm(items):
if not items:
return [[]]
else:
dir = 1
new_items = []
this = [items.pop()]
for item in l_perm(items):
lenitem = len(item)
try:
# Never insert 'this' above any other 'this' in the item
maxinsert = item.index(this[0])
except ValueError:
maxinsert = lenitem
if dir == 1:
# step down
new_items += [item[:i] + this + item[i:] for i in range(lenitem, -1, -1)
if i <= maxinsert]
else:
# step up
new_items += [item[:i] + this + item[i:] for i in range(lenitem + 1)
if i <= maxinsert]
dir *= -1
return new_items
if __name__ == '__main__':
n = [1, 1, 1, 2, 2, 2]
print '\nLexicograpic Permutations of %i items: %r' % (len(n), n)
for i, x in enumerate(l_perm(n)):
print('%3i %r' % (i, x))
End bit.
I'm tired and I want to go to bed so turning l-perm into a RC task or sticking it on Activestate can wait. I am wondering if their is anyone else interested in this as I just pottered around the less deep bits.
- Paddy.
