diff --git a/M2/Macaulay2/d/interface.dd b/M2/Macaulay2/d/interface.dd index b25a030e876..5fd8782f82f 100644 --- a/M2/Macaulay2/d/interface.dd +++ b/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 diff --git a/M2/Macaulay2/e/interface/matrix.cpp b/M2/Macaulay2/e/interface/matrix.cpp index aac6c3d5b64..374543f6d55 100644 --- a/M2/Macaulay2/e/interface/matrix.cpp +++ b/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); diff --git a/M2/Macaulay2/e/interface/matrix.h b/M2/Macaulay2/e/interface/matrix.h index fe006c39c09..1ef6d1cc426 100644 --- a/M2/Macaulay2/e/interface/matrix.h +++ b/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 */ diff --git a/M2/Macaulay2/e/matrix.hpp b/M2/Macaulay2/e/matrix.hpp index 30a82163436..29016a1eadf 100644 --- a/M2/Macaulay2/e/matrix.hpp +++ b/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); diff --git a/M2/Macaulay2/e/pfaff.cpp b/M2/Macaulay2/e/pfaff.cpp index 6289036387b..6116a47cfcb 100644 --- a/M2/Macaulay2/e/pfaff.cpp +++ b/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 diff --git a/M2/Macaulay2/e/pfaff.hpp b/M2/Macaulay2/e/pfaff.hpp index fd9534e3e86..a9481cd3fa6 100644 --- a/M2/Macaulay2/e/pfaff.hpp +++ b/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 diff --git a/M2/Macaulay2/m2/exports.m2 b/M2/Macaulay2/m2/exports.m2 index 397d91e09a0..f0434fa4f46 100644 --- a/M2/Macaulay2/m2/exports.m2 +++ b/M2/Macaulay2/m2/exports.m2 @@ -1000,6 +1000,7 @@ export { "peek'", "permanents", "permutations", + "pfaffian", "pfaffians", "pi", "pivots", diff --git a/M2/Macaulay2/m2/multilin.m2 b/M2/Macaulay2/m2/multilin.m2 index 0621c4e6108..0462e7529d0 100644 --- a/M2/Macaulay2/m2/multilin.m2 +++ b/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 ----------------------------------------------------------------------------- diff --git a/M2/Macaulay2/packages/Macaulay2Doc/functions.m2 b/M2/Macaulay2/packages/Macaulay2Doc/functions.m2 index 7f760685903..0698e319e32 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/functions.m2 +++ b/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" diff --git a/M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffian-doc.m2 b/M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffian-doc.m2 new file mode 100644 index 00000000000..24903903b34 --- /dev/null +++ b/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 +/// diff --git a/M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffians-doc.m2 b/M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffians-doc.m2 index e09985ee536..b7e262c8f35 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffians-doc.m2 +++ b/M2/Macaulay2/packages/Macaulay2Doc/functions/pfaffians-doc.m2 @@ -45,5 +45,5 @@ document { "The skew symmetry of ", TT "M", " is not checked, but the algorithm proceeds as if it were, with somewhat unpredictable results!" }, - SeeAlso => {det, "matrices"} + SeeAlso => {pfaffian, det, "matrices"} } diff --git a/M2/Macaulay2/packages/Macaulay2Doc/ov_matrices.m2 b/M2/Macaulay2/packages/Macaulay2Doc/ov_matrices.m2 index b1ddffd91c6..eebde98ff60 100644 --- a/M2/Macaulay2/packages/Macaulay2Doc/ov_matrices.m2 +++ b/M2/Macaulay2/packages/Macaulay2Doc/ov_matrices.m2 @@ -26,7 +26,6 @@ document { "determinants and related computations", TO "rank of a matrix", TO "determinants and minors", - TO "Pfaffians", TO "exterior power of a matrix", "display of matrices and saving matrices to a file", TO "format and display of matrices in Macaulay2", @@ -635,16 +634,11 @@ document { TO determinant, TO permanents, TO pfaffians, + TO pfaffian, TO fittingIdeal, }, } -document { - Key => "Pfaffians", - - SUBSECTION "pfaffians" - } - document { Key => "exterior power of a matrix", "Since the ", TT "i", "-th exterior power is a functor, it applies diff --git a/M2/Macaulay2/tests/normal/pfaffians.m2 b/M2/Macaulay2/tests/normal/pfaffians.m2 index 5273361c44a..c61b19ae5b9 100644 --- a/M2/Macaulay2/tests/normal/pfaffians.m2 +++ b/M2/Macaulay2/tests/normal/pfaffians.m2 @@ -6,7 +6,9 @@ assert( pfaffians(1,m) == ideal(0_R) ) assert( pfaffians(2,m) == ideal(a,b,c,d,e,f) ) assert( pfaffians(3,m) == ideal(0_R) ) assert( pfaffians(4,m) == ideal(c*d-b*e+a*f) ) +assert( pfaffian m == c*d-b*e+a*f ) +assert( pfaffian matrix {} == 1 ) R = QQ[vars(0..9)] M = genericSkewMatrix(R,5)