pfaffian

[smataigne]

first pfaffian.jl file

[smataigne]

It seems that LinearAlgebra.exactdiv creates problems for some tests of git

[stevengj]

You could just use div since we restricted this to integers.

But the most general thing is to do:

if isdefined(LA, :exactdiv) # julia ≥ 1.7
    const exactdiv = LA.exactdiv
else
    exactdiv(a, b) = a / b
    exactdiv(a::Integer, b::Integer) = div(a, b)
end

The exactdiv function was introduced in Julia 1.7 for JuliaLang/julia#40868, which added an exact determinant algorithm for arbitrary rings via the Bareiss algorithm. However, I'm not 100% sure that the algorithm here works for arbitrary rings, and in particular it might implicitly assume that * commutative. By default we should only use this for BigInt (even Int might overflow), but we can leave the function generic (but unexported) in case someone wants to call it (at their own risk), e.g. for Rational{BigInt}.

[stevengj]

I would also split it into multiple functions in order to support SkewHermitian matrices more efficiently:

function _exactpfaffian!(A::AbstractMatrix)
    ... same as above, but with no argument checks ...
end

function exactpfaffian!(A::AbstractMatrix)
    LinearAlgebra.require_one_based_indexing(A)
    isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-Hermitian matrix"))
    return _pfaffian!(A)
end
exactpfaffian!(A::SkewHermitian) = _pfaffian!(A.data)

pfaffian!(A::AbstractMatrix{<:BigInt}) = exactpfaffian!(A)

Note also that you will want an additional method somewhere for pfaffian(A::SkewHermTridiagonal{<:Real}) that just multiplies every other off-diagonal element.

[stevengj]

I would just call this pfaffian! and have it call hessenberg!

Then have additional methods

pfaffian(A::AbstractMatrix{<:Real}) = pfaffian!(copy(A))

function pfaffian!(A::AbstractMatrix{<:Real})
    isskewhermitian(A) || throw(...)
    return pfaffian!(SkewHermitian(A))
end

[stevengj]

change to commutative ring

[stevengj]

Now you can check

Abig = BigInt.(A)
@test pfaffian(A) ≈ pfaffian(Abig) == pfaffian(SkewHermitian(Abig))
@test pfaffian(Abig)^2 == det(Abig)

and this will call the pfaffian! methods etcetera.

[stevengj]

I would also try a couple of random integer matrices with another package for the Pfaffian, and check that you can reproduce those (just hardcode the matrices and the results into the test) to make sure that you get the right sign.

[stevengj]

I would also check an odd case.

[stevengj]

Suggested change

	A=rand(1:10,n,n)
	for i=1:n
	A[i,i]=0
	for j=i+1:n
	A[i,j]=-A[j,i]
	end
	end
	A = skewhermitian!(rand(-10:10,n,n)*2)

[stevengj]

Looks like this PR can be closed, I guess it was superseded by another PR #47?

In general, instead of opening new PRs, you should update an existing PR with git push (you can do a git push --force if the branch was rebased). That way the discussion history is preserved and associated with the PR that was merged.