axiom-math
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Axiom-math] Pfaffian


From: Martin Rubey
Subject: [Axiom-math] Pfaffian
Date: 27 Sep 2007 22:00:21 +0200
User-agent: Gnus/5.09 (Gnus v5.9.0) Emacs/21.4

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)





reply via email to

[Prev in Thread] Current Thread [Next in Thread]