Numerical calculation of pfaffian

Question

I tried to calculate numerical the pfaffian of a skew symmetric matrix by the recursive definition (from Wikipedia):

$$ \text{pf}\left(A\right) = \sum_{j=2}^{2n}\left(-1\right)^{j}a_{1j}\text{pf}\left(A_{\hat{1}\hat{j}}\right) $$

where $A_{\hat{1}\hat{j}}$ denotes the matrix A with both the i-th and j-th rows and columns removed and the pfaffian of the $0\times 0$ matrix is equal to one. My c++ is looking like

double pfaffian(arma::mat Test, int n)
{
     double pfa = 0.0f;

      if(n == 0)
      {
           pfa *= 1;
           return pfa;
      }
      else{
      for(int i = 1; i < n; i++)
      {
          arma::mat temp = Test;
          temp.shed_row(0);
          temp.shed_col(i);

          pfa +=std::pow((-1),i)*Test(0,i)*pfaffian(temp,n-1);
      }
   }

   return pfa;
   }

However, this code is not working for a simple matrix

$$\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} $$

I got zero and not -1. Also for bigger matrices it is not working.

Answer 1

You can try something like this

double pfaffian(arma::mat Test, int n)
{
  if(n == 0)
  {
       return 1.0f;
  }
  else if(n==1)
  { 
       return 0.0f;
  }
  else
  {
      double pfa = 0.0f;
      for(int i = 1; i < n; i++)
      {
          arma::mat temp = Test;

          temp.shed_col(i);
          temp.shed_row(i);
          temp.shed_row(0);
          if (n>2) temp.shed_col(0);
          pfa +=std::pow((-1),i+1)*Test(0,i)*pfaffian(temp,n-2);
      }
      return pfa;
   }
}