def space(self):
        r'''
        Calculates the space of cocyles modulo coboundaries, as a Z-module.

        TESTS:

        sage: from darmonpoints.sarithgroup import *
        sage: from darmonpoints.cohomology_abstract import *
        sage: from darmonpoints.ocmodule import *
        sage: GS = BigArithGroup(5, 6,1,use_shapiro=False,outfile='/tmp/darmonpoints.tmp') #  optional - magma
        sage: G = GS.large_group() #  optional - magma
        sage: V = OCVn(5,1)     #  optional - magma
        sage: Coh = CohomologyGroup(G,V,trivial_action = False) #  optional - magma
        '''
        verb = get_verbose()
        set_verbose(0)
        V = self.coefficient_module()
        R = V.base_ring()
        Vdim = V.dimension()
        G = self.group()
        gens = G.gens()
        ambient = R**(Vdim * len(gens))
        # Now find the subspace of cocycles
        A = Matrix(R, Vdim * len(gens), 0)
        for r in G.get_relation_words():
            Alist = self.fox_gradient(r)
            newA = block_matrix(Alist, nrows = 1)
            A = A.augment(newA.transpose())
        A = A.transpose()
        cocycles = ambient.submodule([ambient(o) for o in A.right_kernel_matrix().rows()])
        gmat = block_matrix([self._gen_pows[i][1] - 1 for i in range(len(G.gens()))], nrows = len(G.gens()))
        coboundaries = cocycles.submodule([ambient(o) for o in gmat.columns()])
        ans = cocycles.quotient(coboundaries)
        set_verbose(verb)
        return ans

def Hadamard3Design(n):
    """
    Return the Hadamard 3-design with parameters `3-(n, \\frac n 2, \\frac n 4 - 1)`.

    This is the unique extension of the Hadamard `2`-design (see
    :meth:`HadamardDesign`).  We implement the description from pp. 12 in
    [CvL]_.

    INPUT:

    - ``n`` (integer) -- a multiple of 4 such that `n>4`.

    EXAMPLES::

        sage: designs.Hadamard3Design(12)
        Incidence structure with 12 points and 22 blocks

    We verify that any two blocks of the Hadamard `3`-design `3-(8, 4, 1)`
    design meet in `0` or `2` points. More generally, it is true that any two
    blocks of a Hadamard `3`-design meet in `0` or `\\frac{n}{4}` points (for `n
    > 4`).

    ::

        sage: D = designs.Hadamard3Design(8)
        sage: N = D.incidence_matrix()
        sage: N.transpose()*N
        [4 2 2 2 2 2 2 2 2 2 2 2 2 0]
        [2 4 2 2 2 2 2 2 2 2 2 2 0 2]
        [2 2 4 2 2 2 2 2 2 2 2 0 2 2]
        [2 2 2 4 2 2 2 2 2 2 0 2 2 2]
        [2 2 2 2 4 2 2 2 2 0 2 2 2 2]
        [2 2 2 2 2 4 2 2 0 2 2 2 2 2]
        [2 2 2 2 2 2 4 0 2 2 2 2 2 2]
        [2 2 2 2 2 2 0 4 2 2 2 2 2 2]
        [2 2 2 2 2 0 2 2 4 2 2 2 2 2]
        [2 2 2 2 0 2 2 2 2 4 2 2 2 2]
        [2 2 2 0 2 2 2 2 2 2 4 2 2 2]
        [2 2 0 2 2 2 2 2 2 2 2 4 2 2]
        [2 0 2 2 2 2 2 2 2 2 2 2 4 2]
        [0 2 2 2 2 2 2 2 2 2 2 2 2 4]


    REFERENCES:

    .. [CvL] P. Cameron, J. H. van Lint, Designs, graphs, codes and
      their links, London Math. Soc., 1991.
    """
    if n == 1 or n == 4:
        raise ValueError("The Hadamard design with n = %s does not extend to a three design." % n)
    from sage.combinat.matrices.hadamard_matrix import hadamard_matrix
    from sage.matrix.constructor import matrix, block_matrix
    H = hadamard_matrix(n) #assumed to be normalised.
    H1 = H.matrix_from_columns(range(1, n))
    J = matrix(ZZ, n, n-1, [1]*(n-1)*n)
    A1 = (H1+J)/2
    A2 = (J-H1)/2
    A = block_matrix(1, 2, [A1, A2]) #the incidence matrix of the design.
    return IncidenceStructure(incidence_matrix=A, name="HadamardThreeDesign")

def RSHCD_324(e):
    r"""
    Return a size 324x324 Regular Symmetric Hadamard Matrix with Constant Diagonal.

    We build the matrix `M` for the case `n=324`, `\epsilon=1` directly from
    :meth:`JankoKharaghaniTonchevGraph
    <sage.graphs.graph_generators.GraphGenerators.JankoKharaghaniTonchevGraph>`
    and for the case `\epsilon=-1` from the "twist" `M'` of `M`, using Lemma 11
    in [HX10]_. Namely, it turns out that the matrix

    .. MATH::

        M'=\begin{pmatrix} M_{12} & M_{11}\\ M_{11}^\top & M_{21} \end{pmatrix},
        \quad\text{where}\quad
        M=\begin{pmatrix} M_{11} & M_{12}\\ M_{21} & M_{22} \end{pmatrix},

    and the `M_{ij}` are 162x162-blocks, also RSHCD, its diagonal blocks having zero row
    sums, as needed by [loc.cit.]. Interestingly, the corresponding
    `(324,152,70,72)`-strongly regular graph
    has a vertex-transitive automorphism group of order 2592, twice the order of the
    (intransitive) automorphism group of the graph corresponding to `M`. Cf. [CP16]_.

    INPUT:

    - ``e`` -- one of `-1` or `+1`, equal to the value of `\epsilon`

    TESTS::

        sage: from sage.combinat.matrices.hadamard_matrix import RSHCD_324, is_hadamard_matrix
        sage: for e in [1,-1]:  # long time
        ....:     M = RSHCD_324(e)
        ....:     print("{} {} {}".format(M==M.T,is_hadamard_matrix(M),all([M[i,i]==1 for i in range(324)])))
        ....:     print(set(map(sum,M)))
        True True True
        set([18])
        True True True
        set([-18])

    REFERENCE:

    .. [CP16] \N. Cohen, D. Pasechnik,
       Implementing Brouwer's database of strongly regular graphs,
       Designs, Codes, and Cryptography, 2016
       :doi:`10.1007/s10623-016-0264-x`
    """

    from sage.graphs.generators.smallgraphs import JankoKharaghaniTonchevGraph as JKTG
    M = JKTG().adjacency_matrix()
    M = J(324) - 2*M
    if e==-1:
        M1=M[:162].T
        M2=M[162:].T
        M11=M1[:162]
        M12=M1[162:].T
        M21=M2[:162].T
        M=block_matrix([[M12,-M11],[-M11.T,M21]])
    return M

def williamson_goethals_seidel_skew_hadamard_matrix(a, b, c, d, check=True):
    r"""
    Williamson-Goethals-Seidel construction of a skew Hadamard matrix

    Given `n\times n` (anti)circulant matrices `A`, `B`, `C`, `D` with 1,-1 entries, 
    and satisfying `A+A^\top = 2I`, `AA^\top + BB^\top + CC^\top + DD^\top = 4nI`,
    one can construct a skew Hadamard matrix of order `4n`, cf. [GS70s]_.

    INPUT:

    - ``a`` -- 1,-1 list specifying the 1st row of `A`

    - ``b`` -- 1,-1 list specifying the 1st row of `B`

    - ``d`` -- 1,-1 list specifying the 1st row of `C`

    - ``c`` -- 1,-1 list specifying the 1st row of `D`

    EXAMPLES::

        sage: from sage.combinat.matrices.hadamard_matrix import williamson_goethals_seidel_skew_hadamard_matrix as WGS
        sage: a=[ 1,  1, 1, -1,  1, -1,  1, -1, -1]
        sage: b=[ 1, -1, 1,  1, -1, -1,  1,  1, -1]
        sage: c=[-1, -1]+[1]*6+[-1]
        sage: d=[ 1,  1, 1, -1,  1,  1, -1,  1,  1]
        sage: M=WGS(a,b,c,d,check=True)

    REFERENCES:

    .. [GS70s] \J.M. Goethals and J. J. Seidel,
      A skew Hadamard matrix of order 36,
      J. Aust. Math. Soc. 11(1970), 343-344
    .. [Wall71] \J. Wallis,
      A skew-Hadamard matrix of order 92,
      Bull. Aust. Math. Soc. 5(1971), 203-204
    .. [KoSt08] \C. Koukouvinos, S. Stylianou
      On skew-Hadamard matrices,
      Discrete Math. 308(2008) 2723-2731


    """
    n = len(a)
    R = matrix(ZZ, n, n, lambda i,j: 1 if i+j==n-1 else 0)
    A,B,C,D=map(matrix.circulant, [a,b,c,d])
    if check:
        assert A*A.T+B*B.T+C*C.T+D*D.T==4*n*I(n)
        assert A+A.T==2*I(n)

    M = block_matrix([[   A,    B*R,    C*R,    D*R],
                      [-B*R,      A, -D.T*R,  C.T*R],
                      [-C*R,  D.T*R,      A, -B.T*R],
                      [-D*R, -C.T*R,  B.T*R,      A]])
    if check:
        assert is_hadamard_matrix(M, normalized=False, skew=True)
    return M

    def cone_points_iter(self):
        """
        Iterate over the open torus orbits and yield distinct points.
        
        OUTPUT:

        For each open torus orbit (cone): A triple consisting of the
        cone, the nonzero homogeneous coordinates in that orbit (list
        of integers), and the nonzero log coordinates of distinct
        points as a cokernel.

        EXAMPLES::

            sage: fan = NormalFan(ReflexivePolytope(2, 0))
            sage: X = ToricVariety(fan, base_ring=GF(7))
            sage: point_set = X.point_set()
            sage: ffe = point_set._finite_field_enumerator()
            sage: cpi = ffe.cone_points_iter()
            sage: cone, nonzero_points, cokernel = list(cpi)[5]
            sage: cone
            1-d cone of Rational polyhedral fan in 2-d lattice N
            sage: cone.ambient_ray_indices()
            (2,)
            sage: nonzero_points
            [0, 1]
            sage: cokernel
            Finitely generated module V/W over Integer Ring with invariants (2)
            sage: list(cokernel)
            [(0), (1)]
            sage: [p.lift() for p in cokernel]
            [(0, 0), (0, 1)]
        """
        from sage.matrix.constructor import matrix, block_matrix, identity_matrix
        from sage.rings.all import ZZ
        nrays = len(self.rays())
        N = self.multiplicative_group_order()
        # Want cokernel of the log rescalings in (ZZ/N)^(#rays). But
        # ZZ/N is not a integral domain. Instead: work over ZZ
        log_generators = self.rescaling_log_generators()
        log_relations = block_matrix(2, 1, [
            matrix(ZZ, len(log_generators), nrays, log_generators), 
            N * identity_matrix(ZZ, nrays)])
        for cone in self.cone_iter():
            nrays = self.fan().nrays() + len(self.fan().virtual_rays())
            nonzero_coordinates = [i for i in range(nrays)
                                   if i not in cone.ambient_ray_indices()]
            log_relations_nonzero = log_relations.matrix_from_columns(nonzero_coordinates)
            image = log_relations_nonzero.image()
            cokernel = image.ambient_module().quotient(image)
            yield cone, nonzero_coordinates, cokernel

def hadamard_matrix_paleyII(n):
    """
    Implements the Paley type II construction.

    EXAMPLES::

        sage: sage.combinat.matrices.hadamard_matrix.hadamard_matrix_paleyI(12).det()
        2985984
        sage: 12^6
        2985984
    """
    N = Integer(n/2)
    p = N - 1
    if not(is_prime(p) and (p % 4 == 1)):
        raise ValueError("The order %s is not covered by the Paley type II construction." % n)
    S = matrix(ZZ, [[H2(i, j, p) for i in range(N)] for j in range(N)])
    return block_matrix([[S + 1, S - 1], [S - 1, -S - 1]])

    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 _repr_(self):
        """
        Return a string representation.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.double_description_inhomogeneous import Hrep2Vrep
            sage: H = Hrep2Vrep(QQ, 1, [(1,2)], [])
            sage: H._repr_()
            '[-1/2| 1/2|]'
        """
        from sage.matrix.constructor import block_matrix

        def make_matrix(rows):
             return matrix(self.base_ring, len(rows), self.dim, rows).transpose()
        V = make_matrix(self.vertices)
        R = make_matrix(self.rays)
        L = make_matrix(self.lines)
        return str(block_matrix([[V, R, L]]))

def matrice_systeme(systeme, variables):
    """
    Renvoie une matrice par block représentant un programme linéaire sous forme standard.

    INPUT::

    - ``systeme`` -- Un programme linéaire sous forme standard
    - ``variables`` -- La liste des variables du système

    EXAMPLES::

        sage: x = x1,x2,x3 = var('x1,x2,x3')
        sage: Chvatal13 = [[2*x1 + 3*x2 +   x3 <=  5,
        ....:               4*x1 +   x2 + 2*x3 <= 11,
        ....:               3*x1 + 4*x2 + 2*x3 <=  8],
        ....:               5*x1 + 4*x2 + 3*x3]

        sage: m = matrice_systeme(Chvatal13, x); m
        [ z|s1 s2 s3|x1 x2 x3| 0]
        [--+--------+--------+--]
        [ 1| 0  0  0|-5 -4 -3| 0]
        [--+--------+--------+--]
        [ 0| 1  0  0| 2  3  1| 5]
        [ 0| 0  1  0| 4  1  2|11]
        [ 0| 0  0  1| 3  4  2| 8]
    """

    def liste_coeffs(expression):
        return [expression.coeff(v) for v in variables]

    inequations = systeme[0]
    m = matrix([liste_coeffs(ineq.lhs()) for ineq in inequations])
    rhs = vector(ineq.rhs() for ineq in inequations).column()
    slack = SR.var(",".join("s%s" % i for i in range(1, len(inequations) + 1)))
    z = SR.var("z")
    return block_matrix(
        [
            [z, matrix([slack]), matrix([variables]), ZZ(0)],
            [ZZ(1), ZZ(0), -matrix([liste_coeffs(systeme[1])]), ZZ(0)],
            [ZZ(0), ZZ(1), m, rhs],
        ]
    )

def hadamard_matrix_paleyII(n):
    """
    Implements the Paley type II construction.

    The Paley type II case corresponds to the case `p \cong 1 \mod{4}` for a
    prime `p` (see [Hora]_).

    EXAMPLES::

        sage: sage.combinat.matrices.hadamard_matrix.hadamard_matrix_paleyII(12).det()
        2985984
        sage: 12^6
        2985984

    We note that the method returns a normalised Hadamard matrix ::

        sage: sage.combinat.matrices.hadamard_matrix.hadamard_matrix_paleyII(12)
        [ 1  1  1  1  1  1| 1  1  1  1  1  1]
        [ 1  1  1 -1 -1  1|-1 -1  1 -1 -1  1]
        [ 1  1  1  1 -1 -1|-1  1 -1  1 -1 -1]
        [ 1 -1  1  1  1 -1|-1 -1  1 -1  1 -1]
        [ 1 -1 -1  1  1  1|-1 -1 -1  1 -1  1]
        [ 1  1 -1 -1  1  1|-1  1 -1 -1  1 -1]
        [-----------------+-----------------]
        [ 1 -1 -1 -1 -1 -1|-1  1  1  1  1  1]
        [ 1 -1  1 -1 -1  1| 1 -1 -1  1  1 -1]
        [ 1  1 -1  1 -1 -1| 1 -1 -1 -1  1  1]
        [ 1 -1  1 -1  1 -1| 1  1 -1 -1 -1  1]
        [ 1 -1 -1  1 -1  1| 1  1  1 -1 -1 -1]
        [ 1  1 -1 -1  1 -1| 1 -1  1  1 -1 -1]
    """
    N = Integer(n/2)
    p = N - 1
    if not(is_prime(p) and (p % 4 == 1)):
        raise ValueError("The order %s is not covered by the Paley type II construction." % n)
    S = matrix(ZZ, [[H2(i, j, p) for i in range(N)] for j in range(N)])
    H = block_matrix([[S + 1, S - 1], [1 - S, S + 1]])
    # normalising H so that first row and column have only +1 entries.
    return normalise_hadamard(H)

#.........这里部分代码省略.........
        deg = [coefficients[0].degree(), coefficients[1].degree(),
                coefficients[2].degree()]
        # definitions as in [HC2006] and [ACKERMANS2016]
        A = ((deg[1] + deg[2]) / 2).ceil() - case
        B = ((deg[2] + deg[0]) / 2).ceil() - case
        C = ((deg[0] + deg[1]) / 2).ceil() - case
        
        # For all roots as calculated by has_rational_point(), we create
        # a system of linear equations. As in [ACKERMANS2016], we do this
        # by calculating the matrices for all phi_p, with basis consisting
        # of monomials of x, y and z in the space V of potential solutions:
        # t^0, ..., t^A, t^0, ..., t^B and t^0, ..., t^C.
        phi = []
        for (i, p) in enumerate(supports[0]):
            # lift to F[t] and map to R, with R as defined above
            if roots[0][i].parent().is_finite():
                root = roots[0][i].polynomial()
            else:
                root = roots[0][i].lift()
            alpha = root.parent().hom([t])(root)
            d = p.degree()
            # Calculate y - alpha*z mod p for all basis vectors
            phi_p = [[] for i in range(A+B+C+4)]
            phi_p[0:A+1] = [vector(F, d)] * (A+1)
            phi_p[A+1] = vector(F, d, {0: F(1)})
            lastpoly = F(1)
            for n in range(B):
                lastpoly = (lastpoly * t) % p
                phi_p[A+2+n] = vector(F, d, lastpoly.dict())
            lastpoly = -alpha % p
            phi_p[A+B+2] = vector(F, d, lastpoly.dict())
            for n in range(C):
                lastpoly = (lastpoly * t) % p
                phi_p[A+B+3+n] = vector(F, d, lastpoly.dict())
            phi_p[A+B+C+3] = vector(F, d)
            phi.append(matrix(phi_p).transpose())
        for (i, p) in enumerate(supports[1]):
            if roots[1][i].parent().is_finite():
                root = roots[1][i].polynomial()
            else:
                root = roots[1][i].lift()
            alpha = root.parent().hom([t])(root)
            d = p.degree()
            # Calculate z - alpha*x mod p for all basis vectors
            phi_p = [[] for i in range(A+B+C+4)]
            phi_p[A+1:A+B+2] = [vector(F, d)] * (B+1)
            phi_p[A+B+2] = vector(F, d, {0: F(1)})
            lastpoly = F(1)
            for n in range(C):
                lastpoly = (lastpoly * t) % p
                phi_p[A+B+3+n] = vector(F, d, lastpoly.dict())
            lastpoly = -alpha % p
            phi_p[0] = vector(F, d, lastpoly.dict())
            for n in range(A):
                lastpoly = (lastpoly * t) % p
                phi_p[1+n] = vector(F, d, lastpoly.dict())
            phi_p[A+B+C+3] = vector(F, d)
            phi.append(matrix(phi_p).transpose())
        for (i, p) in enumerate(supports[2]):
            if roots[2][i].parent().is_finite():
                root = roots[2][i].polynomial()
            else:
                root = roots[2][i].lift()
            alpha = root.parent().hom([t])(root)
            d = p.degree()
            # Calculate x - alpha*y mod p for all basis vectors
            phi_p = [[] for i in range(A+B+C+4)]
            phi_p[A+B+2:A+B+C+3] = [vector(F, d)] * (C+1)
            phi_p[0] = vector(F, d, {0: F(1)})
            lastpoly = F(1)
            for n in range(A):
                lastpoly = (lastpoly * t) % p
                phi_p[1+n] = vector(F, d, lastpoly.dict())
            lastpoly = -alpha % p
            phi_p[A+1] = vector(F, d, lastpoly.dict())
            for n in range(B):
                lastpoly = (lastpoly * t) % p
                phi_p[A+2+n] = vector(F, d, lastpoly.dict())
            phi_p[A+B+C+3] = vector(F, d)
            phi.append(matrix(phi_p).transpose())
        if case == 0:
            # We need three more equations
            lx = Ft(solution[0]).leading_coefficient()
            ly = Ft(solution[1]).leading_coefficient()
            lz = Ft(solution[2]).leading_coefficient()
            phi.append(matrix([vector(F, A+B+C+4, {A:1, A+B+C+3:-lx}),
                vector(F, A+B+C+4, {A+B+1:1, A+B+C+3:-ly}),
                vector(F, A+B+C+4, {A+B+C+2: 1, A+B+C+3:-lz})]))
        # Create the final matrix which we will solve
        M = block_matrix(phi, ncols = 1, subdivide = False)
        solution_space = M.right_kernel()
        for v in solution_space.basis():
            if v[:A+B+C+3] != 0: # we don't want to return a trivial solution
                X = Ft(list(v[:A+1]))
                Y = Ft(list(v[A+1:A+B+2]))
                Z = Ft(list(v[A+B+2:A+B+C+3]))
                return self.point([X,Y,Z])

        raise RuntimeError("No solution has been found: possibly incorrect\
 solubility certificate.")

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

#.........这里部分代码省略.........
        sage: sage.crypto.gen_lattice(m=10, q=11, seed=42, lattice=True)
        Free module of degree 10 and rank 10 over Integer Ring
        User basis matrix:
        [ 0  0  1  1  0 -1 -1 -1  1  0]
        [-1  1  0  1  0  1  1  0  1  1]
        [-1  0  0  0 -1  1  1 -2  0  0]
        [-1 -1  0  1  1  0  0  1  1 -1]
        [ 1  0 -1  0  0  0 -2 -2  0  0]
        [ 2 -1  0  0  1  0  1  0  0 -1]
        [-1  1 -1  0  1 -1  1  0 -1 -2]
        [ 0  0 -1  3  0  0  0 -1 -1 -1]
        [ 0 -1  0 -1  2  0 -1  0  0  2]
        [ 0  1  1  0  1  1 -2  1 -1 -2]

    REFERENCES:

    .. [A96] Miklos Ajtai.
      Generating hard instances of lattice problems (extended abstract).
      STOC, pp. 99--108, ACM, 1996.

    .. [GM02] Daniel Goldstein and Andrew Mayer.
      On the equidistribution of Hecke points.
      Forum Mathematicum, 15:2, pp. 165--189, De Gruyter, 2003.

    .. [LM06] Vadim Lyubashevsky and Daniele Micciancio.
      Generalized compact knapsacks are collision resistant.
      ICALP, pp. 144--155, Springer, 2006.

    .. [R05] Oded Regev.
      On lattices, learning with errors, random linear codes, and cryptography.
      STOC, pp. 84--93, ACM, 2005.
    """
    from sage.rings.finite_rings.integer_mod_ring import IntegerModRing
    from sage.matrix.constructor import identity_matrix, block_matrix
    from sage.matrix.matrix_space import MatrixSpace
    from sage.rings.integer_ring import IntegerRing
    if seed is not None:
        from sage.misc.randstate import set_random_seed
        set_random_seed(seed)

    if type == 'random':
        if n != 1: raise ValueError('random bases require n = 1')

    ZZ = IntegerRing()
    ZZ_q = IntegerModRing(q)
    A = identity_matrix(ZZ_q, n)

    if type == 'random' or type == 'modular':
        R = MatrixSpace(ZZ_q, m-n, n)
        A = A.stack(R.random_element())

    elif type == 'ideal':
        if quotient is None:
            raise ValueError('ideal bases require a quotient polynomial')
        try:
            quotient = quotient.change_ring(ZZ_q)
        except (AttributeError, TypeError):
            quotient = quotient.polynomial(base_ring=ZZ_q)

        P = quotient.parent()
        # P should be a univariate polynomial ring over ZZ_q
        if not is_PolynomialRing(P):
            raise TypeError("quotient should be a univariate polynomial")
        assert P.base_ring() is ZZ_q

        if quotient.degree() != n:

    def _find_cyclic_isomorphism_matching_edge(self, polytope,
                                               polytope_origin, p_ray_left,
                                               p_ray_right):
        """
        Helper to find an isomorphism of polygons

        INPUT:

        - ``polytope`` -- the lattice polytope to compare to.

        - ``polytope_origin`` -- `\ZZ`-vector. a vertex of ``polytope``

        - ``p_ray_left`` - vector. the vector from ``polytope_origin``
          to one of its neighboring vertices.

        - ``p_ray_right`` - vector. the vector from
          ``polytope_origin`` to the other neighboring vertices.

        OUTPUT:

        The element of the lattice Euclidean group that maps ``self``
        to ``polytope`` with given origin and left/right neighboring
        vertex. A
        :class:`~sage.geometry.polyhedron.lattice_euclidean_group_element.LatticePolytopesNotIsomorphicError`
        is raised if no such isomorphism exists.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.ppl_lattice_polytope import LatticePolytope_PPL
            sage: L1 = LatticePolytope_PPL((1,0),(0,1),(0,0))
            sage: L2 = LatticePolytope_PPL((1,0,3),(0,1,0),(0,0,1))
            sage: v0, v1, v2 = L2.vertices()
            sage: L1._find_cyclic_isomorphism_matching_edge(L2, v0, v1-v0, v2-v0)
            The map A*x+b with A=
            [ 0  1]
            [-1 -1]
            [ 1  3]
            b =
            (0, 1, 0)
        """
        from sage.geometry.polyhedron.lattice_euclidean_group_element import \
            LatticePolytopesNotIsomorphicError
        polytope_matrix = block_matrix(1, 2, [p_ray_left.column(),
                                              p_ray_right.column()])
        self_vertices = self.ordered_vertices()
        for i in range(len(self_vertices)):
            # three consecutive vertices
            v_left = self_vertices[(i+0) % len(self_vertices)]
            v_origin = self_vertices[(i+1) % len(self_vertices)]
            v_right = self_vertices[(i+2) % len(self_vertices)]
            r_left = v_left-v_origin
            r_right = v_right-v_origin
            self_matrix = block_matrix(1, 2, [r_left.column(),
                                              r_right.column()])
            A = self_matrix.solve_left(polytope_matrix)
            b = polytope_origin - A*v_origin
            try:
                A = matrix(ZZ, A)
                b = vector(ZZ, b)
            except TypeError:
                continue
            if A.elementary_divisors()[0:2] != [1, 1]:
                continue
            hom = LatticeEuclideanGroupElement(A, b)
            if hom(self) == polytope:
                return hom
        raise LatticePolytopesNotIsomorphicError('different polygons')

#.........这里部分代码省略.........
        sage: skew_hadamard_matrix(10,existence=True)
        False
        sage: skew_hadamard_matrix(12,existence=True)
        True
        sage: skew_hadamard_matrix(784,existence=True)
        True
        sage: skew_hadamard_matrix(10)
        Traceback (most recent call last):
        ...
        ValueError: A skew Hadamard matrix of order 10 does not exist
        sage: skew_hadamard_matrix(36)
        36 x 36 dense matrix over Integer Ring...
        sage: skew_hadamard_matrix(36)==skew_hadamard_matrix(36,skew_normalize=False)
        False
        sage: skew_hadamard_matrix(52)
        52 x 52 dense matrix over Integer Ring...
        sage: skew_hadamard_matrix(92)
        92 x 92 dense matrix over Integer Ring...
        sage: skew_hadamard_matrix(816)     # long time
        816 x 816 dense matrix over Integer Ring...
        sage: skew_hadamard_matrix(100)
        Traceback (most recent call last):
        ...
        ValueError: A skew Hadamard matrix of order 100 is not yet implemented.
        sage: skew_hadamard_matrix(100,existence=True)
        Unknown

    REFERENCES:

    .. [Ha83] \M. Hall,
      Combinatorial Theory,
      2nd edition,
      Wiley, 1983
    """
    def true():
        _skew_had_cache[n]=True
        return True
    M = None
    if existence and n in _skew_had_cache:
        return True
    if not(n % 4 == 0) and (n > 2):
        if existence:
            return False
        raise ValueError("A skew Hadamard matrix of order %s does not exist" % n)
    if n == 2:
        if existence:
            return true()
        M = matrix([[1, 1], [-1, 1]])
    elif n == 1:
        if existence:
            return true()
        M = matrix([1])
    elif is_prime_power(n - 1) and ((n - 1) % 4 == 3):
        if existence:
            return true()
        M = hadamard_matrix_paleyI(n, normalize=False)

    elif n % 8 == 0:
        if skew_hadamard_matrix(n//2,existence=True): # (Lemma 14.1.6 in [Ha83]_)
            if existence:
                return true()
            H = skew_hadamard_matrix(n//2,check=False)
            M = block_matrix([[H,H], [-H.T,H.T]])

        else: # try Williamson construction (Lemma 14.1.5 in [Ha83]_)
            for d in divisors(n)[2:-2]: # skip 1, 2, n/2, and n
                n1 = n//d
                if is_prime_power(d - 1) and (d % 4 == 0) and (n1 % 4 == 0)\
                    and skew_hadamard_matrix(n1,existence=True):
                    if existence:
                        return true()
                    H = skew_hadamard_matrix(n1, check=False)-I(n1)
                    U = matrix(ZZ, d, lambda i, j: -1 if i==j==0 else\
                                        1 if i==j==1 or (i>1 and j-1==d-i)\
                                          else 0)
                    A = block_matrix([[matrix([0]), matrix(ZZ,1,d-1,[1]*(d-1))],
                                      [ matrix(ZZ,d-1,1,[-1]*(d-1)),
                                        _helper_payley_matrix(d-1,zero_position=0)]])+I(d)
                    M = A.tensor_product(I(n1))+(U*A).tensor_product(H)
                    break
    if M is None: # try Williamson-Goethals-Seidel construction
        if GS_skew_hadamard_smallcases(n, existence=True):
            if existence:
                return true()
            M = GS_skew_hadamard_smallcases(n)

        else:
            if existence:
                return Unknown
            raise ValueError("A skew Hadamard matrix of order %s is not yet implemented." % n)
    if skew_normalize:
        dd = diagonal_matrix(M[0])
        M = dd*M*dd
    if check:
        assert is_hadamard_matrix(M, normalized=False, skew=True)
        if skew_normalize:
            from sage.modules.free_module_element import vector
            assert M[0]==vector([1]*n)
    _skew_had_cache[n]=True
    return M

    def cohomology_complex(self, m):
        r"""
        Return the "cohomology complex" `C^*(m)`

        See [Klyachko]_, equation 4.2.

        INPUT:

        - ``m`` -- tuple of integers or `M`-lattice point. A point in
          the dual lattice of the fan. Must be immutable.

        OUTPUT:

        The "cohomology complex" as a chain complex over the
        :meth:`base_ring`.

        EXAMPLES::

            sage: P3 = toric_varieties.P(3)
            sage: rays = [(1,0,0), (0,1,0), (0,0,1)]
            sage: F1 = FilteredVectorSpace(rays, {0:[0], 1:[2], 2:[1]})
            sage: F2 = FilteredVectorSpace(rays, {0:[1,2], 1:[0]})
            sage: r = P3.fan().rays()
            sage: V = P3.sheaves.Klyachko({r[0]:F1, r[1]:F2, r[2]:F2, r[3]:F2})
            sage: tau = Cone([(1,0,0), (0,1,0)])
            sage: sigma = Cone([(1, 0, 0)])
            sage: M = P3.fan().dual_lattice()
            sage: m = M(1, 1, 0); m.set_immutable()
            sage: V.cohomology_complex(m)
            Chain complex with at most 2 nonzero terms over Rational Field

            sage: F = CyclotomicField(3)
            sage: P3 = toric_varieties.P(3).change_ring(F)
            sage: V = P3.sheaves.Klyachko({r[0]:F1, r[1]:F2, r[2]:F2, r[3]:F2})
            sage: V.cohomology_complex(m)
            Chain complex with at most 2 nonzero terms over Cyclotomic
            Field of order 3 and degree 2
        """
        fan = self._variety.fan()
        C = fan.complex()
        CV = []
        F = self.base_ring()
        for dim in range(1,fan.dim()+1):
            codim = fan.dim() - dim
            d_C = C.differential(codim)
            d_V = []
            for j in range(0, d_C.ncols()):
                tau = fan(dim)[j]
                d_V_row = []
                for i in range(0, d_C.nrows()):
                    sigma = fan(dim-1)[i]
                    if sigma.is_face_of(tau):
                        pr = self.E_quotient_projection(sigma, tau, m)
                        d = d_C[i,j] * pr.matrix().transpose()
                    else:
                        E_sigma = self.E_quotient(sigma, m)
                        E_tau = self.E_quotient(tau, m)
                        d = zero_matrix(F, E_tau.dimension(), E_sigma.dimension())
                    d_V_row.append(d)
                d_V.append(d_V_row)
            d_V = block_matrix(d_V, ring=F)
            CV.append(d_V)
        from sage.homology.chain_complex import ChainComplex
        return ChainComplex(CV, base_ring=self.base_ring())

def rshcd_from_close_prime_powers(n):
    r"""
    Return a `(n^2,1)`-RSHCD when `n-1` and `n+1` are odd prime powers and `n=0\pmod{4}`.

    The construction implemented here appears in Theorem 4.3 from [GS70]_.

    Note that the authors of [SWW72]_ claim in Corollary 5.12 (page 342) to have
    proved the same result without the `n=0\pmod{4}` restriction with a *very*
    similar construction. So far, however, I (Nathann Cohen) have not been able
    to make it work.

    INPUT:

    - ``n`` -- an integer congruent to `0\pmod{4}`

    .. SEEALSO::

        :func:`regular_symmetric_hadamard_matrix_with_constant_diagonal`

    EXAMPLES::

        sage: from sage.combinat.matrices.hadamard_matrix import rshcd_from_close_prime_powers
        sage: rshcd_from_close_prime_powers(4)
        [-1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1]
        [-1 -1  1  1 -1 -1 -1 -1 -1  1  1 -1 -1  1 -1  1]
        [ 1  1 -1  1  1 -1 -1 -1 -1 -1  1 -1 -1 -1  1 -1]
        [-1  1  1 -1  1  1 -1 -1 -1 -1 -1  1 -1 -1 -1  1]
        [ 1 -1  1  1 -1  1  1 -1 -1 -1 -1 -1  1 -1 -1 -1]
        [-1 -1 -1  1  1 -1  1  1 -1 -1 -1  1 -1  1 -1 -1]
        [-1 -1 -1 -1  1  1 -1 -1  1 -1  1 -1  1  1 -1 -1]
        [ 1 -1 -1 -1 -1  1 -1 -1 -1  1 -1  1 -1  1  1 -1]
        [-1 -1 -1 -1 -1 -1  1 -1 -1 -1  1  1  1 -1  1  1]
        [ 1  1 -1 -1 -1 -1 -1  1 -1 -1 -1 -1  1  1 -1  1]
        [-1  1  1 -1 -1 -1  1 -1  1 -1 -1 -1 -1  1  1 -1]
        [-1 -1 -1  1 -1  1 -1  1  1 -1 -1 -1 -1 -1  1  1]
        [ 1 -1 -1 -1  1 -1  1 -1  1  1 -1 -1 -1 -1 -1  1]
        [-1  1 -1 -1 -1  1  1  1 -1  1  1 -1 -1 -1 -1 -1]
        [ 1 -1  1 -1 -1 -1 -1  1  1 -1  1  1 -1 -1 -1 -1]
        [-1  1 -1  1 -1 -1 -1 -1  1  1 -1  1  1 -1 -1 -1]

    REFERENCE:

    .. [SWW72] A Street, W. Wallis, J. Wallis,
      Combinatorics: Room squares, sum-free sets, Hadamard matrices.
      Lecture notes in Mathematics 292 (1972).
    """
    if n%4:
        raise ValueError("n(={}) must be congruent to 0 mod 4")

    a,b = sorted([n-1,n+1],key=lambda x:-x%4)
    Sa  = _helper_payley_matrix(a)
    Sb  = _helper_payley_matrix(b)
    U   = matrix(a,[[int(i+j == a-1) for i in range(a)] for j in range(a)])

    K = (U*Sa).tensor_product(Sb) + U.tensor_product(J(b)-I(b)) - J(a).tensor_product(I(b))

    F = lambda x:diagonal_matrix([-(-1)**i for i in range(x)])
    G = block_diagonal_matrix([J(1),I(a).tensor_product(F(b))])
    e = matrix(a*b,[1]*(a*b))
    H = block_matrix(2,[-J(1),e.transpose(),e,K])

    HH = G*H*G
    assert len(set(map(sum,HH))) == 1
    assert HH**2 == n**2*I(n**2)
    return HH

def hadamard_matrix_paleyII(n):
    """
    Implements the Paley type II construction.

    The Paley type II case corresponds to the case `p \cong 1 \mod{4}` for a
    prime `p` (see [Hora]_).

    EXAMPLES::

        sage: sage.combinat.matrices.hadamard_matrix.hadamard_matrix_paleyII(12).det()
        2985984
        sage: 12^6
        2985984

    We note that the method returns a normalised Hadamard matrix ::

        sage: sage.combinat.matrices.hadamard_matrix.hadamard_matrix_paleyII(12)
        [ 1  1| 1  1| 1  1| 1  1| 1  1| 1  1]
        [ 1 -1|-1  1|-1  1|-1  1|-1  1|-1  1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1| 1 -1| 1  1|-1 -1|-1 -1| 1  1]
        [ 1  1|-1 -1| 1 -1|-1  1|-1  1| 1 -1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1| 1  1| 1 -1| 1  1|-1 -1|-1 -1]
        [ 1  1| 1 -1|-1 -1| 1 -1|-1  1|-1  1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1|-1 -1| 1  1| 1 -1| 1