def BIBD_5q_5_for_q_prime_power(q):
    r"""
    Return a `(5q,5,1)`-BIBD with `q\equiv 1\pmod 4` a prime power.

    See Theorem 24 [ClaytonSmith]_.

    INPUT:

    - ``q`` (integer) -- a prime power such that `q\equiv 1\pmod 4`.

    EXAMPLES::

        sage: from sage.combinat.designs.bibd import BIBD_5q_5_for_q_prime_power
        sage: for q in [25, 45, 65, 85, 125, 145, 185, 205, 305, 405, 605]: # long time
        ....:     _ = BIBD_5q_5_for_q_prime_power(q/5)                      # long time
    """
    from sage.rings.finite_rings.constructor import FiniteField

    if q%4 != 1 or not is_prime_power(q):
        raise ValueError("q is not a prime power or q%4!=1.")

    d = (q-1)/4
    B = []
    F = FiniteField(q,'x')
    a = F.primitive_element()
    L = {b:i for i,b in enumerate(F)}
    for b in L:
        B.append([i*q + L[b] for i in range(5)])
        for i in range(5):
            for j in range(d):
                B.append([        i*q + L[b          ],
                          ((i+1)%5)*q + L[ a**j+b    ],
                          ((i+1)%5)*q + L[-a**j+b    ],
                          ((i+4)%5)*q + L[ a**(j+d)+b],
                          ((i+4)%5)*q + L[-a**(j+d)+b],
                          ])

    return B

def hadamard_matrix_paleyI(n, normalize=True):
    """
    Implements the Paley type I construction.

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

    INPUT:

    - ``n`` -- the matrix size 

    - ``normalize`` (boolean) -- whether to normalize the result.

    EXAMPLES:

    We note that this method by default returns a normalised Hadamard matrix ::

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

    Otherwise, it returns a skew Hadamard matrix `H`, i.e. `H=S+I`, with
    `S=-S^\top`  ::

        sage: M=hadamard_matrix_paleyI(4, normalize=False); M
        [ 1  1  1  1]
        [-1  1  1 -1]
        [-1 -1  1  1]
        [-1  1 -1  1]
        sage: S=M-identity_matrix(4); -S==S.T
        True

    TESTS::

        sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix
        sage: test_cases = [x+1 for x in range(100) if is_prime_power(x) and x%4==3]
        sage: all(is_hadamard_matrix(hadamard_matrix_paleyI(n),normalized=True,verbose=True)
        ....:     for n in test_cases)
        True
        sage: all(is_hadamard_matrix(hadamard_matrix_paleyI(n,normalize=False),verbose=True)
        ....:     for n in test_cases)
        True
    """
    p = n - 1
    if not(is_prime_power(p) and (p % 4 == 3)):
        raise ValueError("The order %s is not covered by the Paley type I construction." % n)

    from sage.rings.finite_rings.finite_field_constructor import FiniteField
    K = FiniteField(p,'x')
    K_list = list(K)
    K_list.insert(0,K.zero())
    H = matrix(ZZ, [[(1 if (x-y).is_square() else -1)
                     for x in K_list]
                    for y in K_list])
    for i in range(n):
        H[i,0] = -1
        H[0,i] =  1
    if normalize:
        for i in range(n):
            H[i,i] = -1
        H = normalise_hadamard(H)
    return H

#.........这里部分代码省略.........
        144 x 144 dense matrix over Integer Ring

    REFERENCE:

    .. [BH12] \A. Brouwer and W. Haemers,
      Spectra of graphs,
      Springer, 2012,
      http://homepages.cwi.nl/~aeb/math/ipm/ipm.pdf

    .. [HX10] \W. Haemers and Q. Xiang,
      Strongly regular graphs with parameters `(4m^4,2m^4+m^2,m^4+m^2,m^4+m^2)` exist for all `m>1`,
      European Journal of Combinatorics,
      Volume 31, Issue 6, August 2010, Pages 1553-1559,
      http://dx.doi.org/10.1016/j.ejc.2009.07.009.
    """
    if existence and (n,e) in _rshcd_cache:
        return _rshcd_cache[n,e]

    from sage.graphs.strongly_regular_db import strongly_regular_graph

    def true():
        _rshcd_cache[n,e] = True
        return True

    M = None
    if abs(e) != 1:
        raise ValueError
    if n<0:
        if existence:
            return False
        raise ValueError
    elif n == 4:
        if existence:
            return true()
        if e == 1:
            M = J(4)-2*matrix(4,[[int(i+j == 3) for i in range(4)] for j in range(4)])
        else:
            M = -J(4)+2*I(4)
    elif n ==  36:
        if existence:
            return true()
        if e == 1:
            M = strongly_regular_graph(36, 15, 6, 6).adjacency_matrix()
            M = J(36) - 2*M
        else:
            M = strongly_regular_graph(36,14,4,6).adjacency_matrix()
            M =  -J(36) + 2*M + 2*I(36)
    elif n == 100:
        if existence:
            return true()
        if e == -1:
            M = strongly_regular_graph(100,44,18,20).adjacency_matrix()
            M = 2*M - J(100) + 2*I(100)
        else:
            M = strongly_regular_graph(100,45,20,20).adjacency_matrix()
            M = J(100) - 2*M
    elif n == 196 and e == 1:
        if existence:
            return true()
        M = strongly_regular_graph(196,91,42,42).adjacency_matrix()
        M = J(196) - 2*M
    elif n == 324:
        if existence:
            return true()
        M = RSHCD_324(e)
    elif (  e  == 1                 and
          n%16 == 0                 and
          is_square(n)              and
          is_prime_power(sqrt(n)-1) and
          is_prime_power(sqrt(n)+1)):
        if existence:
            return true()
        M = -rshcd_from_close_prime_powers(int(sqrt(n)))

    # Recursive construction: the kronecker product of two RSHCD is a RSHCD
    else:
        from itertools import product
        for n1,e1 in product(divisors(n)[1:-1],[-1,1]):
            e2 = e1*e
            n2 = n//n1
            if (regular_symmetric_hadamard_matrix_with_constant_diagonal(n1,e1,existence=True) and
                regular_symmetric_hadamard_matrix_with_constant_diagonal(n2,e2,existence=True)):
                if existence:
                    return true()
                M1 = regular_symmetric_hadamard_matrix_with_constant_diagonal(n1,e1)
                M2 = regular_symmetric_hadamard_matrix_with_constant_diagonal(n2,e2)
                M  = M1.tensor_product(M2)
                break

    if M is None:
        from sage.misc.unknown import Unknown
        _rshcd_cache[n,e] = Unknown
        if existence:
            return Unknown
        raise ValueError("I do not know how to build a {}-RSHCD".format((n,e)))

    assert M*M.transpose() == n*I(n)
    assert set(map(sum,M)) == {e*sqrt(n)}

    return M

def hadamard_matrix(n,existence=False, check=True):
    r"""
    Tries to construct a Hadamard matrix using a combination of Paley
    and Sylvester constructions.

    INPUT:

    - ``n`` (integer) -- dimension of the matrix

    - ``existence`` (boolean) -- whether to build the matrix or merely query if
      a construction is available in Sage. When set to ``True``, the function
      returns:

        - ``True`` -- meaning that Sage knows how to build the matrix

        - ``Unknown`` -- meaning that Sage does not know how to build the
          matrix, although the matrix may exist (see :mod:`sage.misc.unknown`).

        - ``False`` -- meaning that the matrix does not exist.

    - ``check`` (boolean) -- whether to check that output is correct before
      returning it. As this is expected to be useless (but we are cautious
      guys), you may want to disable it whenever you want speed. Set to ``True``
      by default.

    EXAMPLES::

        sage: hadamard_matrix(12).det()
        2985984
        sage: 12^6
        2985984
        sage: hadamard_matrix(1)
        [1]
        sage: hadamard_matrix(2)
        [ 1  1]
        [ 1 -1]
        sage: hadamard_matrix(8) # random
        [ 1  1  1  1  1  1  1  1]
        [ 1 -1  1 -1  1 -1  1 -1]
        [ 1  1 -1 -1  1  1 -1 -1]
        [ 1 -1 -1  1  1 -1 -1  1]
        [ 1  1  1  1 -1 -1 -1 -1]
        [ 1 -1  1 -1 -1  1 -1  1]
        [ 1  1 -1 -1 -1 -1  1  1]
        [ 1 -1 -1  1 -1  1  1 -1]
        sage: hadamard_matrix(8).det() == 8^4
        True

    We note that the method `hadamard_matrix()` returns a normalised Hadamard matrix
    (the entries in the first row and column are all +1) ::

        sage: hadamard_matrix(12) # random
        [ 1  1| 1  1| 1  1| 1  1| 1  1| 1  1]
        [ 1 -1|-1  1|-1  1|-1  1|-1  1|-1  1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1| 1 -1| 1  1|-1 -1|-1 -1| 1  1]
        [ 1  1|-1 -1| 1 -1|-1  1|-1  1| 1 -1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1| 1  1| 1 -1| 1  1|-1 -1|-1 -1]
        [ 1  1| 1 -1|-1 -1| 1 -1|-1  1|-1  1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1|-1 -1| 1  1| 1 -1| 1  1|-1 -1]
        [ 1  1|-1  1| 1 -1|-1 -1| 1 -1|-1  1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1|-1 -1|-1 -1| 1  1| 1 -1| 1  1]
        [ 1  1|-1  1|-1  1| 1 -1|-1 -1| 1 -1]
        [-----+-----+-----+-----+-----+-----]
        [ 1 -1| 1  1|-1 -1|-1 -1| 1  1| 1 -1]
        [ 1  1| 1 -1|-1  1|-1  1| 1 -1|-1 -1]

    TESTS::

        sage: matrix.hadamard(10,existence=True)
        False
        sage: matrix.hadamard(12,existence=True)
        True
        sage: matrix.hadamard(92,existence=True)
        Unknown
        sage: matrix.hadamard(10)
        Traceback (most recent call last):
        ...
        ValueError: The Hadamard matrix of order 10 does not exist
    """
    if not(n % 4 == 0) and (n > 2):
        if existence:
            return False
        raise ValueError("The 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//2 - 1) and (n//2 - 1) % 4 == 1:
        if existence:
            return True
        M = hadamard_matrix_paleyII(n)
    elif n == 4 or n % 8 == 0:
#.........这里部分代码省略.........

#.........这里部分代码省略.........
         2 2
         3 3
         4 4
         5 5
         6 2
         7 7
         8 8
         9 9
        10 2
        11 11
        12 6
        13 13
        14 2
        15 3
        16 16
        17 17
        18 2
        19 19
        20 4
        21 6
        22 2
        23 23
        24 8
        25 25
        26 2
        27 27
        28 6
        29 29

    TESTS::

        sage: designs.difference_matrix(10,12,1,existence=True)
        False
        sage: designs.difference_matrix(10,12,1)
        Traceback (most recent call last):
        ...
        EmptySetError: No (10,12,1)-Difference Matrix exists as k(=12)>g(=10)
        sage: designs.difference_matrix(10,9,1,existence=True)
        Unknown
        sage: designs.difference_matrix(10,9,1)
        Traceback (most recent call last):
        ...
        NotImplementedError: I don't know how to build a (10,9,1)-Difference Matrix!
    """

    if lmbda == 1 and k is not None and k > g:
        if existence:
            return False
        raise EmptySetError("No ({},{},{})-Difference Matrix exists as k(={})>g(={})".format(g, k, lmbda, k, g))

    # Prime powers
    elif lmbda == 1 and is_prime_power(g):
        if k is None:
            if existence:
                return g
            else:
                k = g
        elif existence:
            return True
        F = FiniteField(g, "x")
        F_set = list(F)
        F_k_set = F_set[:k]

        G = F
        M = [[x * y for y in F_k_set] for x in F_set]

    # Treat the case k=None
    # (find the max k such that there exists a DM)
    elif k is None:
        i = 2
        while difference_matrix(g=g, k=i, lmbda=lmbda, existence=True):
            i += 1
        return i - 1

    # From the database
    elif (g, lmbda) in DM_constructions and DM_constructions[g, lmbda][0] >= k:
        if existence:
            return True
        _, f = DM_constructions[g, lmbda]
        G, M = f()
        M = [R[:k] for R in M]

    # Product construction
    elif find_product_decomposition(g, k, lmbda):
        if existence:
            return True
        (g1, lmbda1), (g2, lmbda2) = find_product_decomposition(g, k, lmbda)
        G1, M1 = difference_matrix(g1, k, lmbda1)
        G2, M2 = difference_matrix(g2, k, lmbda2)
        G, M = difference_matrix_product(k, M1, G1, lmbda1, M2, G2, lmbda2, check=False)

    else:
        if existence:
            return Unknown
        raise NotImplementedError("I don't know how to build a ({},{},{})-Difference Matrix!".format(g, k, lmbda))

    if check and not is_difference_matrix(M, G, k, lmbda, 1):
        raise RuntimeError("Sage built something which is not a ({},{},{})-DM!".format(g, k, lmbda))

    return G, M

def AhrensSzekeresGeneralizedQuadrangleGraph(q, dual=False):
    r"""
    Return the collinearity graph of the generalized quadrangle `AS(q)`, or of its dual

    Let `q` be an odd prime power.  `AS(q)` is a generalized quadrangle [GQwiki]_ of
    order `(q-1,q+1)`, see 3.1.5 in [PT09]_. Its points are elements
    of `F_q^3`, and lines are sets of size `q` of the form

    * `\{ (\sigma, a, b) \mid \sigma\in F_q \}`
    * `\{ (a, \sigma, b) \mid \sigma\in F_q \}`
    * `\{ (c \sigma^2 - b \sigma + a, -2 c \sigma + b, \sigma) \mid \sigma\in F_q \}`,

    where `a`, `b`, `c` are arbitrary elements of `F_q`.

    INPUT:

    - ``q`` -- a power of an odd prime number

    - ``dual`` -- if ``False`` (default), return the collinearity graph of `AS(q)`.
      Otherwise return the collinearity graph of the dual `AS(q)`.

    EXAMPLES::

        sage: g=graphs.AhrensSzekeresGeneralizedQuadrangleGraph(5); g
        AS(5); GQ(4, 6): Graph on 125 vertices
        sage: g.is_strongly_regular(parameters=True)
        (125, 28, 3, 7)
        sage: g=graphs.AhrensSzekeresGeneralizedQuadrangleGraph(5,dual=True); g
        AS(5)*; GQ(6, 4): Graph on 175 vertices
        sage: g.is_strongly_regular(parameters=True)
        (175, 30, 5, 5)

    REFERENCE:

    .. [GQwiki] `Generalized quadrangle
      <http://en.wikipedia.org/wiki/Generalized_quadrangle>`__

    .. [PT09] \S. Payne, J. A. Thas.
      Finite generalized quadrangles.
      European Mathematical Society,
      2nd edition, 2009.
    """
    from sage.combinat.designs.incidence_structures import IncidenceStructure
    p, k = is_prime_power(q,get_data=True)
    if k==0 or p==2:
       raise ValueError('q must be an odd prime power')
    F = FiniteField(q, 'a')
    L = []
    for a in F:
        for b in F:
            L.append(tuple(map(lambda s: (s, a, b), F)))
            L.append(tuple(map(lambda s: (a, s, b), F)))
            for c in F:
                L.append(tuple(map(lambda s: (c*s**2 - b*s + a, -2*c*s + b, s), F)))
    if dual:
        G = IncidenceStructure(L).intersection_graph()
        G.name('AS('+str(q)+')*; GQ'+str((q+1,q-1)))
    else:
        G = IncidenceStructure(L).dual().intersection_graph()
        G.name('AS('+str(q)+'); GQ'+str((q-1,q+1)))
    return G

def NonisotropicUnitaryPolarGraph(m, q):
    r"""
    Returns the Graph `NU(m,q)`.

    Returns the graph on nonisotropic, with respect to a nondegenerate
    Hermitean form, points of the `(m-1)`-dimensional projective space over `F_q`,
    with points adjacent whenever they lie on a tangent (to the set of isotropic points)
    line.
    For more information, see Sect. 9.9 of [BH12]_ and series C14 in [Hu75]_.

    INPUT:

    - ``m,q`` (integers) -- `q` must be a prime power.

    EXAMPLES::

        sage: g=graphs.NonisotropicUnitaryPolarGraph(5,2); g
        NU(5, 2): Graph on 176 vertices
        sage: g.is_strongly_regular(parameters=True)
        (176, 135, 102, 108)

    TESTS::

        sage: graphs.NonisotropicUnitaryPolarGraph(4,2).is_strongly_regular(parameters=True)
        (40, 27, 18, 18)
        sage: graphs.NonisotropicUnitaryPolarGraph(4,3).is_strongly_regular(parameters=True) # long time
        (540, 224, 88, 96)
        sage: graphs.NonisotropicUnitaryPolarGraph(6,6)
        Traceback (most recent call last):
        ...
        ValueError: q must be a prime power

    REFERENCE:

    .. [Hu75] \X. L. Hubaut.
      Strongly regular graphs.
      Disc. Math. 13(1975), pp 357--381.
      http://dx.doi.org/10.1016/0012-365X(75)90057-6
    """
    p, k = is_prime_power(q,get_data=True)
    if k==0:
       raise ValueError('q must be a prime power')
    from sage.libs.gap.libgap import libgap
    from itertools import combinations
    F=libgap.GF(q**2)  # F_{q^2}
    W=libgap.FullRowSpace(F, m)  # F_{q^2}^m
    B=libgap.Elements(libgap.Basis(W))      # the standard basis of W
    if m % 2 != 0:
        point = B[(m-1)/2]
    else:
        if p==2:
            point = B[m/2] + F.PrimitiveRoot()*B[(m-2)/2]
        else:
            point = B[(m-2)/2] + B[m/2]
    g = libgap.GeneralUnitaryGroup(m,q)
    V = libgap.Orbit(g,point,libgap.OnLines) # orbit on nonisotropic points
    gp = libgap.Action(g,V,libgap.OnLines)  # make a permutation group

    s = libgap.Subspace(W,[point, point+B[0]]) # a tangent line on point

    # and the points there
    sp = [libgap.Elements(libgap.Basis(x))[0] for x in libgap.Elements(s.Subspaces(1))]
    h = libgap.Set(map(lambda x: libgap.Position(V, x), libgap.Intersection(V,sp))) # indices
    L = libgap.Orbit(gp, h, libgap.OnSets) # orbit on the tangent lines
    G = Graph()
    for x in L: # every pair of points in the subspace is adjacent to each other in G
        G.add_edges(combinations(x, 2))
    G.relabel()
    G.name("NU" + str((m, q)))
    return G

def CossidentePenttilaGraph(q):
    r"""
    Cossidente-Penttila `((q^3+1)(q+1)/2,(q^2+1)(q-1)/2,(q-3)/2,(q-1)^2/2)`-strongly regular graph

    For each odd prime power `q`, one can partition the points of the `O_6^-(q)`-generalized
    quadrange `GQ(q,q^2)` into two parts, so that on any of them the induced subgraph of
    the point graph of the GQ has parameters as above [CP05]_.

    Directly follwing the construction in [CP05]_ is not efficient,
    as one then needs to construct the dual `GQ(q^2,q)`. Thus we
    describe here a more efficient approach that we came up with, following a suggestion by
    T.Penttila. Namely, this partition is invariant
    under the subgroup `H=\Omega_3(q^2)<O_6^-(q)`. We build the appropriate `H`, which
    leaves the form `B(X,Y,Z)=XY+Z^2` invariant, and
    pick up two orbits of `H` on the `F_q`-points. One them is `B`-isotropic, and we
    take the representative `(1:0:0)`. The other one corresponds to the points of
    `PG(2,q^2)` that have all the lines on them either missing the conic specified by `B`, or
    intersecting the conic in two points. We take `(1:1:e)` as the representative. It suffices
    to pick `e` so that `e^2+1` is not a square in `F_{q^2}`. Indeed,
    The conic can be viewed as the union of `\{(0:1:0)\}` and `\{(1:-t^2:t) | t \in F_{q^2}\}`.
    The coefficients of a generic line on `(1:1:e)` are `[1:-1-eb:b]`, for `-1\neq eb`.
    Thus, to make sure the intersection with the conic is always even, we need that the
    discriminant of `1+(1+eb)t^2+tb=0` never vanishes, and this is if and only if
    `e^2+1` is not a square. Further, we need to adjust `B`, by multiplying it by appropriately
    chosen `\nu`, so that `(1:1:e)` becomes isotropic under the relative trace norm
    `\nu B(X,Y,Z)+(\nu B(X,Y,Z))^q`. The latter is used then to define the graph.

    INPUT:

    - ``q`` -- an odd prime power.

    EXAMPLES:

    For `q=3` one gets Sims-Gewirtz graph. ::

        sage: G=graphs.CossidentePenttilaGraph(3)    # optional - gap_packages (grape)
        sage: G.is_strongly_regular(parameters=True) # optional - gap_packages (grape)
        (56, 10, 0, 2)

    For `q>3` one gets new graphs. ::

        sage: G=graphs.CossidentePenttilaGraph(5)    # optional - gap_packages (grape)
        sage: G.is_strongly_regular(parameters=True) # optional - gap_packages (grape)
        (378, 52, 1, 8)

    TESTS::

        sage: G=graphs.CossidentePenttilaGraph(7)    # optional - gap_packages (grape) # long time
        sage: G.is_strongly_regular(parameters=True) # optional - gap_packages (grape) # long time
        (1376, 150, 2, 18)
        sage: graphs.CossidentePenttilaGraph(2)
        Traceback (most recent call last):
        ...
        ValueError: q(=2) must be an odd prime power

    REFERENCES:

    .. [CP05] \A.Cossidente and T.Penttila
       Hemisystems on the Hermitian surface
       Journal of London Math. Soc. 72(2005), 731-741
    """
    p, k = is_prime_power(q,get_data=True)
    if k==0 or p==2:
        raise ValueError('q(={}) must be an odd prime power'.format(q))

    from sage.libs.gap.libgap import libgap
    from sage.misc.package import is_package_installed, PackageNotFoundError

    if not is_package_installed('gap_packages'):
        raise PackageNotFoundError('gap_packages')

    adj_list=libgap.function_factory("""function(q)
        local z, e, so, G, nu, G1, G0, B, T, s, O1, O2, x;
        LoadPackage("grape");
        G0:=SO(3,q^2);
        so:=GeneratorsOfGroup(G0);
        G1:=Group(Comm(so[1],so[2]),Comm(so[1],so[3]),Comm(so[2],so[3]));
        B:=InvariantBilinearForm(G0).matrix;
        z:=Z(q^2); e:=z; sqo:=(q^2-1)/2;
        if IsInt(sqo/Order(e^2+z^0)) then
            e:=z^First([2..q^2-2], x-> not IsInt(sqo/Order(z^(2*x)+z^0)));
        fi;
        nu:=z^First([0..q^2-2], x->z^x*(e^2+z^0)+(z^x*(e^2+z^0))^q=0*z);
        T:=function(x)
            local r;
            r:=nu*x*B*x;
            return r+r^q;
        end;
        s:=Group([Z(q)*IdentityMat(3,GF(q))]);
        O1:=Orbit(G1, Set(Orbit(s,z^0*[1,0,0])), OnSets);
        O2:=Orbit(G1, Set(Orbit(s,z^0*[1,1,e])), OnSets);
        G:=Graph(G1,Concatenation(O1,O2),OnSets,
            function(x,y) return x<>y and 0*z=T(x[1]+y[1]); end);
        return List([1..OrderGraph(G)],x->Adjacency(G,x));
        end;""")

    adj = adj_list(q) # for each vertex, we get the list of vertices it is adjacent to
    G = Graph(((i,int(j-1))
               for i,ni in enumerate(adj) for j in ni),
               format='list_of_edges', multiedges=False)
#.........这里部分代码省略.........

def rshcd_from_prime_power_and_conference_matrix(n):
    r"""
    Return a `((n-1)^2,1)`-RSHCD if `n` is prime power, and symmetric `(n-1)`-conference matrix exists

    The construction implemented here is Theorem 16 (and Corollary 17) from [WW72]_.

    In [SWW72]_ this construction (Theorem 5.15 and Corollary 5.16)
    is reproduced with a typo. Note that [WW72]_ refers to [Sz69]_ for the construction,
    provided by :func:`szekeres_difference_set_pair`,
    of complementary difference sets, and the latter has a typo.

    From a :func:`symmetric_conference_matrix`, we only need the Seidel
    adjacency matrix of the underlying strongly regular conference (i.e. Paley
    type) graph, which we construct directly.

    INPUT:

    - ``n`` -- an integer

    .. SEEALSO::

        :func:`regular_symmetric_hadamard_matrix_with_constant_diagonal`

    EXAMPLES:

    A 36x36 example ::

        sage: from sage.combinat.matrices.hadamard_matrix import rshcd_from_prime_power_and_conference_matrix
        sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix
        sage: H = rshcd_from_prime_power_and_conference_matrix(7); H
        36 x 36 dense matrix over Integer Ring (use the '.str()' method to see the entries)
        sage: H==H.T and is_hadamard_matrix(H) and H.diagonal()==[1]*36 and list(sum(H))==[6]*36
        True

    Bigger examples, only provided by this construction ::

        sage: H = rshcd_from_prime_power_and_conference_matrix(27)  # long time
        sage: H == H.T and is_hadamard_matrix(H)                    # long time
        True
        sage: H.diagonal()==[1]*676 and list(sum(H))==[26]*676      # long time
        True

    In this example the conference matrix is not Paley, as 45 is not a prime power ::

        sage: H = rshcd_from_prime_power_and_conference_matrix(47)  # not tested (long time)

    REFERENCE:

    .. [WW72] \J. Wallis and A.L. Whiteman,
      Some classes of Hadamard matrices with constant diagonal,
      Bull. Austral. Math. Soc. 7(1972), 233-249
    """
    from sage.graphs.strongly_regular_db import strongly_regular_graph as srg
    if is_prime_power(n) and 2==(n-1)%4:
        try:
            M = srg(n-2,(n-3)//2,(n-7)//4)
        except ValueError:
            return
        m = (n-3)//4
        Q,X,Y = szekeres_difference_set_pair(m)
        B = typeI_matrix_difference_set(Q,X)
        A = -typeI_matrix_difference_set(Q,Y) # must be symmetric
        W = M.seidel_adjacency_matrix()
        f = J(1,4*m+1)
        e = J(1,2*m+1)
        JJ = J(2*m+1, 2*m+1)
        II = I(n-2)
        Ib = I(2*m+1)
        J4m = J(4*m+1,4*m+1)
        H34 = -(B+Ib).tensor_product(W)+Ib.tensor_product(J4m)+(Ib-JJ).tensor_product(II)
        A_t_W = A.tensor_product(W)
        e_t_f = e.tensor_product(f)
        H = block_matrix([
            [J(1,1),                 f,                      e_t_f,                  -e_t_f],
            [f.T,                  J4m,     e.tensor_product(W-II),  e.tensor_product(W+II)],
            [ e_t_f.T, (e.T).tensor_product(W-II), A_t_W+JJ.tensor_product(II),         H34],
            [-e_t_f.T, (e.T).tensor_product(W+II), H34.T,      -A_t_W+JJ.tensor_product(II)]])
        return H

def GDD_4_2(q, existence=False, check=True):
    r"""
    Return a `(2q,\{4\},\{2\})`-GDD for `q` a prime power with `q\equiv 1\pmod{6}`.

    This method implements Lemma VII.5.17 from [BJL99] (p.495).

    INPUT:

    - ``q`` (integer)

    - ``existence`` (boolean) -- instead of building the design, return:

        - ``True`` -- meaning that Sage knows how to build the design

        - ``Unknown`` -- meaning that Sage does not know how to build the
          design, but that the design may exist (see :mod:`sage.misc.unknown`).

        - ``False`` -- meaning that the design does not exist.

    - ``check`` -- (boolean) Whether to check that output is correct before
      returning it. As this is expected to be useless (but we are cautious
      guys), you may want to disable it whenever you want speed. Set to ``True``
      by default.

    EXAMPLE::

        sage: from sage.combinat.designs.group_divisible_designs import GDD_4_2
        sage: GDD_4_2(7,existence=True)
        True
        sage: GDD_4_2(7)
        Group Divisible Design on 14 points of type 2^7
        sage: GDD_4_2(8,existence=True)
        Unknown
        sage: GDD_4_2(8)
        Traceback (most recent call last):
        ...
        NotImplementedError
    """
    if q <= 1 or q % 6 != 1 or not is_prime_power(q):
        if existence:
            return Unknown
        raise NotImplementedError
    if existence:
        return True

    from sage.rings.finite_rings.finite_field_constructor import FiniteField as GF

    G = GF(q, "x")
    w = G.primitive_element()
    e = w ** ((q - 1) // 3)

    # A first parallel class is defined. G acts on it, which yields all others.
    first_class = [[(0, 0), (1, w ** i), (1, e * w ** i), (1, e * e * w ** i)] for i in range((q - 1) // 6)]

    label = {p: i for i, p in enumerate(G)}
    classes = [[[2 * label[x[1] + g] + (x[0] + j) % 2 for x in S] for S in first_class] for g in G for j in range(2)]

    return GroupDivisibleDesign(
        2 * q,
        groups=[[i, i + 1] for i in range(0, 2 * q, 2)],
        blocks=sum(classes, []),
        K=[4],
        G=[2],
        check=check,
        copy=False,
    )

def projective_plane(n, check=True, existence=False):
    r"""
    Return a projective plane of order ``n`` as a 2-design.

    A finite projective plane is a 2-design with `n^2+n+1` lines (or blocks) and
    `n^2+n+1` points. For more information on finite projective planes, see the
    :wikipedia:`Projective_plane#Finite_projective_planes`.

    If no construction is possible, then the function raises a ``EmptySetError``
    whereas if no construction is available the function raises a
    ``NotImplementedError``.

    INPUT:

    - ``n`` -- the finite projective plane's order

    EXAMPLES::

        sage: designs.projective_plane(2)
        (7,3,1)-Balanced Incomplete Block Design
        sage: designs.projective_plane(3)
        (13,4,1)-Balanced Incomplete Block Design
        sage: designs.projective_plane(4)
        (21,5,1)-Balanced Incomplete Block Design
        sage: designs.projective_plane(5)
        (31,6,1)-Balanced Incomplete Block Design
        sage: designs.projective_plane(6)
        Traceback (most recent call last):
        ...
        EmptySetError: By the Bruck-Ryser theorem, no projective plane of order 6 exists.
        sage: designs.projective_plane(10)
        Traceback (most recent call last):
        ...
        EmptySetError: No projective plane of order 10 exists by C. Lam, L. Thiel and S. Swiercz "The nonexistence of finite projective planes of order 10" (1989), Canad. J. Math.
        sage: designs.projective_plane(12)
        Traceback (most recent call last):
        ...
        NotImplementedError: If such a projective plane exists, we do not know how to build it.
        sage: designs.projective_plane(14)
        Traceback (most recent call last):
        ...
        EmptySetError: By the Bruck-Ryser theorem, no projective plane of order 14 exists.

    TESTS::

        sage: designs.projective_plane(2197, existence=True)
        True
        sage: designs.projective_plane(6, existence=True)
        False
        sage: designs.projective_plane(10, existence=True)
        False
        sage: designs.projective_plane(12, existence=True)
        Unknown
    """
    from sage.rings.sum_of_squares import is_sum_of_two_squares_pyx

    if n <= 1:
        if existence:
            return False
        raise EmptySetError("There is no projective plane of order <= 1")

    if n == 10:
        if existence:
            return False
        ref = ("C. Lam, L. Thiel and S. Swiercz \"The nonexistence of finite "
               "projective planes of order 10\" (1989), Canad. J. Math.")
        raise EmptySetError("No projective plane of order 10 exists by %s"%ref)

    if (n%4) in [1,2] and not is_sum_of_two_squares_pyx(n):
        if existence:
            return False
        raise EmptySetError("By the Bruck-Ryser theorem, no projective"
                         " plane of order {} exists.".format(n))

    if not is_prime_power(n):
        if existence:
            return Unknown
        raise NotImplementedError("If such a projective plane exists, we do "
                                  "not know how to build it.")

    if existence:
        return True
    else:
        return DesarguesianProjectivePlaneDesign(n, point_coordinates=False, check=check)

def v_4_1_rbibd(v,existence=False):
    r"""
    Return a `(v,4,1)`-RBIBD.

    INPUT:

    - `n` (integer)

    - ``existence`` (boolean; ``False`` by default) -- whether to build the
      design or only answer whether it exists.

    .. SEEALSO::

        - :meth:`IncidenceStructure.is_resolvable`
        - :func:`resolvable_balanced_incomplete_block_design`

    .. NOTE::

        A resolvable `(v,4,1)`-BIBD exists whenever `1\equiv 4\pmod(12)`. This
        function, however, only implements a construction of `(v,4,1)`-BIBD such
        that `v=3q+1\equiv 1\pmod{3}` where `q` is a prime power (see VII.7.5.a
        from [BJL99]_).

    EXAMPLE::

        sage: rBIBD = designs.resolvable_balanced_incomplete_block_design(28,4)
        sage: rBIBD.is_resolvable()
        True
        sage: rBIBD.is_t_design(return_parameters=True)
        (True, (2, 28, 4, 1))

    TESTS::

        sage: for q in prime_powers(2,30):
        ....:     if (3*q+1)%12 == 4:
        ....:         _ = designs.resolvable_balanced_incomplete_block_design(3*q+1,4) # indirect doctest
    """
    # Volume 1, VII.7.5.a from [BJL99]_
    if v%3 != 1 or not is_prime_power((v-1)//3):
        if existence:
            return Unknown
        raise NotImplementedError("I don't know how to build a ({},{},1)-RBIBD!".format(v,4))
    from sage.rings.finite_rings.finite_field_constructor import FiniteField as GF
    q = (v-1)//3
    nn = (q-1)//4
    G = GF(q,'x')
    w = G.primitive_element()
    e = w**(nn)
    assert e**2 == -1

    first_class = [[(w**i,j),(-w**i,j),(e*w**i,j+1),(-e*w**i,j+1)]
                   for i in range(nn) for j in range(3)]

    first_class.append([(0,0),(0,1),(0,2),'inf'])

    label = {p:i for i,p in enumerate(G)}

    classes = [[[v-1 if x=='inf' else (x[1]%3)*q+label[x[0]+g] for x in S]
                for S in first_class]
               for g in G]

    BIBD = BalancedIncompleteBlockDesign(v,
                                         blocks = sum(classes,[]),
                                         k=4,
                                         check=True,
                                         copy=False)
    BIBD._classes = classes
    assert BIBD.is_resolvable()
    return BIBD

def kirkman_triple_system(v,existence=False):
    r"""
    Return a Kirkman Triple System on `v` points.

    A Kirkman Triple System `KTS(v)` is a resolvable Steiner Triple System. It
    exists if and only if `v\equiv 3\pmod{6}`.

    INPUT:

    - `n` (integer)

    - ``existence`` (boolean; ``False`` by default) -- whether to build the
      `KTS(n)` or only answer whether it exists.

    .. SEEALSO::

        :meth:`IncidenceStructure.is_resolvable`

    EXAMPLES:

    A solution to Kirkmman's original problem::

        sage: kts = designs.kirkman_triple_system(15)
        sage: classes = kts.is_resolvable(1)[1]
        sage: names = '0123456789abcde'
        sage: to_name = lambda (r,s,t): ' '+names[r]+names[s]+names[t]+' '
        sage: rows = ['   '.join(('Day {}'.format(i) for i in range(1,8)))]
        sage: rows.extend('   '.join(map(to_name,row)) for row in zip(*classes))
        sage: print '\n'.join(rows)
        Day 1   Day 2   Day 3   Day 4   Day 5   Day 6   Day 7
         07e     18e     29e     3ae     4be     5ce     6de
         139     24a     35b     46c     05d     167     028
         26b     03c     14d     257     368     049     15a
         458     569     06a     01b     12c     23d     347
         acd     7bd     78c     89d     79a     8ab     9bc

    TESTS::

        sage: for i in range(3,300,6):
        ....:     _ = designs.kirkman_triple_system(i)
    """
    if v%6 != 3:
        if existence:
            return False
        raise ValueError("There is no KTS({}) as v!=3 mod(6)".format(v))

    if existence:
        return False

    elif v == 3:
        return BalancedIncompleteBlockDesign(3,[[0,1,2]],k=3,lambd=1)

    elif v == 9:
        classes = [[[0, 1, 5], [2, 6, 7], [3, 4, 8]],
                   [[1, 6, 8], [3, 5, 7], [0, 2, 4]],
                   [[1, 4, 7], [0, 3, 6], [2, 5, 8]],
                   [[4, 5, 6], [0, 7, 8], [1, 2, 3]]]
        KTS = BalancedIncompleteBlockDesign(v,[tr for cl in classes for tr in cl],k=3,lambd=1,copy=False)
        KTS._classes = classes
        return KTS

    # Construction 1.1 from [Stinson91] (originally Theorem 6 from [RCW71])
    #
    # For all prime powers q=1 mod 6, there exists a KTS(2q+1)
    elif ((v-1)//2)%6 == 1 and is_prime_power((v-1)//2):
        from sage.rings.finite_rings.finite_field_constructor import FiniteField as GF
        q = (v-1)//2
        K = GF(q,'x')
        a = K.primitive_element()
        t = (q-1)/6

        # m is the solution of a^m=(a^t+1)/2
        from sage.groups.generic import discrete_log
        m = discrete_log((a**t+1)/2, a)
        assert 2*a**m == a**t+1

        # First parallel class
        first_class = [[(0,1),(0,2),'inf']]
        b0 = K.one(); b1 = a**t; b2 = a**m
        first_class.extend([(b0*a**i,1),(b1*a**i,1),(b2*a**i,2)]
                            for i in range(t)+range(2*t,3*t)+range(4*t,5*t))
        b0 = a**(m+t); b1=a**(m+3*t); b2=a**(m+5*t)
        first_class.extend([[(b0*a**i,2),(b1*a**i,2),(b2*a**i,2)]
                            for i in range(t)])

        # Action of K on the points
        action = lambda v,x : (v+x[0],x[1]) if len(x) == 2 else x

        # relabel to integer
        relabel = {(p,x): i+(x-1)*q
                   for i,p in enumerate(K)
                   for x in [1,2]}
        relabel['inf'] = 2*q

        classes = [[[relabel[action(p,x)] for x in tr] for tr in first_class]
                   for p in K]

        KTS = BalancedIncompleteBlockDesign(v,[tr for cl in classes for tr in cl],k=3,lambd=1,copy=False)

        KTS._classes = classes
#.........这里部分代码省略.........

def skew_hadamard_matrix(n,existence=False, skew_normalize=True, check=True):
    r"""
    Tries to construct a skew Hadamard matrix

    A Hadamard matrix `H` is called skew if `H=S-I`, for `I` the identity matrix
    and `-S=S^\top`. Currently constructions from Section 14.1 of [Ha83]_ and few
    more exotic ones are implemented.

    INPUT:

    - ``n`` (integer) -- dimension of the matrix

    - ``existence`` (boolean) -- whether to build the matrix or merely query if
      a construction is available in Sage. When set to ``True``, the function
      returns:

        - ``True`` -- meaning that Sage knows how to build the matrix

        - ``Unknown`` -- meaning that Sage does not know how to build the
          matrix, but that the design may exist (see :mod:`sage.misc.unknown`).

        - ``False`` -- meaning that the matrix does not exist.

    - ``skew_normalize`` (boolean) -- whether to make the 1st row all-one, and
      adjust the 1st column accordingly. Set to ``True`` by default.

    - ``check`` (boolean) -- whether to check that output is correct before
      returning it. As this is expected to be useless (but we are cautious
      guys), you may want to disable it whenever you want speed. Set to ``True``
      by default.

    EXAMPLES::

        sage: from sage.combinat.matrices.hadamard_matrix import skew_hadamard_matrix
        sage: skew_hadamard_matrix(12).det()
        2985984
        sage: 12^6
        2985984
        sage: skew_hadamard_matrix(1)
        [1]
        sage: skew_hadamard_matrix(2)
        [ 1  1]
        [-1  1]

    TESTS::

        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:
#.........这里部分代码省略.........

def BIBD_from_arc_in_desarguesian_projective_plane(n,k,existence=False):
    r"""
    Returns a `(n,k,1)`-BIBD from a maximal arc in a projective plane.

    This function implements a construction from Denniston [Denniston69]_, who
    describes a maximal :meth:`arc
    <sage.combinat.designs.bibd.BalancedIncompleteBlockDesign.arc>` in a
    :func:`Desarguesian Projective Plane
    <sage.combinat.designs.block_design.DesarguesianProjectivePlaneDesign>` of
    order `2^k`. From two powers of two `n,q` with `n<q`, it produces a
    `((n-1)(q+1)+1,n,1)`-BIBD.

    INPUT:

    - ``n,k`` (integers) -- must be powers of two (among other restrictions).

    - ``existence`` (boolean) -- whether to return the BIBD obtained through
      this construction (default), or to merely indicate with a boolean return
      value whether this method *can* build the requested BIBD.

    EXAMPLES:

    A `(232,8,1)`-BIBD::

        sage: from sage.combinat.designs.bibd import BIBD_from_arc_in_desarguesian_projective_plane
        sage: from sage.combinat.designs.bibd import BalancedIncompleteBlockDesign
        sage: D = BIBD_from_arc_in_desarguesian_projective_plane(232,8)
        sage: BalancedIncompleteBlockDesign(232,D)
        (232,8,1)-Balanced Incomplete Block Design

    A `(120,8,1)`-BIBD::

        sage: D = BIBD_from_arc_in_desarguesian_projective_plane(120,8)
        sage: BalancedIncompleteBlockDesign(120,D)
        (120,8,1)-Balanced Incomplete Block Design

    Other parameters::

        sage: all(BIBD_from_arc_in_desarguesian_projective_plane(n,k,existence=True)
        ....:     for n,k in
        ....:       [(120, 8), (232, 8), (456, 8), (904, 8), (496, 16),
        ....:        (976, 16), (1936, 16), (2016, 32), (4000, 32), (8128, 64)])
        True

    Of course, not all can be built this way::

        sage: BIBD_from_arc_in_desarguesian_projective_plane(7,3,existence=True)
        False
        sage: BIBD_from_arc_in_desarguesian_projective_plane(7,3)
        Traceback (most recent call last):
        ...
        ValueError: This function cannot produce a (7,3,1)-BIBD

    REFERENCE:

    .. [Denniston69] R. H. F. Denniston,
       Some maximal arcs in finite projective planes.
       Journal of Combinatorial Theory 6, no