Ticket #15245: trac_15245-pfaffian-dg.patch

File trac_15245-pfaffian-dg.patch, 25.7 KB (added by Darij Grinberg, 9 years ago)

new version, separating skew-symmetry from alternatingness systematically

  • sage/combinat/perfect_matching.py

    # HG changeset patch
    # User darij grinberg <darijgrinberg@gmail.com>
    # Date 1381097450 25200
    # Node ID 9af1470d1868007c3ddc6444403211f7e86ea351
    # Parent  53809c21d659953c6ae884df767a4e3c237aa41f
    trac #15245: stupid-but-still-implementation of Pfaffian
    
    diff --git a/sage/combinat/perfect_matching.py b/sage/combinat/perfect_matching.py
    a b class PerfectMatching(ElementWrapper): 
    205205        # executed and we do not have an infinite loop.
    206206        return PerfectMatchings(objects)(data)
    207207
     208    def __iter__(self):
     209        r"""
     210        Iterate over the edges of the matching ``self``.
     211
     212        The edges are yielded as 2-tuples. Neither the elements of these
     213        tuples nor the tuples are necessarily sorted in any predictable
     214        way.
     215
     216        EXAMPLES::
     217
     218            sage: list(PerfectMatching([('a','e'),('b','c'),('d','f')]))
     219            [('a', 'e'), ('b', 'c'), ('d', 'f')]
     220            sage: list(PerfectMatchings(2)[0])
     221            [(1, 2)]
     222            sage: list(PerfectMatching([3,8,1,7,6,5,4,2]))
     223            [(1, 3), (2, 8), (4, 7), (5, 6)]
     224        """
     225        return self.value.__iter__()
     226
    208227    def _repr_(self):
    209228        r"""
    210         returns the name of the object
     229        Return a string representation of the matching ``self``.
    211230
    212231        EXAMPLES::
    213232
  • sage/matrix/matrix0.pyx

    diff --git a/sage/matrix/matrix0.pyx b/sage/matrix/matrix0.pyx
    a b cdef class Matrix(sage.structure.element 
    34543454
    34553455    def is_skew_symmetric(self):
    34563456        """
    3457         Returns True if this is a skew symmetric matrix.
     3457        Return ``True`` if ``self`` is a skew-symmetric matrix.
     3458
     3459        Here, "skew-symmetric matrix" means a square matrix `A`
     3460        satisfying `A^T = -A`. It does not require that the
     3461        diagonal entries of `A` are `0` (although this
     3462        automatically follows from `A^T = -A` when `2` is
     3463        invertible in the ground ring over which the matrix is
     3464        considered). Skew-symmetric matrices `A` whose diagonal
     3465        entries are `0` are said to be "alternating", and this
     3466        property is checked by the :meth:`is_alternating`
     3467        method.
    34583468
    34593469        EXAMPLES::
    34603470       
    cdef class Matrix(sage.structure.element 
    34643474            sage: m = matrix(QQ, [[1,2], [2,1]])
    34653475            sage: m.is_skew_symmetric()
    34663476            False
     3477
     3478        Skew-symmetric is not the same as alternating when
     3479        `2` is a zero-divisor in the ground ring::
     3480
     3481            sage: n = matrix(Zmod(4), [[0, 1], [-1, 2]])
     3482            sage: n.is_skew_symmetric()
     3483            True
     3484
     3485        but yet the diagonal cannot be completely
     3486        arbitrary in this case::
     3487
     3488            sage: n = matrix(Zmod(4), [[0, 1], [-1, 3]])
     3489            sage: n.is_skew_symmetric()
     3490            False
    34673491        """
    34683492        if self._ncols != self._nrows: return False
    34693493        # could be bigger than an int on a 64-bit platform, this
    cdef class Matrix(sage.structure.element 
    34763500                    return False
    34773501        return True
    34783502
     3503    def is_alternating(self):
     3504        """
     3505        Return ``True`` if ``self`` is an alternating matrix.
     3506
     3507        Here, "alternating matrix" means a square matrix `A`
     3508        satisfying `A^T = -A` and such that the diagonal entries
     3509        of `0`. Notice that the condition that the diagonal
     3510        entries be `0` is not redundant for matrices over
     3511        arbitrary ground rings (but it is redundant when `2` is
     3512        invertible in the ground ring). A square matrix `A` only
     3513        required to satisfy `A^T = -A` is said to be
     3514        "skew-symmetric", and this property is checked by the
     3515        :meth:`is_skew_symmetric` method.
     3516
     3517        EXAMPLES::
     3518       
     3519            sage: m = matrix(QQ, [[0,2], [-2,0]])
     3520            sage: m.is_alternating()
     3521            True
     3522            sage: m = matrix(QQ, [[1,2], [2,1]])
     3523            sage: m.is_alternating()
     3524            False
     3525
     3526        In contrast to the property of being skew-symmetric, the
     3527        property of being alternating does not tolerate nonzero
     3528        entries on the diagonal even if they are their own
     3529        negatives::
     3530
     3531            sage: n = matrix(Zmod(4), [[0, 1], [-1, 2]])
     3532            sage: n.is_alternating()
     3533            False
     3534        """
     3535        if self._ncols != self._nrows: return False
     3536        # could be bigger than an int on a 64-bit platform, this
     3537        #  is the type used for indexing.
     3538        cdef Py_ssize_t i,j
     3539
     3540        for i from 0 <= i < self._nrows:
     3541            for j from 0 <= j < i:
     3542                if self.get_unsafe(i,j) != -self.get_unsafe(j,i):
     3543                    return False
     3544            if self.get_unsafe(i,i) != 0:
     3545                return False
     3546        return True
     3547
    34793548    def is_symmetrizable(self, return_diag=False, positive=True):
    34803549        r"""
    34813550        This function takes a square matrix over an *ordered integral domain* and checks if it is symmetrizable.
  • sage/matrix/matrix2.pyx

    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b cdef class Matrix(matrix1.Matrix): 
    681681
    682682    def permanent(self):
    683683        r"""
    684         Calculate and return the permanent of this `m \times n`
    685         matrix using Ryser's algorithm.
    686        
     684        Calculate and return the permanent of the `m \times n`
     685        matrix ``self`` using Ryser's algorithm.
     686
    687687        Let `A = (a_{i,j})` be an `m \times n` matrix over
    688688        any commutative ring, with `m \le n`. The permanent of
    689689        `A` is
    690        
    691         .. math::
    692        
    693            \text{per}(A) = \sum_\pi a_{1,\pi(1)}a_{2,\pi(2)} \cdots a_{m,\pi(m)}         
    694        
    695        
     690
     691        .. MATH::
     692
     693           \mathrm{per}(A)
     694           = \sum_\pi a_{1,\pi(1)} a_{2,\pi(2)} \cdots a_{m,\pi(m)}         
     695
    696696        where the summation extends over all one-to-one functions
    697697        `\pi` from `\{1, \ldots, m\}` to
    698698        `\{1, \ldots, n\}`.
    699        
     699
    700700        The product
    701         `a_{1,\pi(1)}a_{2,\pi(2)} \cdots a_{m,\pi(m)}` is
     701        `a_{1,\pi(1)} a_{2,\pi(2)} \cdots a_{m,\pi(m)}` is
    702702        called diagonal product. So the permanent of an
    703703        `m \times n` matrix `A` is the sum of all the
    704704        diagonal products of `A`.
    705        
     705
     706        INPUT:
     707
     708        - ``A`` -- matrix of size `m \times n` with `m \leq n`
     709
     710        OUTPUT:
     711
     712        permanent of the matrix `A`
     713
     714        ALGORITHM:
     715
    706716        Modification of theorem 7.1.1. from Brualdi and Ryser:
    707717        Combinatorial Matrix Theory. Instead of deleting columns from
    708718        `A`, we choose columns from `A` and calculate the
    709719        product of the row sums of the selected submatrix.
    710        
    711         INPUT:
    712        
    713        
    714         -  ``A`` - matrix of size m x n with m = n
    715        
    716        
    717         OUTPUT: permanent of matrix A
    718        
    719         EXAMPLES::
    720        
     720
     721        EXAMPLES::
     722
    721723            sage: M = MatrixSpace(ZZ,4,4)
    722724            sage: A = M([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1])
    723725            sage: A.permanent()
    724726            24
    725        
    726         ::
    727        
     727
     728        ::
     729
    728730            sage: M = MatrixSpace(QQ,3,6)
    729731            sage: A = M([1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1])
    730732            sage: A.permanent()
    731733            36
    732        
    733         ::
    734        
     734
     735        ::
     736
    735737            sage: M = MatrixSpace(RR,3,6)
    736738            sage: A = M([1.0,1.0,1.0,1.0,0,0,0,1.0,1.0,1.0,1.0,0,0,0,1.0,1.0,1.0,1.0])
    737739            sage: A.permanent()
    738740            36.0000000000000
    739        
     741
    740742        See Sloane's sequence OEIS A079908(3) = 36, "The Dancing School
    741743        Problems"
    742        
    743         ::
    744        
     744
     745        ::
     746
    745747            sage: print sloane_sequence(79908)                # optional (internet connection)
    746748            Searching Sloane's online database...
    747749            [79908, 'Solution to the Dancing School Problem with 3 girls and n+3 boys: f(3,n).', [1, 4, 14, 36, 76, 140, 234, 364, 536, 756, 1030, 1364, 1764, 2236, 2786, 3420, 4144, 4964, 5886, 6916, 8060, 9324, 10714, 12236, 13896, 15700, 17654, 19764, 22036, 24476, 27090, 29884, 32864, 36036, 39406, 42980, 46764, 50764, 54986, 59436]]
    748        
    749         ::
    750        
     750
     751        ::
     752
    751753            sage: M = MatrixSpace(ZZ,4,5)
    752754            sage: A = M([1,1,0,1,1,0,1,1,1,1,1,0,1,0,1,1,1,0,1,0])
    753755            sage: A.permanent()
    754756            32
    755        
     757
    756758        See Minc: Permanents, Example 2.1, p. 5.
    757        
    758         ::
    759        
     759
     760        ::
     761
    760762            sage: M = MatrixSpace(QQ,2,2)
    761763            sage: A = M([1/5,2/7,3/2,4/5])
    762764            sage: A.permanent()
    763765            103/175
    764        
    765         ::
    766        
     766
     767        ::
     768
    767769            sage: R.<a> = PolynomialRing(ZZ)
    768770            sage: A = MatrixSpace(R,2)([[a,1], [a,a+1]])
    769771            sage: A.permanent()
    770772            a^2 + 2*a
    771        
    772         ::
    773        
     773
     774        ::
     775
    774776            sage: R.<x,y> = PolynomialRing(ZZ,2)
    775777            sage: A = MatrixSpace(R,2)([x, y, x^2, y^2])
    776778            sage: A.permanent()
    777779            x^2*y + x*y^2
    778        
     780
    779781        AUTHORS:
    780782
    781783        - Jaap Spies (2006-02-16)
    782        
     784
    783785        - Jaap Spies (2006-02-21): added definition of permanent
    784786        """
    785787        cdef Py_ssize_t m, n, r
    cdef class Matrix(matrix1.Matrix): 
    808810
    809811    def permanental_minor(self, Py_ssize_t k):
    810812        r"""
    811         Calculates the permanental `k`-minor of a
    812         `m \times n` matrix.
    813        
     813        Return the permanental `k`-minor of an `m \times n` matrix.
     814
    814815        This is the sum of the permanents of all possible `k` by
    815816        `k` submatrices of `A`.
    816        
     817
    817818        See Brualdi and Ryser: Combinatorial Matrix Theory, p. 203. Note
    818819        the typo `p_0(A) = 0` in that reference! For applications
    819820        see Theorem 7.2.1 and Theorem 7.2.4.
    820        
     821
    821822        Note that the permanental `m`-minor equals
    822         `per(A)`.
    823        
     823        `\mathrm{per}(A)` if `m = n`.
     824
    824825        For a (0,1)-matrix `A` the permanental `k`-minor
    825826        counts the number of different selections of `k` 1's of
    826         `A` with no two of the 1's on the same line.
    827        
     827        `A` with no two of the 1's on the same row and no two of the
     828        1's on the same column.
     829
    828830        INPUT:
    829        
    830        
    831         -  ``self`` - matrix of size m x n with m = n
    832        
    833        
    834         OUTPUT: permanental k-minor of matrix A
    835        
    836         EXAMPLES::
    837        
     831
     832        - ``self`` -- matrix of size `m \times n` with `m \leq n`
     833
     834        OUTPUT:
     835
     836        The permanental `k`-minor of the matrix ``self``.
     837
     838        EXAMPLES::
     839
    838840            sage: M = MatrixSpace(ZZ,4,4)
    839841            sage: A = M([1,0,1,0,1,0,1,0,1,0,10,10,1,0,1,1])
    840842            sage: A.permanental_minor(2)
    841843            114
    842        
    843         ::
    844        
     844
     845        ::
     846
    845847            sage: M = MatrixSpace(ZZ,3,6)
    846848            sage: A = M([1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1])
    847849            sage: A.permanental_minor(0)
    cdef class Matrix(matrix1.Matrix): 
    852854            40
    853855            sage: A.permanental_minor(3)
    854856            36
    855        
    856         Note that if k == m the permanental k-minor equals per(A)
    857        
    858         ::
    859        
     857
     858        Note that if `k = m = n`, the permanental `k`-minor equals
     859        `\mathrm{per}(A)`::
     860
    860861            sage: A.permanent()
    861862            36
    862        
    863         ::
    864        
     863
     864        ::
     865
    865866            sage: A.permanental_minor(5)
    866867            0
    867        
    868         For C the "complement" of A::
    869        
     868
     869        For `C` the "complement" of `A`::
     870
    870871            sage: M = MatrixSpace(ZZ,3,6)
    871872            sage: C = M([0,0,0,0,1,1,1,0,0,0,0,1,1,1,0,0,0,0])
    872873            sage: m, n = 3, 6
    873874            sage: sum([(-1)^k * C.permanental_minor(k)*factorial(n-k)/factorial(n-m) for k in range(m+1)])
    874875            36
    875        
    876         See Theorem 7.2.1 of Brualdi: and Ryser: Combinatorial Matrix
     876
     877        See Theorem 7.2.1 of Brualdi and Ryser: Combinatorial Matrix
    877878        Theory: per(A)
    878        
     879
    879880        AUTHORS:
    880881
    881882        - Jaap Spies (2006-02-19)
    cdef class Matrix(matrix1.Matrix): 
    883884        m = self._nrows
    884885        n = self._ncols
    885886        if not m <= n:
    886             raise ValueError, "must have m <= n, but m (=%s) and n (=%s)"%(m,n)
     887            raise ValueError("must have m <= n, but m (=%s) and n (=%s)"%(m,n))
    887888
    888889        R = self._base_ring
    889890        if k == 0:
    890             return R(1)
     891            return R.one()
    891892        if k > m:
    892             return R(0)
     893            return R.zero()
    893894
    894895        pm = 0
    895896        for cols in _choose(n,k):
    cdef class Matrix(matrix1.Matrix): 
    900901
    901902    def rook_vector(self, check = False):
    902903        r"""
    903         Returns rook vector of this matrix.
    904        
    905         Let `A` be a general `m` by `n`
    906         (0,1)-matrix with `m \le n`. We identify `A` with a
    907         chessboard where rooks can be placed on the fields corresponding
    908         with `a_{ij} = 1`. The number
    909         `r_k =  p_k(A)` (the permanental
    910         `k`-minor) counts the number of ways to place `k`
    911         rooks on this board so that no two rooks can attack another.
    912        
    913         The rook vector is the list consisting of
     904        Return the rook vector of the matrix ``self``.
     905
     906        Let `A` be an `m` by `n` (0,1)-matrix with `m \le n`. We identify
     907        `A` with a chessboard where rooks can be placed on the fields
     908        `(i, j)` with `a_{i,j} = 1`. The number
     909        `r_k = p_k(A)` (the permanental `k`-minor) counts the number of
     910        ways to place `k` rooks on this board so that no rook can attack
     911        another.
     912
     913        The rook vector of the matrix `A` is the list consisting of
    914914        `r_0, r_1, \ldots, r_m`.
    915        
     915
    916916        The rook polynomial is defined by
    917917        `r(x) = \sum_{k=0}^m r_k x^k`.
    918        
     918
    919919        INPUT:
    920        
    921        
    922         -  ``self`` - m by n matrix with m = n
    923        
    924         -  ``check`` - True or False (default), optional
    925        
    926        
    927         OUTPUT: rook vector
    928        
    929         EXAMPLES::
    930        
     920
     921        - ``self`` -- an `m` by `n` (0,1)-matrix with `m \le n`
     922
     923        - ``check`` -- Boolean (default: ``False``) determining whether
     924          to check that ``self`` is a (0,1)-matrix.
     925
     926        OUTPUT:
     927
     928        The rook vector of the matrix ``self``.
     929
     930        EXAMPLES::
     931
    931932            sage: M = MatrixSpace(ZZ,3,6)
    932933            sage: A = M([1,1,1,1,0,0,0,1,1,1,1,0,0,0,1,1,1,1])
    933934            sage: A.rook_vector()
    934935            [1, 12, 40, 36]
    935        
    936         ::
    937        
     936
     937        ::
     938
    938939            sage: R.<x> = PolynomialRing(ZZ)
    939940            sage: rv = A.rook_vector()
    940941            sage: rook_polynomial = sum([rv[k] * x^k for k in range(len(rv))])
    941942            sage: rook_polynomial
    942943            36*x^3 + 40*x^2 + 12*x + 1
    943        
     944
    944945        AUTHORS:
    945946
    946947        - Jaap Spies (2006-02-24)
    cdef class Matrix(matrix1.Matrix): 
    964965        return tmp
    965966
    966967    def minors(self,k):
    967         """
    968         Return the list of all k-minors of self.
    969 
    970         Let A be an m x n matrix and k an integer with 0 < k, k <= m, and
    971         k <= n. A k x k minor of A is the determinant of a k x k matrix
    972         obtained from A by deleting m - k rows and n - k columns.
     968        r"""
     969        Return the list of all `k \times k` minors of self.
     970
     971        Let `A` be an `m \times n` matrix and `k` an integer with
     972        `0 \leq k`, `k \leq m` and `k \leq n`.
     973        A `k \times k` minor of `A` is the determinant of a
     974        `k \times k` matrix obtained from `A` by deleting `m - k`
     975        rows and `n - k` columns.
     976        There are no `k \times k` minors of `A` if `k` is larger
     977        than either `m` or `n`.
    973978
    974979        The returned list is sorted in lexicographical row major ordering,
    975         e.g., if A is a 3 x 3 matrix then the minors returned are with for
    976         these rows/columns: [ [0, 1]x[0, 1], [0, 1]x[0, 2], [0, 1]x[1, 2],
     980        e.g., if A is a `3 \times 3` matrix then the minors returned are
     981        with these rows/columns: [ [0, 1]x[0, 1], [0, 1]x[0, 2], [0, 1]x[1, 2],
    977982        [0, 2]x[0, 1], [0, 2]x[0, 2], [0, 2]x[1, 2], [1, 2]x[0, 1], [1,
    978983        2]x[0, 2], [1, 2]x[1, 2] ].
    979        
     984
    980985        INPUT:
    981        
    982        
    983         -  ``k`` - integer
    984        
    985        
    986         EXAMPLE::
    987        
     986
     987        - ``k`` -- integer
     988
     989        EXAMPLES::
     990
    988991            sage: A = Matrix(ZZ,2,3,[1,2,3,4,5,6]); A
    989992            [1 2 3]
    990993            [4 5 6]
    991994            sage: A.minors(2)
    992995            [-3, -6, -3]
    993        
    994         ::
    995        
     996            sage: A.minors(1)
     997            [1, 2, 3, 4, 5, 6]
     998            sage: A.minors(0)
     999            [1]
     1000            sage: A.minors(5)
     1001            []
     1002
     1003        ::
     1004
    9961005            sage: k = GF(37)
    9971006            sage: P.<x0,x1,x2> = PolynomialRing(k)
    9981007            sage: A = Matrix(P,2,3,[x0*x1, x0, x1, x2, x2 + 16, x2 + 5*x1 ])
    cdef class Matrix(matrix1.Matrix): 
    12421251                self.swap_rows(level, i)
    12431252            return d       
    12441253
     1254    def pfaffian(self, algorithm=None, check=True):
     1255        r"""
     1256        Return the Pfaffian of ``self``, assuming that ``self`` is an
     1257        alternating matrix.
     1258
     1259        INPUT:
     1260
     1261        - ``algorithm`` -- string, the algorithm to use; currently the
     1262          following algorithms have been implemented:
     1263
     1264          * ``'definition'`` - using the definition given by perfect
     1265            matchings
     1266
     1267        - ``check`` (default: ``True``) -- Boolean determining whether to
     1268          check ``self`` for alternatingness and squareness. This has to
     1269          be set to ``False`` if ``self`` is defined over a non-discrete
     1270          ring.
     1271
     1272        The Pfaffian of an alternating matrix is defined as follows:
     1273
     1274        Let `A` be an alternating `k \times k` matrix over a commutative
     1275        ring. (Here, "alternating" means that `A^T = -A` and that the
     1276        diagonal entries of `A` are zero.)
     1277        If `k` is odd, then the Pfaffian of the matrix `A` is
     1278        defined to be `0`. Let us now define it when `k` is even. In this
     1279        case, set `n = k/2` (this is an integer). For every `i` and `j`,
     1280        we denote the `(i, j)`-th entry of `A` by `a_{i, j}`. Let `M`
     1281        denote the set of all perfect matchings of the set
     1282        `\{ 1, 2, \ldots, 2n \}` (see
     1283        :class:`sage.combinat.perfect_matching.PerfectMatchings` ).
     1284        For every matching `m \in M`, define the sign `\mathrm{sign}(m)`
     1285        of `m` by writing `m` as `\{ \{ i_1, j_1 \}, \{ i_2, j_2 \},
     1286        \ldots, \{ i_n, j_n \} \}` with `i_k < j_k` for all `k`, and
     1287        setting `\mathrm{sign}(m)` to be the sign of the permutation
     1288        `( i_1, j_1, i_2, j_2, \ldots, i_n, j_n )` (written here in
     1289        one-line notation). For every matching `m \in M`, define the
     1290        weight `w(m)` of `m` by writing `m` as
     1291        `\{ \{ i_1, j_1 \}, \{ i_2, j_2 \}, \ldots, \{ i_n, j_n \} \}`
     1292        with `i_k < j_k` for all `k`, and setting
     1293        `w(m) = a_{i_1, j_1} a_{i_2, j_2} \cdots a_{i_n, j_n}`. Now, the
     1294        Pfaffian of the matrix `A` is defined to be the sum
     1295
     1296        .. MATH::
     1297
     1298            \sum_{m \in M} \mathrm{sign}(m) w(m).
     1299
     1300        The Pfaffian of `A` is commonly denoted by `\mathrm{Pf}(A)`. It
     1301        is well-known that `(\mathrm{Pf}(A))^2 = \det A` for every
     1302        alternating matrix `A`, and that
     1303        `\mathrm{Pf} (U^T A U) = \det U \cdot \mathrm{Pf}(A)` for any
     1304        `n \times n` matrix `U` and any alternating `n \times n`
     1305        matrix `A`.
     1306
     1307        See [Kn95]_, [DW95]_ and [Rote2001]_, just to name three
     1308        sources, for further properties of Pfaffians.
     1309
     1310        ALGORITHM:
     1311
     1312        The current implementation uses the definition given above.
     1313        It checks alternatingness of the matrix ``self`` only if
     1314        ``check`` is ``True`` (this is important because even if ``self``
     1315        is alternating, a non-discrete base ring might prevent Sage
     1316        from being able to check this).
     1317
     1318        REFERENCES:
     1319
     1320        .. [Kn95] Donald E. Knuth, *Overlapping Pfaffians*,
     1321           :arxiv:`math/9503234v1`.
     1322
     1323        .. [Rote2001] Gunter Rote,
     1324           *Division-Free Algorithms for the Determinant and the
     1325           Pfaffian: Algebraic and Combinatorial Approaches*,
     1326           H. Alt (Ed.): Computational Discrete Mathematics, LNCS
     1327           2122, pp. 119–135, 2001.
     1328           http://page.mi.fu-berlin.de/rote/Papers/pdf/Division-free+algorithms.pdf
     1329
     1330        .. [DW95] Andreas W.M. Dress, Walter Wenzel,
     1331           *A Simple Proof of an Identity Concerning Pfaffians of
     1332           Skew Symmetric Matrices*,
     1333           Advances in Mathematics, volume 112, Issue 1, April
     1334           1995, pp. 120-134.
     1335           http://www.sciencedirect.com/science/article/pii/S0001870885710298
     1336
     1337        .. TODO::
     1338
     1339            Implement faster algorithms, including a division-free one.
     1340            Does [Rote2001]_, section 3.3 give one?
     1341
     1342            Check the implementation of the matchings used here for
     1343            performance?
     1344
     1345        EXAMPLES:
     1346
     1347        A `3 \times 3` alternating matrix has Pfaffian 0 independently
     1348        of its entries::
     1349
     1350            sage: MSp = MatrixSpace(Integers(27), 3)
     1351            sage: A = MSp([0, 2, -3,  -2, 0, 8,  3, -8, 0])
     1352            sage: A.pfaffian()
     1353            0
     1354            sage: parent(A.pfaffian())
     1355            Ring of integers modulo 27
     1356
     1357        The Pfaffian of a `2 \times 2` alternating matrix is just its
     1358        northeast entry::
     1359
     1360            sage: MSp = MatrixSpace(QQ, 2)
     1361            sage: A = MSp([0, 4,  -4, 0])
     1362            sage: A.pfaffian()
     1363            4
     1364            sage: parent(A.pfaffian())
     1365            Rational Field
     1366
     1367        The Pfaffian of a `0 \times 0` alternating matrix is `1`::
     1368
     1369            sage: MSp = MatrixSpace(ZZ, 0)
     1370            sage: A = MSp([])
     1371            sage: A.pfaffian()
     1372            1
     1373            sage: parent(A.pfaffian())
     1374            Integer Ring
     1375
     1376        Let us compute the Pfaffian of a generic `4 \times 4`
     1377        alternating matrix::
     1378       
     1379            sage: R = PolynomialRing(QQ, 'x12,x13,x14,x23,x24,x34')
     1380            sage: x12, x13, x14, x23, x24, x34 = R.gens()
     1381            sage: A = matrix(R, [[   0,  x12,  x13,  x14],
     1382            ....:                [-x12,    0,  x23,  x24],
     1383            ....:                [-x13, -x23,    0,  x34],
     1384            ....:                [-x14, -x24, -x34,    0]])
     1385            sage: A.pfaffian()
     1386            x14*x23 - x13*x24 + x12*x34
     1387            sage: parent(A.pfaffian())
     1388            Multivariate Polynomial Ring in x12, x13, x14, x23, x24, x34 over Rational Field
     1389
     1390        The Pfaffian of an alternating matrix squares to its
     1391        determinant::
     1392
     1393            sage: A = [[0] * 6 for i in range(6)]
     1394            sage: for i in range(6):
     1395            ....:     for j in range(i):
     1396            ....:         u = floor(random() * 10)
     1397            ....:         A[i][j] = u
     1398            ....:         A[j][i] = -u
     1399            ....:     A[i][i] = 0
     1400            sage: AA = Matrix(ZZ, A)
     1401            sage: AA.pfaffian() ** 2 == AA.det()
     1402            True
     1403
     1404        AUTHORS:
     1405
     1406        - Darij Grinberg (1 Oct 2013): first (slow)
     1407          implementation.
     1408        """
     1409        k = self._nrows
     1410
     1411        if check:
     1412            if k != self._ncols:
     1413                raise ValueError("self must be a square matrix")
     1414            if not self.is_alternating():
     1415                raise ValueError("self must be alternating, which includes the diagonal entries being 0")
     1416
     1417        R = self.base_ring()
     1418
     1419        if k % 2 == 1:
     1420            return R.zero()
     1421
     1422        if k == 0:
     1423            return R.one()
     1424
     1425        n = k // 2
     1426
     1427        res = R.zero()
     1428
     1429        from sage.combinat.perfect_matching import PerfectMatchings
     1430        from sage.misc.flatten import flatten
     1431        from sage.misc.misc_c import prod
     1432        from sage.combinat.permutation import Permutation
     1433        for m in PerfectMatchings(k):
     1434            # We only need each edge of the matching to be sorted,
     1435            # not the totality of edges.
     1436            edges = [sorted(edge) for edge in list(m)]
     1437            sgn = Permutation(flatten(edges)).sign()
     1438            # Subtract 1 from everything for indexing purposes:
     1439            edges2 = [[i-1 for i in edge] for edge in edges]
     1440            # Product without base case since k == 0 case has
     1441            # already been dealt with:
     1442            res += sgn * prod([self.get_unsafe(edge[0], edge[1]) for edge in edges2])
     1443
     1444        return res
    12451445
    12461446    # shortcuts
    12471447    def det(self, *args, **kwds):
    cdef class Matrix(matrix1.Matrix): 
    13221522            sage: factor(A.minpoly('y'))
    13231523            (y + 1) * (y + 2)^2
    13241524
    1325         We can take the minimal polynomail of symbolic matrices::
     1525        We can take the minimal polynomial of symbolic matrices::
    13261526
    13271527            sage: t = var('t')
    13281528            sage: m = matrix(2,[1,2,4,t])
    cdef class Matrix(matrix1.Matrix): 
    1297513175          be returned, where each list contains coefficients of a
    1297613176          polynomial associated with a companion matrix.
    1297713177
    12978         - `subdivide` - default: 'True' - if 'True' and a matrix is
     13178        - ``subdivide`` - default: 'True' - if 'True' and a matrix is
    1297913179          returned, then it contains subdivisions delineating the
    1298013180          companion matrices along the diagonal.
    1298113181
  • sage/misc/sagedoc.py

    diff --git a/sage/misc/sagedoc.py b/sage/misc/sagedoc.py
    a b def format(s, embedded=False): 
    408408
    409409        sage: from sage.misc.sagedoc import format
    410410        sage: identity_matrix(2).rook_vector.__doc__[115:184]
    411         'Let `A` be a general `m` by `n`\n        (0,1)-matrix with `m \\le n`. '
     411        '   Let `A` be an `m` by `n` (0,1)-matrix with `m \\le n`. We identify\n'
    412412        sage: format(identity_matrix(2).rook_vector.__doc__[115:184])
    413         'Let A be a general m by n\n   (0,1)-matrix with m <= n.\n'
     413        '   Let A be an m by n (0,1)-matrix with m <= n. We identify\n'
    414414
    415415    If the first line of the string is 'nodetex', remove 'nodetex' but
    416416    don't modify any TeX commands::