Add "pfaffian" for finding Pfaffians of a skew symmetric matrices

We already have "pfaffians" for returning an ideal containing all the
sub-Pfaffians of a particular size, but we may just want the Pfaffian
itself as a number or polynomial, cf. "minors" v. "det".

Author: d-torrance

M2/Macaulay2/d/interface.dd

@@ -2262,6 +2262,13 @@ export rawPfaffians(e:Expr):Expr := (
      );
 setupfun("rawPfaffians",rawPfaffians);
 
+export rawPfaffian(e:Expr):Expr := (
+    when e
+    is M:RawMatrixCell do toExpr(Ccode(RawRingElementOrNull,
+	    "IM2_Matrix_pfaffian(", M.p, ")"))
+    else WrongArgMatrix());
+setupfun("rawPfaffian", rawPfaffian);
+
 export rawTensor(e:Expr):Expr := (
      when e is s:Sequence do
      when s.0

M2/Macaulay2/e/interface/matrix.cpp

@@ -556,6 +556,11 @@ const Matrix /* or null */ *IM2_Matrix_pfaffians(int p, const Matrix *M)
   return M->pfaffians(p);
 }
 
+const RingElement /* or null */ *IM2_Matrix_pfaffian(const Matrix *M)
+{
+  return RingElement::make_raw(M->get_ring(), M->pfaffian());
+}
+
 const Matrix /* or null */ *IM2_Matrix_diff(const Matrix *M, const Matrix *N)
 {
   return M->diff(N, 1);

M2/Macaulay2/e/interface/matrix.h

@@ -241,6 +241,9 @@ const Matrix /* or null */ *IM2_Matrix_pfaffians(
     int p,
     const Matrix *M); /* drg: connected rawPfaffians*/
 
+const RingElement /* or null */ *IM2_Matrix_pfaffian(
+    const Matrix *M);
+
 const Matrix *rawMatrixCompress(
     const Matrix *M); /* connected rawMatrixCompress */
 

M2/Macaulay2/e/matrix.hpp

@@ -165,7 +165,10 @@ class Matrix : public EngineObject
                                M2_arrayintOrNull first_col   // possibly NULL
                                ) const;
 
-  Matrix *pfaffians(int p) const;  // in pfaff.cpp
+  // in pfaff.cpp
+  Matrix *pfaffians(int p) const;
+  ring_elem pfaffian() const;
+
   static Matrix *wedge_product(int p, int q, const FreeModule *F);
   //  static Matrix wedge_dual(int p, const FreeModule *F);
 

M2/Macaulay2/e/pfaff.cpp

@@ -30,7 +30,7 @@ int PfaffianComputation::step()
 {
   if (done) return COMP_DONE;
 
-  ring_elem r = calc_pfaff(row_set, p);
+  ring_elem r = calc_pfaff();
   if (!R->is_zero(r))
     pfaffs.append(R->make_vec(0, r));
   else
@@ -54,6 +54,11 @@ int PfaffianComputation::calc(int nsteps)
     }
 }
 
+ring_elem PfaffianComputation::calc_pfaff(void)
+{
+  return calc_pfaff(row_set, p);
+}
+
 ring_elem PfaffianComputation::calc_pfaff(size_t *r, int p2)
 // Compute the pfaffian of the (skew symmetric)
 // minor with rows and columns r[0]..r[p2-1].
@@ -105,6 +110,13 @@ Matrix *Matrix::pfaffians(int p) const
   return d.pfaffians();
 }
 
+ring_elem Matrix::pfaffian() const
+{
+  PfaffianComputation d {this, n_cols()};
+
+  return d.calc_pfaff();
+}
+
 // Local Variables:
 // compile-command: "make -C $M2BUILDDIR/Macaulay2/e "
 // indent-tabs-mode: nil

M2/Macaulay2/e/pfaff.hpp

@@ -42,6 +42,7 @@ class PfaffianComputation : public our_new_delete
 
   Matrix *pfaffians() { return pfaffs.to_matrix(); }
   const Ring *get_ring() const { return R; }
+  ring_elem calc_pfaff(void);
 };
 
 #endif

M2/Macaulay2/m2/exports.m2

@@ -1000,6 +1000,7 @@ export {
 	"peek'",
 	"permanents",
 	"permutations",
+	"pfaffian",
 	"pfaffians",
 	"pi",
 	"pivots",

M2/Macaulay2/m2/multilin.m2

@@ -157,6 +157,16 @@ pfaffians = method(TypicalValue => Ideal)
 pfaffians(ZZ, Matrix) := (j, m) -> (
      ideal(map(ring m, rawPfaffians(j,raw m))))
 
+pfaffian = method()
+pfaffian Matrix := M -> (
+    if (
+	numRows M != numColumns M or
+	not zero(M + transpose M))
+    then error "expected a skew symmetric matrix";
+    if odd numRows M then 0_(ring M)
+    else if numRows M == 0 then 1_(ring M)
+    else promote(rawPfaffian raw M, ring M))
+
 -----------------------------------------------------------------------------
 -- trace and determinant
 -----------------------------------------------------------------------------

M2/Macaulay2/packages/Macaulay2Doc/functions.m2

@@ -203,6 +203,7 @@ load "./functions/pdim-doc.m2"
 load "./functions/peek-doc.m2"
 load "./functions/permanents-doc.m2"
 load "./functions/permutations-doc.m2"
+load "./functions/pfaffian-doc.m2"
 load "./functions/pfaffians-doc.m2"
 load "./functions/pivots-doc.m2"
 load "./functions/poincare-doc.m2"

M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffian-doc.m2

@@ -0,0 +1,35 @@
+doc ///
+  Key
+     pfaffian
+    (pfaffian, Matrix)
+  Headline
+    Pfaffian of a skew-symmetric matrix
+  Usage
+    pfaffian M
+  Inputs
+    M:Matrix -- must be skew-symmetric
+  Outputs
+    :{Number,RingElement} -- the Pfaffian of @VAR "M"@
+  Description
+    Text
+      The @wikipedia "Pfaffian"@ of a $2n\times 2n$ skew-symmetric matrix
+      $M=(m_{ij})$ is
+      $$\frac{1}{2^nn!}\sum_{\sigma\in S_{2n}}\operatorname{sgn}(\sigma)\prod_{i=1}^nm_{\sigma(2i-1),\sigma(2i)},$$
+      where $S_n$ is the symmetric group on $2n$ elements.
+    Example
+      M = matrix {{0, 1, 2, 3}, {-1, 0, 4, 5}, {-2, -4, 0, 6}, {-3, -5, -6, 0}}
+      pfaffian M
+    Text
+      The square of the Pfaffian is the determinant.
+    Example
+      determinant M
+    Text
+      Skew-symmetric matrices with an odd number of rows and columns have
+      Pfaffian zero.
+    Example
+      M = matrix {{0, 1, 2}, {-1, 0, 3}, {-2, -3, 0}}
+      pfaffian M
+ SeeAlso
+   determinant
+   pfaffians
+///