def invariant_form(self):
        """
        Return the quadratic form preserved by the symplectic group.

        OUTPUT:

        A matrix.

        EXAMPLES::

            sage: Sp(4, QQ).invariant_form()
            [ 0  0  0  1]
            [ 0  0  1  0]
            [ 0 -1  0  0]
            [-1  0  0  0]
        """
        if self._invariant_form is not None:
            return self._invariant_form

        R = self.base_ring()
        d = self.degree()
        from sage.matrix.constructor import zero_matrix
        m = zero_matrix(R, d)
        for i in range(d):
            m[i, d-i-1] = 1 if i < d/2 else -1
        m.set_immutable()
        return m

 def iter_positive_forms(self) :
     if self.is_reduced() :
         sub4 = self._calc_iter_reduced_sub4()
         
         for (a0, a1, b01, sub4s) in sub4 :
             t = zero_matrix(ZZ, 4)
             t[0,0] = a0
             t[1,1] = a1
             t[0,1] = b01
             t[1,0] = b01
             
             for (a2, b02, b12, sub4ss) in sub4s :
                 ts = copy(t)
                 ts[2,2] = a2
                 ts[0,2] = b02
                 ts[2,0] = b02
                 ts[1,2] = b12
                 ts[2,1] = b12
                 
                 for (a3, b03, b13, b23) in sub4ss :
                     tss = copy(ts)
                     tss[3,3] = a3
                     tss[0,1] = b01
                     tss[1,0] = b01
                     tss[0,2] = b02
                     tss[2,0] = b02
                     tss[0,3] = b03
                     tss[3,0] = b03
                     tss.set_immutable()
                     
                     yield tss
     else :
         raise NotImplementedError

    def in_degree(self, n):
        """
        The matrix representing this morphism in degree n

        INPUT:

        - ``n`` -- degree

        EXAMPLES::

            sage: C = ChainComplex({0: identity_matrix(ZZ, 1)})
            sage: D = ChainComplex({0: zero_matrix(ZZ, 1), 1: zero_matrix(ZZ, 1)})
            sage: f = Hom(C,D)({0: identity_matrix(ZZ, 1), 1: zero_matrix(ZZ, 1)})
            sage: f.in_degree(0)
            [1]

        Note that if the matrix is not specified in the definition of
        the map, it is assumed to be zero::

            sage: f.in_degree(2)
            []
            sage: f.in_degree(2).nrows(), f.in_degree(2).ncols()
            (1, 0)
            sage: C.free_module(2)
            Ambient free module of rank 0 over the principal ideal domain Integer Ring
            sage: D.free_module(2)
            Ambient free module of rank 1 over the principal ideal domain Integer Ring
        """
        try:
            return self._matrix_dictionary[n]
        except KeyError:
            rows = self.codomain().free_module_rank(n)
            cols = self.domain().free_module_rank(n)
            return zero_matrix(self.domain().base_ring(), rows, cols)

    def __apply_atkin_lehner(self,q,f):
        r"""
        This function applies an Atkin-Lehner involution to a harmonic cocycle

        INPUT:

          - ``q`` - an integer dividing the full level p*Nminus*Nplus

          - ``f`` - a harmonic cocycle

        OUTPUT:

          The harmonic cocycle obtained by hitting f with the Atkin-Lehner at q

        EXAMPLES:
        ::
        """
        R=self._R
        Data=self._X._get_atkin_lehner_data(q)
        p=self._X._p
        tmp=[self._U.element_class(self._U,zero_matrix(self._R,self._k-1,1),quick=True) for jj in range(len(self._E))]
        d1=Data[1]
        mga=self.embed_quaternion(Data[0])
        for jj in range(len(self._E)):
            t=d1[jj]
            tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.embed_quaternion))
        return HarmonicCocycleElement(self,tmp,from_values=True)

    def in_degree(self, n):
        """
        The matrix representing this chain homotopy in degree ``n``.

        INPUT:

        - ``n`` -- degree

        EXAMPLES::

            sage: from sage.homology.chain_homotopy import ChainHomotopy
            sage: C = ChainComplex({1: matrix(ZZ, 0, 2)}) # one nonzero term in degree 1
            sage: D = ChainComplex({0: matrix(ZZ, 0, 1)}) # one nonzero term in degree 0
            sage: f = Hom(C, D)({})
            sage: H = ChainHomotopy({1: matrix(ZZ, 1, 2, (3,1))}, f, f)
            sage: H.in_degree(1)
            [3 1]

        This returns an appropriately sized zero matrix if the chain
        homotopy is not defined in degree n::

            sage: H.in_degree(-3)
            []
        """
        try:
            return self._matrix_dictionary[n]
        except KeyError:
            from sage.matrix.constructor import zero_matrix
            deg = self.domain().degree_of_differential()
            rows = self.codomain().free_module_rank(n-deg)
            cols = self.domain().free_module_rank(n)
            return zero_matrix(self.domain().base_ring(), rows, cols)

    def __apply_hecke_operator(self,l,f):
        r"""
        This function applies a Hecke operator to a harmonic cocycle.

        INPUT:

          - ``l`` - an integer

          - ``f`` - a harmonic cocycle

        OUTPUT:

          A harmonic cocycle which is the result of applying the lth Hecke operator
          to f

        EXAMPLES:
        ::

        """
        R=self._R
        HeckeData,alpha=self._X._get_hecke_data(l)
        if(self.level()%l==0):
            factor=QQ(l**(Integer((self._k-2)/2))/(l+1))
        else:
            factor=QQ(l**(Integer((self._k-2)/2)))
        p=self._X._p
        alphamat=self.embed_quaternion(alpha)
        tmp=[self._U.element_class(self._U,zero_matrix(self._R,self._k-1,1),quick=True) for jj in range(len(self._E))]
        for ii in range(len(HeckeData)):
            d1=HeckeData[ii][1]
            mga=self.embed_quaternion(HeckeData[ii][0])*alphamat
            for jj in range(len(self._E)):
                t=d1[jj]
                tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.embed_quaternion))
        return HarmonicCocycleElement(self,[factor*x for x in tmp],from_values=True)

    def __compute_operator_matrix(self,T):
        r"""
        Compute the matrix of the operator ``T``.

        EXAMPLES::

        """
        R = self._R
        A = self.basis_matrix().transpose()
        basis = self.basis()
        B = zero_matrix(R,len(self._E) * (self._k-1),self.dimension())
        for rr in range(len(basis)):
            g = T(basis[rr])
            B.set_block(0,rr,Matrix(R,len(self._E) * (self._k-1),1,[g._F[e]._val[ii,0]  for e in range(len(self._E)) for ii in range(self._k-1) ]))

        try:
            if not R.is_exact():
                smin = min([a.valuation() for a in A.list()+B.list()])
                A = R.prime()**(-smin)*A
                B = R.prime()**(-smin)*B
                prec = min([a.precision_absolute() for a in A.list()+B.list()])
                A = A.parent()([R(x,absprec = prec) for x in A.list()])
                B = B.parent()([R(x,absprec = prec) for x in B.list()])
            res = (A.solve_right(B)).transpose()
            res.set_immutable()
        except ValueError:
            # save([A,B],'error.log.sobj')
            # print A.nrows(),A.ncols()
            # print B.nrows(),B.ncols()
            raise ValueError,'The hecke operator action is wrong.'
        return res

    def _rank(self, K) :
        
        if K is QQ or K in NumberFields() :
            return len(_jacobi_forms_by_taylor_expansion_coords(self.__index, self.__weight, 0))

            ## This is the formula used by Poor and Yuen in Paramodular cusp forms
            if self.__weight == 2 :
                delta = len(self.__index.divisors()) // 2 - 1
            else :
                delta = 0
                
            return sum( ModularForms(1, self.__weight + 2 * j).dimension() + j**2 // (4 * self.__index)
                        for j in xrange(self.__index + 1) ) \
                   + delta
            

            ## This is the formula given by Skoruppa in 
            ## Jacobi forms of critical weight and Weil representations
            ##FIXME: There is some mistake here
            if self.__weight % 2 != 0 :
                ## Otherwise the space X(i**(n - 2 k)) is different
                ## See: Skoruppa, Jacobi forms of critical weight and Weil representations
                raise NotImplementedError
            
            m = self.__index
            K = CyclotomicField(24 * m, 'zeta')
            zeta = K.gen(0)
            
            quadform = lambda x : 6 * x**2
            bilinform = lambda x,y : quadform(x + y) - quadform(x) - quadform(y)
            
            T = diagonal_matrix([zeta**quadform(i) for i in xrange(2*m)])
            S =   sum(zeta**(-quadform(x)) for x in xrange(2 * m)) / (2 * m) \
                * matrix([[zeta**(-bilinform(j,i)) for j in xrange(2*m)] for i in xrange(2*m)])
            subspace_matrix_1 = matrix( [ [1 if j == i or j == 2*m - i else 0 for j in xrange(m + 1) ]
                                        for i in xrange(2*m)] )
            subspace_matrix_2 = zero_matrix(ZZ, m + 1, 2*m)
            subspace_matrix_2.set_block(0,0,identity_matrix(m+1))
            
            T = subspace_matrix_2 * T * subspace_matrix_1
            S = subspace_matrix_2 * S * subspace_matrix_1
            
            sqrt3 = (zeta**(4*m) - zeta**(-4*m)) * zeta**(-6*m) 
            rank =   (self.__weight - 1/2 - 1) / 2 * (m + 1) \
                   + 1/8 * (   zeta**(3*m * (2*self.__weight - 1)) * S.trace()
                             + zeta**(3*m * (1 - 2*self.__weight)) * S.trace().conjugate() ) \
                   + 2/(3*sqrt3) * (   zeta**(4 * m * self.__weight) * (S*T).trace()
                                     + zeta**(-4 * m * self.__weight) * (S*T).trace().conjugate() ) \
                   - sum((j**2 % (m+1))/(m+1) -1/2 for j in range(0,m+1))
            
            if self.__weight > 5 / 2 :
                return rank
            else :
                raise NotImplementedError
            
        raise NotImplementedError

 def gen(self, i = 0) :
     if i < self.__n :
         t = diagonal_matrix(ZZ, i * [0] + [2] + (self.__n - i - 1) * [0])
         t.set_immutable()
         
         return t
     elif i >= self.__n and i < (self.__n * (self.__n + 1)) // 2 :
         i = i - self.__n
         
         for r in xrange(self.__n) :
             if i >=  self.__n - r - 1 :
                 i = i - (self.__n - r - 1)
                 continue
             
             c = i + r + 1
             break
         
         t = zero_matrix(ZZ, self.__n)
         t[r,c] = 1
         t[c,r] = 1
         
         t.set_immutable()
         return t
     elif not self.__reduced and i >= (self.__n * (self.__n + 1)) // 2 \
          and i < self.__n**2 :
         i = i - (self.__n * (self.__n + 1)) // 2
         
         for r in xrange(self.__n) :
             if i >=  self.__n - r - 1 :
                 i = i - (self.__n - r - 1)
                 continue
             
             c = i + r + 1
             break
         
         t = zero_matrix(ZZ, self.__n)
         t[r,c] = -1
         t[c,r] = -1
         
         t.set_immutable()
         return t
         
     raise ValueError, "Generator not defined"

    def matrix(self):
        """
        Return the standard matrix representation of ``self``.

        .. SEEALSO::

            - :meth:`AffineGroup.linear_space()`

        EXAMPLES::

            sage: G = AffineGroup(3, GF(7))
            sage: g = G([1,2,3,4,5,6,7,8,0], [10,11,12])
            sage: g
                  [1 2 3]     [3]
            x |-> [4 5 6] x + [4]
                  [0 1 0]     [5]
            sage: g.matrix()
            [1 2 3|3]
            [4 5 6|4]
            [0 1 0|5]
            [-----+-]
            [0 0 0|1]
            sage: parent(g.matrix())
            Full MatrixSpace of 4 by 4 dense matrices over Finite Field of size 7
            sage: g.matrix() == matrix(g)
            True

        Composition of affine group elements equals multiplication of
        the matrices::

            sage: g1 = G.random_element()
            sage: g2 = G.random_element()
            sage: g1.matrix() * g2.matrix() == (g1*g2).matrix()
            True
        """
        A = self._A
        b = self._b
        parent = self.parent()
        d = parent.degree()
        from sage.matrix.constructor import matrix, zero_matrix, block_matrix

        zero = zero_matrix(parent.base_ring(), 1, d)
        one = matrix(parent.base_ring(), [[1]])
        m = block_matrix(2, 2, [A, b.column(), zero, one])
        m.set_immutable()
        return m

    def _an_element_(self):
        """
        Construct a sample morphism.

        OUTPUT:

        An element of the homset.

        EXAMPLES::

            sage: P2 = toric_varieties.P2()
            sage: homset = P2.Hom(P2)
            sage: homset.an_element()   # indirect doctest
            Scheme endomorphism of 2-d CPR-Fano toric variety covered by 3 affine patches
              Defn: Defined by sending Rational polyhedral fan in 2-d lattice N to
                    Rational polyhedral fan in 2-d lattice N.
        """
        from sage.matrix.constructor import zero_matrix
        zero = zero_matrix(self.domain().dimension_relative(),
                           self.codomain().dimension_relative())
        return self(zero)

    def invariant_form(self):
        """
        Return the quadratic form preserved by the orthogonal group.

        OUTPUT:

        A matrix.

        EXAMPLES::

            sage: Sp(4, QQ).invariant_form()
            [0 0 0 1]
            [0 0 1 0]
            [0 1 0 0]
            [1 0 0 0]
        """
        from sage.matrix.constructor import zero_matrix
        m = zero_matrix(self.base_ring(), self.degree())
        for i in range(self.degree()):
            m[i, self.degree()-i-1] = 1
        m.set_immutable()
        return m

    def __compute_operator_matrix(self,T):
        r"""
        Compute the matrix of the operator ``T``.

        EXAMPLES:
        ::

        """
        R=self._R
        A=self.basis_matrix().transpose()
        basis=self.basis()
        B=zero_matrix(R,len(self._E)*(self._k-1),self.dimension())
        for rr in range(len(basis)):
            g=T(basis[rr])
            B.set_block(0,rr,Matrix(R,len(self._E)*(self._k-1),1,[g._F[e]._val[ii,0]  for e in range(len(self._E)) for ii in range(self._k-1) ]))

        try:
            res=(A.solve_right(B)).transpose()
            res.set_immutable()
            return res
        except ValueError:
            print A
            print B
            raise ValueError

#.........这里部分代码省略.........
                    1 |--> 1
                    2 |--> 2
            sage: a = x.associated_chain_complex_morphism()
            sage: a
            Chain complex morphism:
              From: Chain complex with at most 2 nonzero terms over Integer Ring
              To:   Chain complex with at most 3 nonzero terms over Integer Ring
            sage: a._matrix_dictionary
            {0: [1 0 0]
             [0 1 0]
             [0 0 1]
             [0 0 0], 1: [1 0 0]
             [0 1 0]
             [0 0 0]
             [0 0 1]
             [0 0 0]
             [0 0 0], 2: []}
            sage: x.associated_chain_complex_morphism(augmented=True)
            Chain complex morphism:
              From: Chain complex with at most 3 nonzero terms over Integer Ring
              To:   Chain complex with at most 4 nonzero terms over Integer Ring
            sage: x.associated_chain_complex_morphism(cochain=True)
            Chain complex morphism:
              From: Chain complex with at most 3 nonzero terms over Integer Ring
              To:   Chain complex with at most 2 nonzero terms over Integer Ring
            sage: x.associated_chain_complex_morphism(augmented=True,cochain=True)
            Chain complex morphism:
              From: Chain complex with at most 4 nonzero terms over Integer Ring
              To:   Chain complex with at most 3 nonzero terms over Integer Ring
            sage: x.associated_chain_complex_morphism(base_ring=GF(11))
            Chain complex morphism:
              From: Chain complex with at most 2 nonzero terms over Finite Field of size 11
              To:   Chain complex with at most 3 nonzero terms over Finite Field of size 11

        Some simplicial maps which reverse the orientation of a few simplices::

            sage: g = {0:1, 1:2, 2:0}
            sage: H(g).associated_chain_complex_morphism()._matrix_dictionary
            {0: [0 0 1]
             [1 0 0]
             [0 1 0]
             [0 0 0], 1: [ 0 -1  0]
             [ 0  0 -1]
             [ 0  0  0]
             [ 1  0  0]
             [ 0  0  0]
             [ 0  0  0], 2: []}
            sage: X = SimplicialComplex([[0, 1]], is_mutable=False)
            sage: Hom(X,X)({0:1, 1:0}).associated_chain_complex_morphism()._matrix_dictionary
            {0: [0 1]
             [1 0], 1: [-1]}
        """
        max_dim = max(self.domain().dimension(),self.codomain().dimension())
        min_dim = min(self.domain().dimension(),self.codomain().dimension())
        matrices = {}
        if augmented is True:
            m = matrix(base_ring,1,1,1)
            if not cochain:
                matrices[-1] = m
            else:
                matrices[-1] = m.transpose()
        for dim in range(min_dim+1):
            X_faces = list(self.domain().n_cells(dim))
            Y_faces = list(self.codomain().n_cells(dim))
            num_faces_X = len(X_faces)
            num_faces_Y = len(Y_faces)
            mval = [0 for i in range(num_faces_X*num_faces_Y)]
            for i in X_faces:
                y, oriented = self(i, orientation=True)
                if y.dimension() < dim:
                    pass
                else:
                    mval[X_faces.index(i)+(Y_faces.index(y)*num_faces_X)] = oriented
            m = matrix(base_ring,num_faces_Y,num_faces_X,mval,sparse=True)
            if not cochain:
                matrices[dim] = m
            else:
                matrices[dim] = m.transpose()
        for dim in range(min_dim+1,max_dim+1):
            try:
                l1 = len(self.codomain().n_cells(dim))
            except KeyError:
                l1 = 0
            try:
                l2 = len(self.domain().n_cells(dim))
            except KeyError:
                l2 = 0
            m = zero_matrix(base_ring,l1,l2,sparse=True)
            if not cochain:
                matrices[dim] = m
            else:
                matrices[dim] = m.transpose()
        if not cochain:
            return ChainComplexMorphism(matrices,\
                    self.domain().chain_complex(base_ring=base_ring,augmented=augmented,cochain=cochain),\
                    self.codomain().chain_complex(base_ring=base_ring,augmented=augmented,cochain=cochain))
        else:
            return ChainComplexMorphism(matrices,\
                    self.codomain().chain_complex(base_ring=base_ring,augmented=augmented,cochain=cochain),\
                    self.domain().chain_complex(base_ring=base_ring,augmented=augmented,cochain=cochain))

#.........这里部分代码省略.........
                    pass

            # One small typo in the published algorithm: it says
            # "if bdry(c_bar) == 0", but should say
            # "if pi(bdry(c_bar)) == 0".
            if not pi_bdry_c_bar:
                # Append c to list of gens.
                gens[dim].append(c)
                # iota(c) = c_bar
                iota_dict[dim][c] = c_bar
                # pi(c) = c
                pi_dict[dim][c_idx] = c_vec
            else:
                # Take any u in gens so that lambda_i = <u, pi(bdry(c_bar))> != 0.
                # u_idx will be the index of the corresponding cell.
                for (u_idx, lambda_i) in pi_bdry_c_bar.iteritems():
                    # Now find the actual cell.
                    u = old_cells[u_idx]
                    if u in gens[dim-1]:
                        break

                # pi(c) = 0: no need to do anything about this.
                for c_j_idx in range(old_rank):
                    # eta_ij = <u, pi(c_j)>.
                    try:
                        eta_ij = pi_dict[dim-1][c_j_idx][u_idx]
                    except (KeyError, IndexError):
                        eta_ij = 0
                    if eta_ij:
                        # Adjust phi(c_j).
                        try:
                            phi_dict[dim-1][c_j_idx] += eta_ij * lambda_i**(-1) * c_bar
                        except KeyError:
                            phi_dict[dim-1][c_j_idx] = eta_ij * lambda_i**(-1) * c_bar
                        # Adjust pi(c_j).
                        try:
                            pi_dict[dim-1][c_j_idx] += -eta_ij * lambda_i**(-1) * pi_bdry_c_bar
                        except KeyError:
                            pi_dict[dim-1][c_j_idx] = -eta_ij * lambda_i**(-1) * pi_bdry_c_bar

                gens[dim-1].remove(u)
                del iota_dict[dim-1][u]
        old_cells = n_cells

    # Now we have constructed the raw data for M, pi, iota, phi, so we
    # have to convert that to data which can be used to construct chain
    # complexes, chain maps, and chain contractions.

    # M_data will contain (trivial) matrices defining the differential
    # on M. Keep track of the sizes using "M_rows" and "M_cols", which are
    # just the ranks of consecutive graded pieces of M.
    M_data = {}
    M_rows = 0
    # pi_data: the matrices defining pi. Similar for iota_data and phi_data.
    pi_data = {}
    iota_data = {}
    phi_data = {}
    for n in range(K.dimension()+1):
        n_cells = K.n_cells(n)
        # Remove zero entries from pi_dict and phi_dict.
        pi_dict[n] = {i: pi_dict[n][i] for i in pi_dict[n] if pi_dict[n][i]}
        phi_dict[n] = {i: phi_dict[n][i] for i in phi_dict[n] if phi_dict[n][i]}
        # Convert gens to data defining the chain complex M with
        # trivial differential.
        M_cols = len(gens[n])
        M_data[n] = zero_matrix(base_ring, M_rows, M_cols)
        M_rows = M_cols
        # Convert the dictionaries for pi, iota, phi to matrices which
        # will define chain maps and chain homotopies.
        pi_cols = []
        phi_cols = []
        for (idx, c) in enumerate(n_cells):
            # First pi:
            if idx in pi_dict[n]:
                column = vector(base_ring, M_rows)
                for (entry, coeff) in pi_dict[n][idx].iteritems():
                    # Translate from cells in n_cells to cells in gens[n].
                    column[gens[n].index(n_cells[entry])] = coeff
            else:
                column = vector(base_ring, M_rows)
            pi_cols.append(column)

            # Now phi:
            try:
                column = phi_dict[n][idx]
            except KeyError:
                column = vector(base_ring, len(K.n_cells(n+1)))
            phi_cols.append(column)
        # Now iota:
        iota_cols = [iota_dict[n][c] for c in gens[n]]

        pi_data[n] = matrix(base_ring, pi_cols).transpose()
        iota_data[n] = matrix(base_ring, len(gens[n]), len(n_cells), iota_cols).transpose()
        phi_data[n] = matrix(base_ring, phi_cols).transpose()

    M = ChainComplex(M_data, base_ring=base_ring, degree=-1)
    pi = ChainComplexMorphism(pi_data, C, M)
    iota = ChainComplexMorphism(iota_data, M, C)
    phi = ChainContraction(phi_data, pi, iota)
    return phi, M

    def _find_isomorphism_degenerate(self, polytope):
        """
        Helper to pick an isomorphism of degenerate polygons

        INPUT:

        - ``polytope`` -- a :class:`LatticePolytope_PPL_class`. The
          polytope to compare with.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL, C_Polyhedron
            sage: L1 = LatticePolytope_PPL(C_Polyhedron(2, 'empty'))
            sage: L2 = LatticePolytope_PPL(C_Polyhedron(3, 'empty'))
            sage: iso = L1.find_isomorphism(L2)   # indirect doctest
            sage: iso(L1) == L2
            True
            sage: iso = L1._find_isomorphism_degenerate(L2)
            sage: iso(L1) == L2
            True

            sage: L1 = LatticePolytope_PPL((-1,4))
            sage: L2 = LatticePolytope_PPL((2,1,5))
            sage: iso = L1.find_isomorphism(L2)
            sage: iso(L1) == L2
            True

            sage: L1 = LatticePolytope_PPL((-1,), (3,))
            sage: L2 = LatticePolytope_PPL((2,1,5), (2,-3,5))
            sage: iso = L1.find_isomorphism(L2)
            sage: iso(L1) == L2
            True

            sage: L1 = LatticePolytope_PPL((-1,-1), (3,-1))
            sage: L2 = LatticePolytope_PPL((2,1,5), (2,-3,5))
            sage: iso = L1.find_isomorphism(L2)
            sage: iso(L1) == L2
            True

            sage: L1 = LatticePolytope_PPL((-1,2), (3,1))
            sage: L2 = LatticePolytope_PPL((1,2,3),(1,2,4))
            sage: iso = L1.find_isomorphism(L2)
            sage: iso(L1) == L2
            True

            sage: L1 = LatticePolytope_PPL((-1,2), (3,2))
            sage: L2 = LatticePolytope_PPL((1,2,3),(1,2,4))
            sage: L1.find_isomorphism(L2)
            Traceback (most recent call last):
            ...
            LatticePolytopesNotIsomorphicError: different number of integral points

            sage: L1 = LatticePolytope_PPL((-1,2), (3,1))
            sage: L2 = LatticePolytope_PPL((1,2,3),(1,2,5))
            sage: L1.find_isomorphism(L2)
            Traceback (most recent call last):
            ...
            LatticePolytopesNotIsomorphicError: different number of integral points
        """
        from sage.geometry.polyhedron.lattice_euclidean_group_element import \
            LatticePolytopesNotIsomorphicError
        polytope_vertices = polytope.vertices()
        self_vertices = self.ordered_vertices()
        # handle degenerate cases
        if self.n_vertices() == 0:
            A = zero_matrix(ZZ, polytope.space_dimension(), self.space_dimension())
            b = zero_vector(ZZ, polytope.space_dimension())
            return LatticeEuclideanGroupElement(A, b)
        if self.n_vertices() == 1:
            A = zero_matrix(ZZ, polytope.space_dimension(), self.space_dimension())
            b = polytope_vertices[0]
            return LatticeEuclideanGroupElement(A, b)
        if self.n_vertices() == 2:
            self_origin = self_vertices[0]
            self_ray = self_vertices[1] - self_origin
            polytope_origin = polytope_vertices[0]
            polytope_ray = polytope_vertices[1] - polytope_origin
            Ds, Us, Vs = self_ray.column().smith_form()
            Dp, Up, Vp = polytope_ray.column().smith_form()
            assert Vs.nrows() == Vs.ncols() == Vp.nrows() == Vp.ncols() == 1
            assert abs(Vs[0, 0]) == abs(Vp[0, 0]) == 1
            A = zero_matrix(ZZ, Dp.nrows(), Ds.nrows())
            A[0, 0] = 1
            A = Up.inverse() * A * Us * (Vs[0, 0] * Vp[0, 0])
            b = polytope_origin - A*self_origin
            try:
                A = matrix(ZZ, A)
                b = vector(ZZ, b)
            except TypeError:
                raise LatticePolytopesNotIsomorphicError('different lattice')
            hom = LatticeEuclideanGroupElement(A, b)
            if hom(self) == polytope:
                return hom
            raise LatticePolytopesNotIsomorphicError('different polygons')

#.........这里部分代码省略.........

        zero_vector = vector(base_ring, rank)
        pi_nrows = pi_old.nrows()

        for c_idx, c in enumerate(n_cells):
            # c_bar = c - phi(bdry(c)):
            # Avoid a bug in matrix-vector multiplication (trac 19378):
            if not diff:
                c_bar = c
                pi_bdry_c_bar = False
            else:
                if base_ring == QQ:
                    c_bar = c - phi_old * (diff * c)
                    pi_bdry_c_bar = conditionally_sparse(pi_old) * (diff * c_bar)
                else:
                    c_bar = c - phi_old * diff * c
                    pi_bdry_c_bar = conditionally_sparse(pi_old) * diff * c_bar

            # One small typo in the published algorithm: it says
            # "if bdry(c_bar) == 0", but should say
            # "if pi(bdry(c_bar)) == 0".
            if not pi_bdry_c_bar:
                # Append c to list of gens.
                gens[dim].append(c_idx)
                # iota(c) = c_bar
                iota_cols[c_idx] = c_bar
                # pi(c) = c
                pi_cols.append(c)
            else:
                # Take any u in gens so that lambda_i = <u, pi(bdry(c_bar))> != 0.
                # u_idx will be the index of the corresponding cell.
                (u_idx, lambda_i) = pi_bdry_c_bar.leading_item()
                for (u_idx, lambda_i) in pi_bdry_c_bar.iteritems():
                    if u_idx not in to_be_deleted:
                        break
                # This element/column needs to be deleted from gens and
                # iota_old. Do that later.
                to_be_deleted.append(u_idx)
                # pi(c) = 0.
                pi_cols.append(zero_vector)
                for c_j_idx, c_j in enumerate(old_cells):
                    # eta_ij = <u, pi(c_j)>.
                    # That is, eta_ij is the u_idx entry in the vector pi_old * c_j:
                    eta_ij = c_j.dot_product(pi_old.row(u_idx))
                    if eta_ij:
                        # Adjust phi(c_j).
                        phi_old_cols[c_j_idx] += eta_ij * lambda_i**(-1) * c_bar
                        # Adjust pi(c_j).
                        pi_cols_old[c_j_idx] -= eta_ij * lambda_i**(-1) * pi_bdry_c_bar

                # The matrices involved have many zero entries. For
                # such matrices, using sparse matrices is faster over
                # the rationals, slower over finite fields.
                phi_old = matrix(base_ring, phi_old_cols, sparse=(base_ring==QQ)).transpose()
                keep = vector(base_ring, pi_nrows, {i:1 for i in range(pi_nrows)
                                                    if i not in to_be_deleted})
                cols = [v.pairwise_product(keep) for v in pi_cols_old]
                pi_old = MS_pi_t.matrix(cols).transpose()

        # Here cols is a temporary storage for the columns of iota.
        cols = [iota_cols_old[i] for i in sorted(iota_cols_old.keys())]
        for r in sorted(to_be_deleted, reverse=True):
            del cols[r]
            del gens[dim-1][r]
        iota_data[dim-1] = matrix(base_ring, len(gens[dim-1]), old_rank, cols).transpose()
        # keep: rows to keep in pi_cols_old. Start with all
        # columns, then delete those in to_be_deleted.
        keep = sorted(set(range(pi_nrows)).difference(to_be_deleted))
        # Now cols is a temporary storage for columns of pi.
        cols = [v.list_from_positions(keep) for v in pi_cols_old]
        pi_data[dim-1] = matrix(base_ring, old_rank, len(gens[dim-1]), cols).transpose()
        phi_data[dim-1] = phi_old

        V_gens = VectorSpace(base_ring, len(gens[dim]))
        if pi_cols:
            cols = []
            for v in pi_cols:
                cols.append(V_gens(v.list_from_positions(gens[dim])))
            pi_cols = cols

    pi_data[dim] = matrix(base_ring, rank, len(gens[dim]), pi_cols).transpose()
    cols = [iota_cols[i] for i in sorted(iota_cols.keys())]
    iota_data[dim] = matrix(base_ring, len(gens[dim]), rank, cols).transpose()

    # M_data will contain (trivial) matrices defining the differential
    # on M. Keep track of the sizes using "M_rows" and "M_cols", which are
    # just the ranks of consecutive graded pieces of M.
    M_data = {}
    M_rows = 0
    for n in range(K.dimension()+1):
        M_cols = len(gens[n])
        M_data[n] = zero_matrix(base_ring, M_rows, M_cols)
        M_rows = M_cols

    M = ChainComplex(M_data, base_ring=base_ring, degree=-1)

    pi = ChainComplexMorphism(pi_data, C, M)
    iota = ChainComplexMorphism(iota_data, M, C)
    phi = ChainContraction(phi_data, pi, iota)
    return phi, M

def SymplecticPolarGraph(d, q, algorithm=None):
    r"""
    Returns the Symplectic Polar Graph `Sp(d,q)`.

    The Symplectic Polar Graph `Sp(d,q)` is built from a projective space of dimension
    `d-1` over a field `F_q`, and a symplectic form `f`. Two vertices `u,v` are
    made adjacent if `f(u,v)=0`.

    See the page `on symplectic graphs on Andries Brouwer's website
    <http://www.win.tue.nl/~aeb/graphs/Sp.html>`_.

    INPUT:

    - ``d,q`` (integers) -- note that only even values of `d` are accepted by
      the function.

    - ``algorithm`` -- if set to 'gap' then the computation is carried via GAP
      library interface, computing totally singular subspaces, which is faster for `q>3`.
      Otherwise it is done directly.

    EXAMPLES:

    Computation of the spectrum of `Sp(6,2)`::

        sage: g = graphs.SymplecticGraph(6,2)
        doctest:...: DeprecationWarning: SymplecticGraph is deprecated. Please use sage.graphs.generators.classical_geometries.SymplecticPolarGraph instead.
        See http://trac.sagemath.org/19136 for details.
        sage: g.is_strongly_regular(parameters=True)
        (63, 30, 13, 15)
        sage: set(g.spectrum()) == {-5, 3, 30}
        True

    The parameters of `Sp(4,q)` are the same as of `O(5,q)`, but they are
    not isomorphic if `q` is odd::

        sage: G = graphs.SymplecticPolarGraph(4,3)
        sage: G.is_strongly_regular(parameters=True)
        (40, 12, 2, 4)
        sage: O=graphs.OrthogonalPolarGraph(5,3)
        sage: O.is_strongly_regular(parameters=True)
        (40, 12, 2, 4)
        sage: O.is_isomorphic(G)
        False
        sage: graphs.SymplecticPolarGraph(6,4,algorithm="gap").is_strongly_regular(parameters=True) # not tested (long time)
        (1365, 340, 83, 85)

    TESTS::

        sage: graphs.SymplecticPolarGraph(4,4,algorithm="gap").is_strongly_regular(parameters=True)
        (85, 20, 3, 5)
        sage: graphs.SymplecticPolarGraph(4,4).is_strongly_regular(parameters=True)
        (85, 20, 3, 5)
        sage: graphs.SymplecticPolarGraph(4,4,algorithm="blah")
        Traceback (most recent call last):
        ...
        ValueError: unknown algorithm!
    """
    if d < 1 or d%2 != 0:
        raise ValueError("d must be even and greater than 2")

    if algorithm == "gap":     # faster for larger (q>3)  fields
        from sage.libs.gap.libgap import libgap
        G = _polar_graph(d, q, libgap.SymplecticGroup(d, q))

    elif algorithm == None:    # faster for small (q<4) fields
        from sage.modules.free_module import VectorSpace
        from sage.schemes.projective.projective_space import ProjectiveSpace
        from sage.matrix.constructor import identity_matrix, block_matrix, zero_matrix

        F = FiniteField(q,"x")
        M = block_matrix(F, 2, 2,
                         [zero_matrix(F,d/2),
                          identity_matrix(F,d/2),
                          -identity_matrix(F,d/2),
                          zero_matrix(F,d/2)])

        V = VectorSpace(F,d)
        PV = list(ProjectiveSpace(d-1,F))
        G = Graph([[tuple(_) for _ in PV], lambda x,y:V(x)*(M*V(y)) == 0], loops = False)

    else:
        raise ValueError("unknown algorithm!")

    G.name("Symplectic Polar Graph Sp("+str(d)+","+str(q)+")")
    G.relabel()
    return G

    def __iter__(self) :
        if self.index() is infinity :
            raise ValueError, "infinity is not a true filter index"

        if self.is_reduced() :
            ## We only iterate positive definite matrices
            ## and later build the semidefinite ones
            ## We first find possible upper left matrices
            
            sub2 = self._calc_iter_reduced_sub2()
            sub3 = self._calc_iter_reduced_sub3()
            sub4 = self._calc_iter_reduced_sub4()

                
            t = zero_matrix(ZZ, 4)
            t.set_immutable()
            yield t
            
            for a0 in xrange(2, 2 * self.index(), 2) :
                t = zero_matrix(ZZ, 4)
                t[3,3] = a0
                t.set_immutable()
                
                yield t
                
            for (a0, a1, b01) in sub2 :
                t = zero_matrix(ZZ, 4)
                t[2,2] = a0
                t[3,3] = a1
                t[2,3] = b01
                t[3,2] = b01
                t.set_immutable()
                
                yield t
                
            for (a0, a1, b01, sub3s) in sub3 :
                t = zero_matrix(ZZ, 4)
                t[1,1] = a0
                t[2,2] = a1
                t[1,2] = b01
                t[2,1] = b01
                
                for (a2, b02, b12) in sub3s :
                    ts = copy(t)
                    ts[3,3] = a2
                    ts[1,3] = b02
                    ts[3,1] = b02
                    ts[2,3] = b12
                    ts[3,2] = b12