[Axiom-math] Pfaffian

[Martin Rubey]

I have hacked together an algorithm to compute a Pfaffian, using an algorithm
of Günter Rote. Currently it's only an .input script, but if it's useful for
somebody else than myself, we could make it a little more proffessional.

Martin


B0 n == matrix [[(if i=j+1 and odd? j then -1 else if i=j-1 and odd? i then 1 _
else 0) for j in 1..n] for i in 1..n]

PfChar(lambda, A) ==
n := nrows A
(n = 2) => lambda^2 + A.(1,2)
M := subMatrix(A, 3, n, 3, n)
r := subMatrix(A, 1, 1, 3, n)
s := subMatrix(A, 3, n, 2, 2)

p := PfChar(lambda, M)
d := degree(p, lambda)

B := B0(n-2)
C := r*B
g := [(C*s).(1,1), A.(1,2), 1]
if d >= 4 then
B := M*B
for i in 4..d by 2 repeat
C := C*B
g := cons((C*s).(1,1), g)
g := reverse! g

res := 0
for i in 0..d by 2 for j in 2..d+2 repeat
c := coefficient(p, lambda, i)
for e in first(g, j) for k in 2..-d by -2 repeat
res := res + c * e * lambda^(k+i)

res

pfaffian A == eval(PfChar(l, A), l=0)

[root]

excellent. i've added it to the test queue.
can you give an input example or two? -- Tim