def ProjectiveGeometryDesign(n, d, F, algorithm=None):
    """
    Returns a projective geometry design.

    A projective geometry design of parameters `n,d,F` has for points the lines
    of `F^{n+1}`, and for blocks the `d+1`-dimensional subspaces of `F^{n+1}`,
    each of which contains `\\frac {|F|^{d+1}-1} {|F|-1}` lines.

    INPUT:

    - ``n`` is the projective dimension

    - ``d`` is the dimension of the subspaces of `P = PPn(F)` which
      make up the blocks.

    - ``F`` is a finite field.

    - ``algorithm`` -- set to ``None`` by default, which results in using Sage's
      own implementation. In order to use GAP's implementation instead (i.e. its
      ``PGPointFlatBlockDesign`` function) set ``algorithm="gap"``. Note that
      GAP's "design" package must be available in this case, and that it can be
      installed with the ``gap_packages`` spkg.

    EXAMPLES:

    The points of the following design are the `\\frac {2^{2+1}-1} {2-1}=7`
    lines of `\mathbb{Z}_2^{2+1}`. It has `7` blocks, corresponding to each
    2-dimensional subspace of `\mathbb{Z}_2^{2+1}`::

        sage: designs.ProjectiveGeometryDesign(2, 1, GF(2))
        Incidence structure with 7 points and 7 blocks
        sage: BD = designs.ProjectiveGeometryDesign(2, 1, GF(2), algorithm="gap") # optional - gap_packages (design package)
        sage: BD.is_block_design()                                     # optional - gap_packages (design package)
        (True, [2, 7, 3, 1])
    """
    q = F.order()
    from sage.interfaces.gap import gap, GapElement
    from sage.sets.set import Set
    if algorithm is None:
        V = VectorSpace(F, n+1)
        points = list(V.subspaces(1))
        flats = list(V.subspaces(d+1))
        blcks = []
        for p in points:
            b = []
            for i in range(len(flats)):
                if p.is_subspace(flats[i]):
                    b.append(i)
            blcks.append(b)
        v = (q**(n+1)-1)/(q-1)
        return BlockDesign(v, blcks, name="ProjectiveGeometryDesign")
    if algorithm == "gap":   # Requires GAP's Design
        gap.load_package("design")
        gap.eval("D := PGPointFlatBlockDesign( %s, %s, %d )"%(n,q,d))
        v = eval(gap.eval("D.v"))
        gblcks = eval(gap.eval("D.blocks"))
        gB = []
        for b in gblcks:
            gB.append([x-1 for x in b])
        return BlockDesign(v, gB, name="ProjectiveGeometryDesign")

def WittDesign(n):
    """
    INPUT:

    - ``n`` is in `9,10,11,12,21,22,23,24`.

    Wraps GAP Design's WittDesign. If ``n=24`` then this function returns the
    large Witt design `W_{24}`, the unique (up to isomorphism) `5-(24,8,1)`
    design. If ``n=12`` then this function returns the small Witt design
    `W_{12}`, the unique (up to isomorphism) `5-(12,6,1)` design.  The other
    values of `n` return a block design derived from these.

    .. NOTE:

        Requires GAP's Design package (included in the gap_packages Sage spkg).

    EXAMPLES::

        sage: BD = designs.WittDesign(9)             # optional - gap_packages (design package)
        sage: BD.is_t_design(return_parameters=True) # optional - gap_packages (design package)
        (True, (2, 9, 3, 1))
        sage: BD                             # optional - gap_packages (design package)
        Incidence structure with 9 points and 12 blocks
        sage: print BD                       # optional - gap_packages (design package)
        Incidence structure with 9 points and 12 blocks
    """
    from sage.interfaces.gap import gap, GapElement
    gap.load_package("design")
    gap.eval("B:=WittDesign(%s)"%n)
    v = eval(gap.eval("B.v"))
    gblcks = eval(gap.eval("B.blocks"))
    gB = []
    for b in gblcks:
       gB.append([x-1 for x in b])
    return BlockDesign(v, gB, name="WittDesign", check=True)

    def dual(self, algorithm=None):
        """
        Return the dual of the incidence structure.

        INPUT:

        - ``algorithm`` -- whether to use Sage's implementation
          (``algorithm=None``, default) or use GAP's (``algorithm="gap"``).

          .. NOTE::

              The ``algorithm="gap"`` option requires GAP's Design package
              (included in the gap_packages Sage spkg).

        EXAMPLES:

        The dual of a projective plane is a projective plane::

            sage: PP = designs.DesarguesianProjectivePlaneDesign(4)
            sage: PP.dual().is_t_design(return_parameters=True)
            (True, (2, 21, 5, 1))

        TESTS::

            sage: D = designs.IncidenceStructure(4, [[0,2],[1,2,3],[2,3]])
            sage: D
            Incidence structure with 4 points and 3 blocks
            sage: D.dual()
            Incidence structure with 3 points and 4 blocks
            sage: print D.dual(algorithm="gap")       # optional - gap_packages
            IncidenceStructure<points=[0, 1, 2], blocks=[[0], [0, 1, 2], [1], [1, 2]]>
            sage: blocks = [[0,1,2],[0,3,4],[0,5,6],[1,3,5],[1,4,6],[2,3,6],[2,4,5]]
            sage: BD = designs.IncidenceStructure(7, blocks, name="FanoPlane");
            sage: BD
            Incidence structure with 7 points and 7 blocks
            sage: print BD.dual(algorithm="gap")         # optional - gap_packages
            IncidenceStructure<points=[0, 1, 2, 3, 4, 5, 6], blocks=[[0, 1, 2], [0, 3, 4], [0, 5, 6], [1, 3, 5], [1, 4, 6], [2, 3, 6], [2, 4, 5]]>
            sage: BD.dual()
            Incidence structure with 7 points and 7 blocks

        REFERENCE:

        - Soicher, Leonard, Design package manual, available at
          http://www.gap-system.org/Manuals/pkg/design/htm/CHAP003.htm
        """
        if algorithm == "gap":
            from sage.interfaces.gap import gap
            gap.load_package("design")
            gD = self._gap_()
            gap.eval("DD:=DualBlockDesign("+gD+")")
            v = eval(gap.eval("DD.v"))
            gblcks = eval(gap.eval("DD.blocks"))
            gB = []
            for b in gblcks:
                gB.append([x-1 for x in b])
            return IncidenceStructure(range(v), gB, name=None, check=False)
        else:
            return IncidenceStructure(
                          incidence_matrix=self.incidence_matrix().transpose(),
                          check=False)

def AffineGeometryDesign(n, d, F):
    r"""
    Return an Affine Geometry Design.

    INPUT:

    - `n` (integer) -- the Euclidean dimension. The number of points is
      `v=|F^n|`.

    - `d` (integer) -- the dimension of the (affine) subspaces of `P = GF(q)^n`
      which make up the blocks.

    - `F` -- a Finite Field (i.e. ``FiniteField(17)``), or a prime power
      (i.e. an integer)

    `AG_{n,d} (F)`, as it is sometimes denoted, is a `2` - `(v, k, \lambda)`
    design of points and `d`- flats (cosets of dimension `n`) in the affine
    geometry `AG_n (F)`, where

    .. math::

             v = q^n,\  k = q^d ,
             \lambda =\frac{(q^{n-1}-1) \cdots (q^{n+1-d}-1)}{(q^{n-1}-1) \cdots (q-1)}.

    Wraps some functions used in GAP Design's ``PGPointFlatBlockDesign``.  Does
    *not* require GAP's Design package.

    EXAMPLES::

        sage: BD = designs.AffineGeometryDesign(3, 1, GF(2))
        sage: BD.is_t_design(return_parameters=True)
        (True, (2, 8, 2, 1))
        sage: BD = designs.AffineGeometryDesign(3, 2, GF(2))
        sage: BD.is_t_design(return_parameters=True)
        (True, (3, 8, 4, 1))

    With an integer instead of a Finite Field::

        sage: BD = designs.AffineGeometryDesign(3, 2, 4)
        sage: BD.is_t_design(return_parameters=True)
        (True, (2, 64, 16, 5))
    """
    try:
        q = int(F)
    except TypeError:
        q = F.order()

    from sage.interfaces.gap import gap, GapElement
    from sage.sets.set import Set
    gap.eval("V:=GaloisField(%s)^%s"%(q,n))
    gap.eval("points:=AsSet(V)")
    gap.eval("Subs:=AsSet(Subspaces(V,%s));"%d)
    gap.eval("CP:=Cartesian(points,Subs)")
    flats = gap.eval("flats:=List(CP,x->Sum(x))") # affine spaces
    gblcks = eval(gap.eval("Set(List(flats,f->Filtered([1..Length(points)],i->points[i] in f)));"))
    v = q**n
    gB = []
    for b in gblcks:
       gB.append([x-1 for x in b])
    return BlockDesign(v, gB, name="AffineGeometryDesign")

    def invariant_quadratic_form(self):
        """
        This wraps GAP's command "InvariantQuadraticForm". From the GAP
        documentation:

        INPUT:

        -  ``self`` - a matrix group G

        OUTPUT:

        -  ``Q`` - the matrix satisfying the property: The
           quadratic form q on the natural vector space V on which G acts is
           given by `q(v) = v Q v^t`, and the invariance under G is
           given by the equation `q(v) = q(v M)` for all
           `v \in V` and `M \in G`.

        EXAMPLES::

            sage: G = GO( 4, GF(7), 1)
            sage: G.invariant_quadratic_form()
            [0 1 0 0]
            [0 0 0 0]
            [0 0 3 0]
            [0 0 0 1]
        """
        F = self.base_ring()
        G = self._gap_init_()
        cmd = "r := InvariantQuadraticForm("+G+")"
        gap.eval(cmd)
        cmd = "r.matrix"
        Q = gap(cmd)
        return Q._matrix_(F)

    def dual_incidence_structure(self, algorithm=None):
        """
        Wraps GAP Design's DualBlockDesign (see [1]). The dual of a block
        design may not be a block design.

        Also can be called with ``dual_design``.

        .. NOTE:

        The algorithm="gap" option requires GAP's Design package (included in
        the gap_packages Sage spkg).

        EXAMPLES::

            sage: from sage.combinat.designs.block_design import BlockDesign
            sage: D = BlockDesign(4, [[0,2],[1,2,3],[2,3]], test=False)
            sage: D
            Incidence structure with 4 points and 3 blocks
            sage: D.dual_design()
            Incidence structure with 3 points and 4 blocks
            sage: print D.dual_design(algorithm="gap")       # optional - gap_design
            IncidenceStructure<points=[0, 1, 2], blocks=[[0], [0, 1, 2], [1], [1, 2]]>
            sage: BD = IncidenceStructure(range(7),[[0,1,2],[0,3,4],[0,5,6],[1,3,5],[1,4,6],[2,3,6],[2,4,5]], name="FanoPlane")
            sage: BD
            Incidence structure with 7 points and 7 blocks
            sage: print BD.dual_design(algorithm="gap")         # optional - gap_design
            IncidenceStructure<points=[0, 1, 2, 3, 4, 5, 6], blocks=[[0, 1, 2], [0, 3, 4], [0, 5, 6], [1, 3, 5], [1, 4, 6], [2, 3, 6], [2, 4, 5]]>
            sage: BD.dual_incidence_structure()
            Incidence structure with 7 points and 7 blocks

        REFERENCE:

        - Soicher, Leonard, Design package manual, available at
          http://www.gap-system.org/Manuals/pkg/design/htm/CHAP003.htm
        """
        from sage.interfaces.gap import gap, GapElement
        from sage.sets.set import Set
        from sage.misc.flatten import flatten
        from sage.combinat.designs.block_design import BlockDesign
        from sage.misc.functional import transpose

        if algorithm == "gap":
            gap.load_package("design")
            gD = self._gap_()
            gap.eval("DD:=DualBlockDesign(" + gD + ")")
            v = eval(gap.eval("DD.v"))
            gblcks = eval(gap.eval("DD.blocks"))
            gB = []
            for b in gblcks:
                gB.append([x - 1 for x in b])
            return IncidenceStructure(range(v), gB, name=None, test=False)
        pts = self.blocks()
        M = transpose(self.incidence_matrix())
        blks = self.points()
        return IncidenceStructure(pts, blks, M, name=None, test=False)

def AffineGeometryDesign(n, d, F):
    r"""
    INPUT:

    - ``n`` is the Euclidean dimension, so the number of points is

    - ``v`` is the number of points `v = |F^n|` (`F = GF(q)`, some `q`)

    - `d` is the dimension of the (affine) subspaces of `P = GF(q)^n` which make
      up the blocks.

    `AG_{n,d} (F)`, as it is sometimes denoted, is a `2` - `(v, k, \lambda)`
    design of points and `d`- flats (cosets of dimension `n`) in the affine
    geometry `AG_n (F)`, where

    .. math::

             v = q^n,\  k = q^d ,
             \lambda =\frac{(q^{n-1}-1) \cdots (q^{n+1-d}-1)}{(q^{n-1}-1) \cdots (q-1)}.

    Wraps some functions used in GAP Design's PGPointFlatBlockDesign.
    Does *not* require GAP's Design.

    EXAMPLES::

        sage: BD = designs.AffineGeometryDesign(3, 1, GF(2))
        sage: BD.parameters()
        (2, 8, 2, 2)
        sage: BD.is_block_design()
        (True, [2, 8, 2, 2])
        sage: BD = designs.AffineGeometryDesign(3, 2, GF(2))
        sage: BD.parameters()
        (2, 8, 4, 12)
        sage: BD.is_block_design()
        (True, [3, 8, 4, 4])
    """
    q = F.order()
    from sage.interfaces.gap import gap, GapElement
    from sage.sets.set import Set
    gap.eval("V:=GaloisField(%s)^%s"%(q,n))
    gap.eval("points:=AsSet(V)")
    gap.eval("Subs:=AsSet(Subspaces(V,%s));"%d)
    gap.eval("CP:=Cartesian(points,Subs)")
    flats = gap.eval("flats:=List(CP,x->Sum(x))") # affine spaces
    gblcks = eval(gap.eval("AsSortedList(List(flats,f->Filtered([1..Length(points)],i->points[i] in f)));"))
    v = q**n
    gB = []
    for b in gblcks:
       gB.append([x-1 for x in b])
    return BlockDesign(v, gB, name="AffineGeometryDesign")

def ProjectiveGeometryDesign(n, d, F, algorithm=None):
    """
    INPUT:

    - ``n`` is the projective dimension

    - ``v`` is the number of points `PPn(GF(q))`

    - ``d`` is the dimension of the subspaces of `P = PPn(GF(q))` which
      make up the blocks

    - ``b`` is the number of `d`-dimensional subspaces of `P`

    Wraps GAP Design's PGPointFlatBlockDesign. Does *not* require
    GAP's Design.

    EXAMPLES::

        sage: designs.ProjectiveGeometryDesign(2, 1, GF(2))
        Incidence structure with 7 points and 7 blocks
        sage: BD = designs.ProjectiveGeometryDesign(2, 1, GF(2), algorithm="gap") # optional - gap_packages (design package)
        sage: BD.is_block_design()                                     # optional - gap_packages (design package)
        (True, [2, 7, 3, 1])
    """
    q = F.order()
    from sage.interfaces.gap import gap, GapElement
    from sage.sets.set import Set
    if algorithm == None:
        V = VectorSpace(F, n+1)
        points = list(V.subspaces(1))
        flats = list(V.subspaces(d+1))
        blcks = []
        for p in points:
            b = []
            for i in range(len(flats)):
                if p.is_subspace(flats[i]):
                    b.append(i)
            blcks.append(b)
        v = (q**(n+1)-1)/(q-1)
        return BlockDesign(v, blcks, name="ProjectiveGeometryDesign")
    if algorithm == "gap":   # Requires GAP's Design
        gap.load_package("design")
        gap.eval("D := PGPointFlatBlockDesign( %s, %s, %d )"%(n,q,d))
        v = eval(gap.eval("D.v"))
        gblcks = eval(gap.eval("D.blocks"))
        gB = []
        for b in gblcks:
            gB.append([x-1 for x in b])
        return BlockDesign(v, gB, name="ProjectiveGeometryDesign")

    def points_from_gap(self):
        """
        Literally pushes this block design over to GAP and returns the
        points of that. Other than debugging, usefulness is unclear.

        REQUIRES: GAP's Design package.

        EXAMPLES::

            sage: from sage.combinat.designs.block_design import BlockDesign
            sage: BD = BlockDesign(7,[[0,1,2],[0,3,4],[0,5,6],[1,3,5],[1,4,6],[2,3,6],[2,4,5]])
            sage: BD.points_from_gap()      # optional - gap_packages (design package)
            doctest:1: DeprecationWarning: Unless somebody protests this method will be removed, as nobody seems to know why it is there.
            See http://trac.sagemath.org/14499 for details.
            [1, 2, 3, 4, 5, 6, 7]
        """
        from sage.misc.superseded import deprecation
        deprecation(14499, ('Unless somebody protests this method will be '
                            'removed, as nobody seems to know why it is there.'))
        from sage.interfaces.gap import gap, GapElement
        from sage.sets.set import Set
        gap.load_package("design")
        gD = self._gap_()
        gP = gap.eval("BlockDesignPoints("+gD+")").replace("..",",")
        return range(eval(gP)[0],eval(gP)[1]+1)

    def _gap_latex_(self):
        r"""
        Return the GAP latex version of this matrix.
        
        AUTHORS:

        - S. Kohl: Wrote the GAP function.
        
        EXAMPLES::
        
            sage: F = GF(3); MS = MatrixSpace(F,2,2)
            sage: gens = [MS([[1,0],[0,1]]),MS([[1,1],[0,1]])]
            sage: G = MatrixGroup(gens)
            sage: g = G([[1, 1], [0, 1]])
            sage: print g._gap_latex_()
            \left(\begin{array}{rr}%
            Z(3)^{0}&Z(3)^{0}\\%
            0*Z(3)&Z(3)^{0}\\%
            \end{array}\right)%
        
        Type view(g._latex_()) to see the object in an xdvi window
        (assuming you have latex and xdvi installed).
        """
        s1 = self.__mat._gap_init_()
        s2 = gap.eval("LaTeX(" + s1 + ")")
        return eval(s2)

 def points_from_gap(self):
     """
     Literally pushes this block design over to GAP and returns the
     points of that. Other than debugging, usefulness is unclear.
     
     REQUIRES: GAP's Design package.
     
     EXAMPLES::
     
         sage: from sage.combinat.designs.block_design import BlockDesign
         sage: BD = BlockDesign(7,[[0,1,2],[0,3,4],[0,5,6],[1,3,5],[1,4,6],[2,3,6],[2,4,5]])
         sage: BD.points_from_gap()      # optional - gap_packages (design package)
         [1, 2, 3, 4, 5, 6, 7]
     """
     from sage.interfaces.gap import gap, GapElement
     from sage.sets.set import Set
     gap.eval('LoadPackage("design")')
     gD = self._gap_()
     gP = gap.eval("BlockDesignPoints("+gD+")").replace("..",",")
     return range(eval(gP)[0],eval(gP)[1]+1)

def WittDesign(n):
    """
    Input: n is in 9,10,11,12,21,22,23,24.
    
    Wraps GAP Design's WittDesign. If n=24 then this function returns
    the large Witt design W24, the unique (up to isomorphism)
    5-(24,8,1) design. If n=12 then this function returns the small
    Witt design W12, the unique (up to isomorphism) 5-(12,6,1) design.
    The other values of n return a block design derived from these.
    
    REQUIRES: GAP's Design package.
    
    EXAMPLES::
    
        sage: BD = WittDesign(9)   # optional - gap_packages (design package)
        sage: BD.parameters()      # optional - gap_packages (design package)
        (2, 9, 3, 1)
        sage: BD                   # optional - gap_packages (design package)
        Incidence structure with 9 points and 12 blocks
        sage: print BD             # optional - gap_packages (design package)
        WittDesign<points=[0, 1, 2, 3, 4, 5, 6, 7, 8], blocks=[[0, 1, 7], [0, 2, 5], [0, 3, 4], [0, 6, 8], [1, 2, 6], [1, 3, 5], [1, 4, 8], [2, 3, 8], [2, 4, 7], [3, 6, 7], [4, 5, 6], [5, 7, 8]]>
        sage: BD = WittDesign(12)  # optional - gap_packages (design package)
        sage: BD.parameters(t=5)   # optional - gap_packages (design package)
        (5, 12, 6, 1)
    """
    from sage.interfaces.gap import gap, GapElement
    gap.eval('LoadPackage("design")')
    gap.eval("B:=WittDesign(%s)"%n)
    v = eval(gap.eval("B.v"))
    gblcks = eval(gap.eval("B.blocks"))
    gB = []
    for b in gblcks:
       gB.append([x-1 for x in b])
    return BlockDesign(v, gB, name="WittDesign", test=True)

 def __init__(self, homset, imgsH, check=True):
     MatrixGroupMorphism.__init__(self, homset)   # sets the parent
     G = homset.domain()
     H = homset.codomain()
     gaplist_gens = [gap(x) for x in G.gens()]
     gaplist_imgs = [gap(x) for x in imgsH]
     genss = '[%s]'%(','.join(str(v) for v in gaplist_gens))
     imgss = '[%s]'%(','.join(str(v) for v in gaplist_imgs))
     args = '%s, %s, %s, %s'%(G._gap_init_(), H._gap_init_(), genss, imgss)
     self._gap_str = 'GroupHomomorphismByImages(%s)'%args
     phi0 = gap(self)
     if gap.eval("IsGroupHomomorphism(%s)"%phi0.name())!="true":
         raise ValueError,"The map "+str(gensG)+"-->"+str(imgsH)+" isn't a homomorphism."

    def pushforward(self, J, *args,**kwds):
        """
        The image of an element or a subgroup.

        INPUT:

        ``J`` -- a subgroup or an element of the domain of ``self``.

        OUTPUT:

        The image of ``J`` under ``self``
        
        NOTE:

        ``pushforward`` is the method that is used when a map is called on
        anything that is not an element of its domain. For historical reasons,
        we keep the alias ``image()`` for this method.

        EXAMPLES::
        
            sage: F = GF(7); MS = MatrixSpace(F,2,2)
            sage: F.multiplicative_generator()
            3
            sage: G = MatrixGroup([MS([3,0,0,1])])
            sage: a = G.gens()[0]^2
            sage: phi = G.hom([a])
            sage: phi.image(G.gens()[0]) # indirect doctest
            [2 0]
            [0 1]
            sage: H = MatrixGroup([MS(a.list())])
            sage: H
            Matrix group over Finite Field of size 7 with 1 generators:
             [[[2, 0], [0, 1]]]

        The following tests against trac ticket #10659::

            sage: phi(H)   # indirect doctestest
            Matrix group over Finite Field of size 7 with 1 generators: 
             [[[4, 0], [0, 1]]]
        """
        phi = gap(self)
        F = self.codomain().base_ring()
        from sage.all import MatrixGroup
        gapJ = gap(J)
        if gap.eval("IsGroup(%s)"%gapJ.name()) == "true":
            return MatrixGroup([x._matrix_(F) for x in phi.Image(gapJ).GeneratorsOfGroup()])
        return phi.Image(gapJ)._matrix_(F)

 def character_table(self):
     """
     Returns the GAP character table as a string. For larger tables you
     may preface this with a command such as
     gap.eval("SizeScreen([120,40])") in order to widen the screen.
     
     EXAMPLES::
     
         sage: print WeylGroup(['A',3]).character_table()
         CT1
         <BLANKLINE>
              2  3  2  2  .  3
              3  1  .  .  1  .
         <BLANKLINE>
                1a 4a 2a 3a 2b
         <BLANKLINE>
         X.1     1 -1 -1  1  1
         X.2     3  1 -1  . -1
         X.3     2  .  . -1  2
         X.4     3 -1  1  . -1
         X.5     1  1  1  1  1
     """
     return gap.eval("Display(CharacterTable(%s))" % gap(self).name())

    def tuple_orbit_reps(self, k, prefix=None):

        if prefix is None:
            prefix = []

        s = len(prefix)
        tp = tuple(prefix)
        if s > k:
            raise ValueError

        if k == 0:
            return 1, {() : 1}

        gens = self._graph.automorphism_group_gens()

        # Pass generators to GAP to create a group for us.

        gen_str = ",".join("(" + "".join(str(cy) for cy in cys) + ")" for cys in gens)
        gap.eval("g := Group(%s);" % gen_str)
        if len(prefix) > 0:
            gap.eval("g := Stabilizer(g, %s, OnTuples);" % list(set(prefix)))

        S = []
        for i in range(1, k - s + 1):
            S.extend([tuple(sorted(list(x))) for x in Subsets(self._graph.n, i)])

        set_orb_reps = {}

        #sys.stdout.write("Calculating orbits")

        while len(S) > 0:

            rep = list(S[0])

            o = gap.new("Orbit(g, %s, OnSets);" % (rep,)).sage()
            o = list(set([tuple(sorted(t)) for t in o]))
            ot = o[0]
            set_orb_reps[ot] = len(o)
            for t in o:
                S.remove(t)
            #sys.stdout.write(".")
            #sys.stdout.flush()

        #sys.stdout.write("\n")

        combs = [tuple(c) for c in Compositions(k - s)]
        factors = []
        for c in combs:
            factor = factorial(k - s)
            for x in c:
                factor /= factorial(x)
            factors.append(factor)

        orb_reps = {}
        total = 0

        for ot, length in set_orb_reps.iteritems():

            ne = len(ot)
            for ci in range(len(combs)):
                c = combs[ci]
                if len(c) == ne:
                    t = tp
                    for i in range(ne):
                        t += c[i] * (ot[i],)
                    weight = factors[ci] * length
                    orb_reps[t] = weight
                    total += weight

        return total, orb_reps

    def module_composition_factors(self, algorithm=None):
        r"""
        Return a list of triples consisting of [base field, dimension,
        irreducibility], for each of the Meataxe composition factors
        modules. The ``algorithm="verbose"`` option returns more information,
        but in Meataxe notation.

        EXAMPLES::

            sage: F=GF(3);MS=MatrixSpace(F,4,4)
            sage: M=MS(0)
            sage: M[0,1]=1;M[1,2]=1;M[2,3]=1;M[3,0]=1
            sage: G = MatrixGroup([M])
            sage: G.module_composition_factors()
            [(Finite Field of size 3, 1, True),
             (Finite Field of size 3, 1, True),
             (Finite Field of size 3, 2, True)]
            sage: F = GF(7); MS = MatrixSpace(F,2,2)
            sage: gens = [MS([[0,1],[-1,0]]),MS([[1,1],[2,3]])]
            sage: G = MatrixGroup(gens)
            sage: G.module_composition_factors()
            [(Finite Field of size 7, 2, True)]

        Type ``G.module_composition_factors(algorithm='verbose')`` to get a
        more verbose version.

        For more on MeatAxe notation, see
        http://www.gap-system.org/Manuals/doc/ref/chap69.html
        """
        from sage.misc.sage_eval import sage_eval
        F = self.base_ring()
        if not(F.is_finite()):
            raise NotImplementedError("Base ring must be finite.")
        q = F.cardinality()
        gens = self.gens()
        n = self.degree()
        MS = MatrixSpace(F,n,n)
        mats = [] # initializing list of mats by which the gens act on self
        W = self.matrix_space().row_space()
        for g in gens:
            p = MS(g.matrix())
            m = p.rows()
            mats.append(m)
        mats_str = str(gap([[list(r) for r in m] for m in mats]))
        gap.eval("M:=GModuleByMats("+mats_str+", GF("+str(q)+"))")
        gap.eval("MCFs := MTX.CompositionFactors( M )")
        N = eval(gap.eval("Length(MCFs)"))
        if algorithm == "verbose":
            print(gap.eval('MCFs') + "\n")
        L = []
        for i in range(1,N+1):
            gap.eval("MCF := MCFs[%s]"%i)
            L.append(tuple([sage_eval(gap.eval("MCF.field")),
                            eval(gap.eval("MCF.dimension")),
                            sage_eval(gap.eval("MCF.IsIrreducible")) ]))
        return sorted(L)

    def as_permutation_group(self, algorithm=None):
        r"""
        Return a permutation group representation for the group.

        In most cases occurring in practice, this is a permutation
        group of minimal degree (the degree begin determined from
        orbits under the group action). When these orbits are hard to
        compute, the procedure can be time-consuming and the degree
        may not be minimal.

        INPUT:

        - ``algorithm`` -- ``None`` or ``'smaller'``. In the latter
          case, try harder to find a permutation representation of
          small degree.

        OUTPUT:

        A permutation group isomorphic to ``self``. The
        ``algorithm='smaller'`` option tries to return an isomorphic
        group of low degree, but is not guaranteed to find the
        smallest one.

        EXAMPLES::
        
            sage: MS = MatrixSpace(GF(2), 5, 5)
            sage: A = MS([[0,0,0,0,1],[0,0,0,1,0],[0,0,1,0,0],[0,1,0,0,0],[1,0,0,0,0]])
            sage: G = MatrixGroup([A])
            sage: G.as_permutation_group()
            Permutation Group with generators [(1,2)]
            sage: MS = MatrixSpace( GF(7), 12, 12)
            sage: GG = gap("ImfMatrixGroup( 12, 3 )")
            sage: GG.GeneratorsOfGroup().Length()
            3
            sage: g1 = MS(eval(str(GG.GeneratorsOfGroup()[1]).replace("\n","")))
            sage: g2 = MS(eval(str(GG.GeneratorsOfGroup()[2]).replace("\n","")))
            sage: g3 = MS(eval(str(GG.GeneratorsOfGroup()[3]).replace("\n","")))
            sage: G = MatrixGroup([g1, g2, g3])
            sage: G.cardinality()
            21499084800
            sage: set_random_seed(0); current_randstate().set_seed_gap()
            sage: P = G.as_permutation_group()
            sage: P.cardinality()
            21499084800
            sage: P.degree()  # random output
            144
            sage: set_random_seed(3); current_randstate().set_seed_gap()
            sage: Psmaller = G.as_permutation_group(algorithm="smaller")
            sage: Psmaller.cardinality()
            21499084800
            sage: Psmaller.degree()  # random output
            108

        In this case, the "smaller" option returned an isomorphic group of
        lower degree. The above example used GAP's library of irreducible
        maximal finite ("imf") integer matrix groups to construct the
        MatrixGroup G over GF(7). The section "Irreducible Maximal Finite
        Integral Matrix Groups" in the GAP reference manual has more
        details.

        TESTS::

            sage: A= matrix(QQ, 2, [0, 1, 1, 0])
            sage: B= matrix(QQ, 2, [1, 0, 0, 1])
            sage: a, b= MatrixGroup([A, B]).as_permutation_group().gens()
            sage: a.order(), b.order()
            (2, 1)
        """
        # Note that the output of IsomorphismPermGroup() depends on
        # memory locations and will change if you change the order of
        # doctests and/or architecture
        from sage.groups.perm_gps.permgroup import PermutationGroup
        if not self.is_finite():
            raise NotImplementedError("Group must be finite.")
        n = self.degree()
        MS = MatrixSpace(self.base_ring(), n, n)
        mats = [] # initializing list of mats by which the gens act on self
        for g in self.gens():
            p = MS(g.matrix())
            m = p.rows()
            mats.append(m)
        mats_str = str(gap([[list(r) for r in m] for m in mats]))
        gap.eval("iso:=IsomorphismPermGroup(Group("+mats_str+"))")
        if algorithm == "smaller":
            gap.eval("small:= SmallerDegreePermutationRepresentation( Image( iso ) );")
            C = gap("Image( small )")
        else:
            C = gap("Image( iso )")
        return PermutationGroup(gap_group=C, canonicalize=False)

#.........这里部分代码省略.........
    - ``d`` is the dimension of the subspaces which make up the blocks.

    - ``F`` -- a finite field or a prime power.

    - ``algorithm`` -- set to ``None`` by default, which results in using Sage's
      own implementation. In order to use GAP's implementation instead (i.e. its
      ``PGPointFlatBlockDesign`` function) set ``algorithm="gap"``. Note that
      GAP's "design" package must be available in this case, and that it can be
      installed with the ``gap_packages`` spkg.

    - ``point_coordinates`` -- ``True`` by default. Ignored and assumed to be ``False`` if
      ``algorithm="gap"``. If ``True``, the ground set is indexed by coordinates
      in `\GF{q}^{n+1}`.  Otherwise the ground set is indexed by integers.

    - ``check`` -- (optional, default to ``True``) whether to check the output.

    EXAMPLES:

    The set of `d`-dimensional subspaces in a `n`-dimensional projective space
    forms `2`-designs (or balanced incomplete block designs)::

        sage: PG = designs.ProjectiveGeometryDesign(4, 2, GF(2))
        sage: PG
        Incidence structure with 31 points and 155 blocks
        sage: PG.is_t_design(return_parameters=True)
        (True, (2, 31, 7, 7))

        sage: PG = designs.ProjectiveGeometryDesign(3, 1, GF(4))
        sage: PG.is_t_design(return_parameters=True)
        (True, (2, 85, 5, 1))

    Check with ``F`` being a prime power::

        sage: PG = designs.ProjectiveGeometryDesign(3, 2, 4)
        sage: PG
        Incidence structure with 85 points and 85 blocks

    Use coordinates::

        sage: PG = designs.ProjectiveGeometryDesign(2, 1, GF(3))
        sage: PG.blocks()[0]
        [(1, 0, 0), (1, 0, 1), (1, 0, 2), (0, 0, 1)]

    Use indexing by integers::

        sage: PG = designs.ProjectiveGeometryDesign(2,1,GF(3),point_coordinates=0)
        sage: PG.blocks()[0]
        [0, 1, 2, 12]

    Check that the constructor using gap also works::

        sage: BD = designs.ProjectiveGeometryDesign(2, 1, GF(2), algorithm="gap") # optional - gap_packages (design package)
        sage: BD.is_t_design(return_parameters=True)                              # optional - gap_packages (design package)
        (True, (2, 7, 3, 1))
    """
    try:
        q = int(F)
    except TypeError:
        q = F.cardinality()
    else:
        from sage.rings.finite_rings.finite_field_constructor import GF
        F = GF(q)

    if algorithm is None:
        from sage.matrix.echelon_matrix import reduced_echelon_matrix_iterator

        points = {p:i for i,p in enumerate(reduced_echelon_matrix_iterator(F,1,n+1,copy=True,set_immutable=True))}
        blocks = []
        for m1 in reduced_echelon_matrix_iterator(F,d+1,n+1,copy=False):
            b = []
            for m2 in reduced_echelon_matrix_iterator(F,1,d+1,copy=False):
                m = m2*m1
                m.echelonize()
                m.set_immutable()
                b.append(points[m])
            blocks.append(b)
        B = BlockDesign(len(points), blocks, name="ProjectiveGeometryDesign", check=check)
        if point_coordinates:
            B.relabel({i:p[0] for p,i in points.iteritems()})

    elif algorithm == "gap":   # Requires GAP's Design
        from sage.interfaces.gap import gap
        gap.load_package("design")
        gap.eval("D := PGPointFlatBlockDesign( %s, %s, %d )"%(n,F.order(),d))
        v = eval(gap.eval("D.v"))
        gblcks = eval(gap.eval("D.blocks"))
        gB = []
        for b in gblcks:
            gB.append([x-1 for x in b])
        B = BlockDesign(v, gB, name="ProjectiveGeometryDesign", check=check)

    if check:
        from sage.combinat.q_analogues import q_binomial
        q = F.cardinality()
        if not B.is_t_design(t=2, v=q_binomial(n+1,1,q),
                                  k=q_binomial(d+1,1,q),
                                  l=q_binomial(n-1, d-1, q)):
            raise RuntimeError("error in ProjectiveGeometryDesign "
                    "construction. Please e-mail [email protected]")
    return B

    def dual_incidence_structure(self, algorithm=None):
        ""&