exact pfaffian bug

[stevengj]

Here is an example of an 8x8 matrix for which our (my) "exact" Pfaffian routine is apparently giving the wrong answer:

julia> A = [0 14 7 -10 0 10 0 -11; -14 0 -10 7 13 -9 -12 -13; -7 10 0 -4 6 -17 -1 18; 10 -7 4 0 -2 -4 0 11; 0 -13 -6 2 0 -8 -18 17; -10 9 17 4 8 0 -8 12; 0 12 1 0 18 8 0 0; 11 13 -18 -11 -17 -12 0 0]
8×8 Matrix{Int64}:
   0   14    7  -10    0   10    0  -11
 -14    0  -10    7   13   -9  -12  -13
  -7   10    0   -4    6  -17   -1   18
  10   -7    4    0   -2   -4    0   11
   0  -13   -6    2    0   -8  -18   17
 -10    9   17    4    8    0   -8   12
   0   12    1    0   18    8    0    0
  11   13  -18  -11  -17  -12    0    0

julia> using LinearAlgebra, SkewLinearAlgebra

julia> pfaffian(big.(A))^2
11244481600

julia> det(big.(A))
14161000000

julia> pfaffian(A)^2
1.4161e10

[stevengj]

In particular, it works for most 8x8 matrices, but occasionally it will fail:

julia> for i = 1:20
           A = skewhermitian(big.(rand(-10:10,8,8)*2))
           @show pfaffian(A)^2 == det(A)
       end
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = false
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true
pfaffian(A) ^ 2 == det(A) = true

So it must be in one of the rare branches of the code.

[stevengj]

Okay, it looks like there is a bug in the swapping code here: https://github.com/smataigne/SkewLinearAlgebra.jl/blob/8bb7028e4c460af4733c591a39adab02be85a36f/src/pfaffian.jl#L27-L34

The failure only occurs in the (relatively) rare cases where row/column swaps are needed.

[stevengj]

Yeah, it seems that I can't swap rows and columns in the same loop without some additional logic. Easier just to swap in separate loops, since this step is O(n) and hence negligible anyway.