diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index 04164b1b0b1..19d4f7c95fa 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -2069,10 +2069,15 @@ REFERENCES: *Ten New Orders for Hadamard Matrices of Skew Type*, Publikacije Elektrotehničkog fakulteta. Serija Matematika 2 (1992): 47-59. -.. [Djo1994] \D. Đoković. +.. [Djo1994a] \D. Đoković. *Five New Orders for Hadamard Matrices of Skew Type*, Australasian Journal of Combinatorics 10 (1994): 259-264. +.. [Djo1994b] \D. Đoković. + *Two Hadamard matrices of order 956 of Goethals-Seidel type*, + Combinatorica 14(3) (1994): 375-377. + :doi:`10.1007/BF01212983` + .. [Djo2008a] \D. Đoković. *Skew-Hadamard matrices of orders 188 and 388 exist*, International Mathematical Forum 3 no.22 (2008): 1063-1068. @@ -2083,6 +2088,11 @@ REFERENCES: Journal of Combinatorial Designs 16 (2008): 493-498. :arxiv:`0706.1973` +.. [Djo2008c] \D. Đoković. + *Hadamard matrices of order 764 exist*, + Combinatorica 28(4) (2008): 487-489. + :doi:`10.1007/s00493-008-2384-z` + .. [Djo2023a] \D. Đoković. *Skew-Hadamard matrices of order 276*. :arxiv:`10.48550/ARXIV.2301.02751` @@ -4141,6 +4151,11 @@ REFERENCES: of a genus 2 Jacobian*, Mathematics of Computation 88 (2019), 889-929. :doi:`10.1090/mcom/3358`. +.. [Lon2013] \S. London, + *Constructing New Turyn Type Sequences, T-Sequences and Hadamard Matrices*. + PhD Thesis, University of Illinois at Chicago, 2013. + https://hdl.handle.net/10027/9916 + .. [LOS2012] \C. Lecouvey, M. Okado, M. Shimozono. "Affine crystals, one-dimensional sums and parabolic Lusztig `q`-analogues". Mathematische Zeitschrift. **271** (2012). Issue 3-4. @@ -4516,6 +4531,11 @@ REFERENCES: .. [Mit2008] \A. Mitra. *On the construction of M-sequences via primitive polynomials with a fast identification method*, International Journal of Electronics and Communication Engineering 2(9) (2008): 1991-1996. +.. [Miy1991] \M. Miyamoto. + *A construction of Hadamard matrices*, + Journal of Combinatorial Theory, Series A 57(1) (1991): 86-108. + :doi:`10.1016/0097-3165(91)90008-5` + .. [MKO1998] Hans Munthe--Kaas and Brynjulf Owren. *Computations in a free Lie algebra*. (1998). `Downloadable from Munthe-Kaas's website diff --git a/src/sage/combinat/designs/difference_family.py b/src/sage/combinat/designs/difference_family.py index 609cdf6542e..95c9104e4b3 100644 --- a/src/sage/combinat/designs/difference_family.py +++ b/src/sage/combinat/designs/difference_family.py @@ -1618,8 +1618,8 @@ def is_supplementary_difference_set(Ks, v, lmbda): EXAMPLES:: - sage: from sage.combinat.designs.difference_family import supplementary_difference_set, is_supplementary_difference_set - sage: S1, S2, S3, S4 = supplementary_difference_set(17) + sage: from sage.combinat.designs.difference_family import supplementary_difference_set_from_rel_diff_set, is_supplementary_difference_set + sage: S1, S2, S3, S4 = supplementary_difference_set_from_rel_diff_set(17) sage: is_supplementary_difference_set([S1, S2, S3, S4], 16, 16) True sage: is_supplementary_difference_set([S1, S2, S3, S4], 16, 14) @@ -1629,7 +1629,7 @@ def is_supplementary_difference_set(Ks, v, lmbda): .. SEEALSO:: - :func:`supplementary_difference_set` + :func:`supplementary_difference_set_from_rel_diff_set` """ from sage.groups.additive_abelian.additive_abelian_group import AdditiveAbelianGroup @@ -1651,7 +1651,7 @@ def is_supplementary_difference_set(Ks, v, lmbda): return True -def supplementary_difference_set(q, existence=False, check=True): +def supplementary_difference_set_from_rel_diff_set(q, existence=False, check=True): r"""Construct `4-\{2v; v, v+1, v, v; 2v\}` supplementary difference sets where `q=2v+1`. The sets are created from relative difference sets as detailed in Theorem 3.3 of [Spe1975]_. this construction @@ -1682,8 +1682,8 @@ def supplementary_difference_set(q, existence=False, check=True): EXAMPLES:: - sage: from sage.combinat.designs.difference_family import supplementary_difference_set - sage: supplementary_difference_set(17) #random + sage: from sage.combinat.designs.difference_family import supplementary_difference_set_from_rel_diff_set + sage: supplementary_difference_set_from_rel_diff_set(17) #random ([0, 2, 5, 6, 8, 10, 13, 14], [0, 1, 2, 6, 7, 9, 10, 14, 15], [0, 1, 2, 6, 11, 12, 13, 15], @@ -1691,31 +1691,31 @@ def supplementary_difference_set(q, existence=False, check=True): If existence is ``True``, the function returns a boolean:: - sage: supplementary_difference_set(7, existence=True) + sage: supplementary_difference_set_from_rel_diff_set(7, existence=True) False - sage: supplementary_difference_set(17, existence=True) + sage: supplementary_difference_set_from_rel_diff_set(17, existence=True) True TESTS:: sage: from sage.combinat.designs.difference_family import is_supplementary_difference_set - sage: is_supplementary_difference_set(supplementary_difference_set(17), 16, 16) + sage: is_supplementary_difference_set(supplementary_difference_set_from_rel_diff_set(17), 16, 16) True - sage: is_supplementary_difference_set(supplementary_difference_set(9), 8, 8) + sage: is_supplementary_difference_set(supplementary_difference_set_from_rel_diff_set(9), 8, 8) True - sage: supplementary_difference_set(7) + sage: supplementary_difference_set_from_rel_diff_set(7) Traceback (most recent call last): ... ValueError: There is no s for which m-1 is an odd prime power - sage: supplementary_difference_set(8) + sage: supplementary_difference_set_from_rel_diff_set(8) Traceback (most recent call last): ... ValueError: q must be an odd prime power - sage: supplementary_difference_set(8, existence=True) + sage: supplementary_difference_set_from_rel_diff_set(8, existence=True) False - sage: supplementary_difference_set(7, existence=True) + sage: supplementary_difference_set_from_rel_diff_set(7, existence=True) False - sage: supplementary_difference_set(1, existence=True) + sage: supplementary_difference_set_from_rel_diff_set(1, existence=True) False .. SEEALSO:: @@ -1823,7 +1823,7 @@ def get_fixed_relative_difference_set(rel_diff_set, as_elements=False): `\{td | d\in R\}= R`. In addition, the set returned by this function will contain the element `0`. This is needed in the - construction of supplementary difference sets (see :func:`supplementary_difference_set`). + construction of supplementary difference sets (see :func:`supplementary_difference_set_from_rel_diff_set`). INPUT: @@ -1929,12 +1929,12 @@ def is_fixed_relative_difference_set(R, q): def skew_supplementary_difference_set(n, existence=False, check=True): r"""Construct `4-\{n; n_1, n_2, n_3, n_4; \lambda\}` supplementary difference sets where `S_1` is skew and `n_1 + n_2 + n_3 + n_4 = n+\lambda`. - These sets are constructed from available data, as described in [Djo1994]_. The set `S_1 \subset G` is + These sets are constructed from available data, as described in [Djo1994a]_. The set `S_1 \subset G` is always skew, i.e. `S_1 \cap (-S_1) = \emptyset` and `S_1 \cup (-S_1) = G \setminus \{0\}`. The data is taken from: - * `n = 103, 151`: [Djo1994]_ + * `n = 103, 151`: [Djo1994a]_ * `n = 67, 113, 127, 157, 163, 181, 241`: [Djo1992a]_ * `n = 37, 43`: [Djo1992b]_ * `n = 39, 49, 65, 93, 121, 129, 133, 217, 219, 267`: [Djo1992c]_ @@ -2155,6 +2155,88 @@ def skew_supplementary_difference_set(n, existence=False, check=True): 267: [1, 4, 16, 64, 67, 91, 97, 121, 217, 223, 256], } + if existence: + return n in indices + + if n not in indices: + raise ValueError(f'Skew SDS of order {n} not yet implemented.') + + S1, S2, S3, S4 = _construction_supplementary_difference_set(n, H_db[n], indices[n], cosets_gens[n], check=False) + + if check: + lmbda = len(S1) + len(S2) + len(S3) + len(S4) - n + assert is_supplementary_difference_set([S1, S2, S3, S4], n, lmbda) + assert _is_skew_set(S1, n) + + return S1, S2, S3, S4 + + +def _construction_supplementary_difference_set(n, H, indices, cosets_gen, check=True): + r"""Construct `4-\{n; n_1, n_2, n_3, n_4; \lambda\}` supplementary difference sets where `n_1 + n_2 + n_3 + n_4 = n+\lambda`. + + This construction is described in [Djo1994a]_. + + Let H be a subgroup of Zmod(n) of order `t`. We construct the `2s = n/t` cosets + `\alpha_0, .., \alpha_{2s-1}` by using the elements contained in a sequence `A`: + `\alpha_{2i} = a_iH` and `\alpha_{2i+1} = -\alpha_{2i}`. + + Then, we use four indices sets `J_1, J_2, J_3, J_4` to construct the four + supplementary difference sets: `S_i = \bigcup_{j\in J__i} \alpha_i`. Note that + if `J_i` contains the value `-1`, this function will add `0` to the set `S_i`. + + To construct a coset `\alpha_{2i}` by listing it directly, replace the `2i`-th + element of the list `A` with the desired set. + + INPUT: + + - ``n`` -- integer, the parameter of the supplementary difference set. + + - ``H`` -- list of integers, the set `H` used to construct the cosets. + + - ``indices`` -- list containing four list of integers, which are the sets + `J_1, J_2, J_3, J_4` described above. + + - ``cosets_gen`` -- list containing integers or list of integers, is the set `A` + described above. + + - ``check`` -- boolean (default True). If true, check that the sets + are supplementary difference sets. Setting this parameter to False may speed + up the computation considerably. + + OUTPUT: + + The function returns the 4 sets (containing integers modulo `n`). + + TESTS:: + + sage: from sage.combinat.designs.difference_family import is_supplementary_difference_set, _construction_supplementary_difference_set + sage: H = [1, 10, -11] + sage: cosets_gen = [1, 2, 3, 5, 6, 9] + sage: indices = [[0, 3, 5, 7, 9, 10], [0, 5, 6, 7, 8], [1, 2, 6, 7, 9], [2, 6, 8, 9, 10]] + sage: _construction_supplementary_difference_set(37, H, indices, cosets_gen) + ([1, 10, 26, 35, 17, 22, 34, 7, 33, 32, 24, 18, 31, 14, 29, 9, 16, 12], + [1, 10, 26, 34, 7, 33, 5, 13, 19, 32, 24, 18, 6, 23, 8], + [36, 27, 11, 2, 20, 15, 5, 13, 19, 32, 24, 18, 31, 14, 29], + [2, 20, 15, 5, 13, 19, 6, 23, 8, 31, 14, 29, 9, 16, 12]) + sage: H = [1, 16, 22] + sage: cosets_gen = [1, 2, 3, 4, 6, 8, [13]] + sage: indices = [[1, 3, 5, 6, 8, 10, 12], [0, 1, 5, 8, 12, 13], [1, 3, 4, 7, 9, 12, 13], [0, 1, 2, 3, 7, 8]] + sage: _construction_supplementary_difference_set(39, H, indices, cosets_gen) + ([38, 23, 17, 37, 7, 34, 36, 30, 12, 4, 25, 10, 6, 18, 15, 8, 11, 20, 13], + [1, 16, 22, 38, 23, 17, 36, 30, 12, 6, 18, 15, 13, 26], + [38, 23, 17, 37, 7, 34, 3, 9, 27, 35, 14, 29, 33, 21, 24, 13, 26], + [1, 16, 22, 38, 23, 17, 2, 32, 5, 37, 7, 34, 35, 14, 29, 6, 18, 15]) + sage: H = [1, 4, 11, 16, 21, -2, -8] + sage: cosets_gen = [1, 3, 7] + sage: indices = [[1, 2, 4], [1, 2, 4], [0, 2, 3], [3, 4, -1]] + sage: sets = _construction_supplementary_difference_set(43, H, indices, cosets_gen, check=False) + sage: is_supplementary_difference_set(sets, 43, 35) + True + + .. SEEALSO:: + + :func:`skew_supplementary_difference_set` + """ def generate_set(index_set, cosets): S = [] for idx in index_set: @@ -2164,17 +2246,11 @@ def generate_set(index_set, cosets): S += cosets[idx] return S - if existence: - return n in indices - - if n not in indices: - raise ValueError(f'Skew SDS of order {n} not yet implemented.') - Z = Zmod(n) - H = list(map(Z, H_db[n])) + H = list(map(Z, H)) cosets = [] - for el in cosets_gens[n]: + for el in cosets_gen: if isinstance(el, list): even_coset = [Z(x) for x in el] else: @@ -2183,22 +2259,118 @@ def generate_set(index_set, cosets): cosets.append(even_coset) cosets.append(odd_coset) - S1 = generate_set(indices[n][0], cosets) - S2 = generate_set(indices[n][1], cosets) - S3 = generate_set(indices[n][2], cosets) - S4 = generate_set(indices[n][3], cosets) + S1 = generate_set(indices[0], cosets) + S2 = generate_set(indices[1], cosets) + S3 = generate_set(indices[2], cosets) + S4 = generate_set(indices[3], cosets) if check: lmbda = len(S1) + len(S2) + len(S3) + len(S4) - n assert is_supplementary_difference_set([S1, S2, S3, S4], n, lmbda) - assert _is_skew_set(S1, n) return S1, S2, S3, S4 + +def supplementary_difference_set(n, existence=False, check=True): + r"""Construct `4-\{n; n_1, n_2, n_3, n_4; \lambda\}` supplementary difference sets where `n_1 + n_2 + n_3 + n_4 = n+\lambda`. + + These sets are constructed from available data, as described in [Djo1994a]_. + + The data for `n=191` is taken from [Djo2008c]_, and data for `n=239` is from [Djo1994b]_. + Additional SDS are constructed using :func:`skew_supplementary_difference_set`. + + INPUT: + + - ``n`` -- integer, the parameter of the supplementary difference set. + + - ``existence`` -- boolean (default False). If true, only check whether the + supplementary difference sets can be constructed. + + - ``check`` -- boolean (default True). If true, check that the sets are + supplementary difference sets before returning them. Setting this parameter + to False may speed up the computation considerably. + + OUTPUT: + + If ``existence`` is false, the function returns the 4 sets (containing integers + modulo `n`), or raises an error if data for the given ``n`` is not available. + If ``existence`` is true, the function returns a boolean representing whether + skew supplementary difference sets can be constructed. + + EXAMPLES:: + + sage: from sage.combinat.designs.difference_family import supplementary_difference_set + sage: S1, S2, S3, S4 = supplementary_difference_set(191) # long time + + If existence is ``True``, the function returns a boolean :: + + sage: supplementary_difference_set(191, existence=True) + True + sage: supplementary_difference_set(17, existence=True) + False + + TESTS:: + + sage: from sage.combinat.designs.difference_family import is_supplementary_difference_set + sage: S1, S2, S3, S4 = supplementary_difference_set(191, check=False) + sage: is_supplementary_difference_set([S1, S2, S3, S4], 191, len(S1)+len(S2)+len(S3)+len(S4)-191) # long time + True + sage: S1, S2, S3, S4 = supplementary_difference_set(37, check=False) + sage: supplementary_difference_set(7) + Traceback (most recent call last): + ... + ValueError: SDS of order 7 not yet implemented. + sage: supplementary_difference_set(7, existence=True) + False + sage: supplementary_difference_set(127, existence=True) + True + """ + + indices = { + 191: [[1, 7, 9, 10, 11, 13, 17, 18, 25, 26, 30, 31, 33, 34, 35, 36, 37], + [1, 4, 7, 9, 11, 12, 13, 14, 19, 21, 22, 23, 24, 25, 26, 29, 36, 37], + [0, 3, 4, 5, 7, 8, 9, 16, 17, 19, 24, 25, 29, 30, 31, 33, 35, 37], + [1, 3, 4, 5, 8, 11, 14, 18, 19, 20, 21, 23, 24, 25, 28, 29, 30, 32, 34, 35]], + 239: [[0, 1, 2, 3, 4, 5, 6, 7, 14, 18, 19, 21, 24, 25, 29, 30], + [0, 1, 3, 7, 9, 12, 15, 18, 20, 22, 26, 28, 29, 30, 31, 32, 33], + [2, 3, 4, 5, 8, 9, 10, 11, 13, 17, 19, 21, 22, 24, 27, 31, 32], + [0, 1, 2, 3, 6, 7, 8, 11, 13, 15, 17, 18, 19, 22, 25, 26, 27, 32, 33]], + } + + cosets_gens = { + 191: [1, 2, 3, 4, 6, 8, 9, 11, 12, 13, 16, 17, 18, 19, 22, 32, 36, 38, 41], + 239: [1, 2, 3, 4, 5, 6, 7, 8, 9, 13, 14, 15, 18, 21, 28, 35, 42], + } + + H_db = { + 191: [1, 39, 184, 109, 49], + 239: [1, 10, 24, 44, 98, 100, 201], + } + + if existence: + return n in indices or skew_supplementary_difference_set(n, existence=True) + + sets = None + if n in indices: + sets = _construction_supplementary_difference_set(n, H_db[n], indices[n], cosets_gens[n], check=False) + elif skew_supplementary_difference_set(n, existence=True): + sets = skew_supplementary_difference_set(n, check=False) + + if sets is None: + raise ValueError(f'SDS of order {n} not yet implemented.') + + S1, S2, S3, S4 = sets + if check: + lmbda = len(S1) + len(S2) + len(S3) + len(S4) - n + assert is_supplementary_difference_set([S1, S2, S3, S4], n, lmbda) + + return S1, S2, S3, S4 + + def _is_skew_set(S, n): r"""Check if `S` is a skew set over the set of integers modulo `n`. - From [Djo1994]_, a set `S \subset G` (where `G` is a finite abelian group of order `n`) is of skew + From [Djo1994a]_, a set `S \subset G` (where `G` is a finite abelian group of order `n`) is of skew type if `S_1 \cap (-S_1) = \emptyset` and `S_1 \cup (-S_1) = G\setminus \{0\}`. INPUT: diff --git a/src/sage/combinat/matrices/hadamard_matrix.py b/src/sage/combinat/matrices/hadamard_matrix.py index bf915399d20..0434d20d8c2 100644 --- a/src/sage/combinat/matrices/hadamard_matrix.py +++ b/src/sage/combinat/matrices/hadamard_matrix.py @@ -179,6 +179,56 @@ def hadamard_matrix_paleyI(n, normalize=True): return H +def symmetric_conference_matrix_paley(n): + r""" + Construct a symmetric conference matrix of order n. + + A conference matrix is an `n\times n` matrix `C` with 0s on the main diagonal + and 1s and -1s elsewhere, satisfying `CC^\top=(n-1)I`. This construction assumes + that `q = n-1` is a prime power, with `q \cong 1 \mod 4`. See [Hora]_ or [Lon2013]_. + + These matrices are used in :func:`hadamard_matrix_paleyII`. + + INPUT: + + - ``n`` -- integer, the order of the symmetric conference matrix to consruct. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import symmetric_conference_matrix_paley + sage: symmetric_conference_matrix_paley(6) + [ 0 1 1 1 1 1] + [ 1 0 1 -1 -1 1] + [ 1 1 0 1 -1 -1] + [ 1 -1 1 0 1 -1] + [ 1 -1 -1 1 0 1] + [ 1 1 -1 -1 1 0] + + TESTS:: + + sage: symmetric_conference_matrix_paley(5) + Traceback (most recent call last): + ... + ValueError: The order 5 is not covered by Paley construction of symmetric conference matrices. + """ + q = n - 1 + if not (is_prime_power(q) and (q % 4 == 1)): + raise ValueError("The order %s is not covered by Paley construction of symmetric conference matrices." % n) + + from sage.rings.finite_rings.finite_field_constructor import FiniteField + K = FiniteField(q, '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[0, i] = 1 + H[i, 0] = 1 + H[i, i] = 0 + return H + + def hadamard_matrix_paleyII(n): r""" Implement the Paley type II construction. @@ -226,17 +276,7 @@ def hadamard_matrix_paleyII(n): if not (n % 2 == 0 and is_prime_power(q) and (q % 4 == 1)): raise ValueError("The order %s is not covered by the Paley type II construction." % n) - from sage.rings.finite_rings.finite_field_constructor import FiniteField - K = FiniteField(q, '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(q+1): - H[0, i] = 1 - H[i, 0] = 1 - H[i, i] = 0 + H = symmetric_conference_matrix_paley(q+1) tr = { 0: matrix(2, 2, [ 1, -1, -1, -1]), 1: matrix(2, 2, [ 1, 1, 1, -1]), @@ -247,6 +287,114 @@ def hadamard_matrix_paleyII(n): return normalise_hadamard(H) +def hadamard_matrix_miyamoto_construction(n, existence=False, check=True): + r""" + Construct Hadamard matrix using Miyamoto construction. + + If `q = n/4` is a prime power, and there exists an Hadamard matrix of order + `q-1`, then a Hadamard matrix of order `n` can be constructed (see [Miy1991]_). + + INPUT: + + - ``n`` -- integer, the order of the matrix to be constructed. + + - ``check`` -- boolean: if True (default), check the the matrix is a Hadamard + before returning. + + - ``existence`` -- boolean (default False): if True, only check if the matrix exists. + + OUTPUT: + + If ``existence`` is false, returns the Hadamard matrix of order `n`. It raises + an error if no data is available to construct the matrix of the given order, + or if `n` does not satisfies the constraints. + If ``existence`` is true, returns a boolean representing whether the matrix + can be constructed or not. + + EXAMPLES: + + By default the function returns the Hadamard matrix :: + + sage: from sage.combinat.matrices.hadamard_matrix import hadamard_matrix_miyamoto_construction + sage: hadamard_matrix_miyamoto_construction(20) + 20 x 20 dense matrix over Integer Ring... + + If ``existence`` is set to True, the function returns a boolean :: + + sage: hadamard_matrix_miyamoto_construction(36, existence=True) + True + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix + sage: is_hadamard_matrix(hadamard_matrix_miyamoto_construction(68, check=False)) + True + sage: hadamard_matrix_miyamoto_construction(64, existence=True) + False + sage: hadamard_matrix_miyamoto_construction(64) + Traceback (most recent call last): + ... + ValueError: The order 64 is not covered by Miyamoto construction. + sage: hadamard_matrix_miyamoto_construction(14) + Traceback (most recent call last): + ... + ValueError: No Hadamard matrix of order 14 exists. + """ + if n < 0 or n % 4 != 0: + raise ValueError(f'No Hadamard matrix of order {n} exists.') + + q = n // 4 + if existence: + return is_prime_power(q) and q % 4 == 1 and hadamard_matrix(q-1, existence=True) is True + + if not (is_prime_power(q) and q % 4 == 1 and hadamard_matrix(q-1, existence=True)): + raise ValueError(f'The order {n} is not covered by Miyamoto construction.') + + m = (q-1) // 2 + + C = symmetric_conference_matrix_paley(q + 1) + + neg = [i for i in range(2, m+2) if C[1, i] == -1] + pos = [i for i in range(m+2, 2*m+2) if C[1, i] == 1] + + for i, j in zip(neg, pos): + C.swap_rows(i, j) + C.swap_columns(i, j) + + C1 = -C.submatrix(row=2, col=2, nrows=m, ncols=m) + C2 = C.submatrix(row=2, col=m+2, nrows=m, ncols=m) + C4 = C.submatrix(row=m+2, col=m+2, nrows=m, ncols=m) + + K = hadamard_matrix(q - 1) + K1 = K.submatrix(row=0, col=0, nrows=(q-1)//2, ncols=(q-1)//2) + K2 = K.submatrix(row=0, col=(q-1)//2, nrows=(q-1)//2, ncols=(q-1)//2) + K3 = -K.submatrix(row=(q-1)//2, col=0, nrows=(q-1)//2, ncols=(q-1)//2) + K4 = K.submatrix(row=(q-1)//2, col=(q-1)//2, nrows=(q-1)//2, ncols=(q-1)//2) + + Zr = zero_matrix(m) + Us = [[C1, C2, Zr, Zr], [C2.T, C4, Zr, Zr], [Zr, Zr, C1, C2], [Zr, Zr, C2.T, C4]] + Vs = [[I(m), Zr, K1, K2], [Zr, I(m), K3, K4], [K1.T, K3.T, I(m), Zr], [K2.T, K4.T, Zr, I(m)]] + + def T(i, j): + return block_matrix([[Us[i][j]+Vs[i][j], Us[i][j]-Vs[i][j]], + [Us[i][j]-Vs[i][j], Us[i][j]+Vs[i][j]]]) + + e = matrix([[1] * (2*m)]) + one = matrix([1]) + H = block_matrix([[ one, -e, one, e, one, e, one, e], + [-e.T, T(0, 0), e.T, T(0, 1), e.T, T(0, 2), e.T, T(0, 3)], + [-one, -e, one, -e, one, e, -one, -e], + [-e.T, -T(1, 0), -e.T, T(1, 1), e.T, T(1, 2), -e.T, -T(1, 3)], + [-one, -e, -one, -e, one, -e, one, e], + [-e.T, -T(2, 0), -e.T, -T(2, 1), -e.T, T(2, 2), e.T, T(2, 3)], + [-one, -e, one, e, -one, -e, one, -e], + [-e.T, -T(3, 0), e.T, T(3, 1), -e.T, -T(3, 2), -e.T, T(3, 3)]]) + + if check: + assert is_hadamard_matrix(H) + return H + + def hadamard_matrix_williamson_type(a, b, c, d, check=True): r""" Construction of Williamson type Hadamard matrix. @@ -319,7 +467,8 @@ def williamson_type_quadruples_smallcases(n, existence=False): Williamson construction of Hadamard matrices. Namely, the function returns the first row of 4 `n\times n` circulant matrices with the properties described in :func:`sage.combinat.matrices.hadamard_matrix.hadamard_matrix_williamson_type`. - The matrices for n=29 and n=43 are given in [Ha83]_. + The matrices for `n = 3, 5, ..., 29, 37, 43` are given in [Ha83]_. The matrices + for `n = 31, 33, 39, 41, 45, 49, 51, 55, 57, 61, 63` are given in [Lon2013]_. INPUT: @@ -354,31 +503,67 @@ def williamson_type_quadruples_smallcases(n, existence=False): ValueError: The Williamson type quadruple of order 123 is not yet implemented. """ db = { - 1: ([1], [1], [1], [1]), - 7: ([1, -1, -1, 1, 1, -1, -1], - [1, -1, 1, -1, -1, 1, -1], - [1, 1, -1, -1, -1, -1, 1], - [1, -1, -1, -1, -1, -1, -1]), - 9: ([1, -1, -1, -1, 1, 1, -1, -1, -1], - [1, -1, -1, 1, -1, -1, 1, -1, -1], - [1, -1, 1, -1, -1, -1, -1, 1, -1], - [1, 1, -1, -1, -1, -1, -1, -1, 1]), - 29: ([1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1], - [1, -1, 1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1], - [1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1], - [1, 1, -1, -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, 1]), - 43: ([1, 1, -1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1], - [1, 1, 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, 1, 1], - [1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, 1], - [1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1]), + 1: ('+', '+', '+', '+'), + 3: ('+++', '+--', '+--', '+--'), + 5: ('+-++-', '++--+', '+----', '+----'), + 7: ('+--++--', '+-+--+-', '++----+', '+------'), + 9: ('+---++---', '+--+--+--', '+-+----+-', '++------+'), + 11: ('++--------+', '++-+-++-+-+', '++-++--++-+', '+-++----++-'), + 13: ('++++-+--+-+++', '+---+-++-+---', '++---+--+---+', '++---+--+---+'), + 15: ('+-+---++++---+-', '++-++------++-+', + '++-++++--++++-+', '++-++-+--+-++-+'), + 17: ('+---+++----+++---', '++-+---+--+---+-+', + '+--+-++++++++-+--', '+-++-+++--+++-++-'), + 19: ('++--+++-+--+-+++--+', '++-++--+-++-+--++-+', + '+-+---++++++++---+-', '++--+-++++++++-+--+'), + 21: ('+--++++---++---++++--', '++++-+---+--+---+-+++', + '++--+-+-++--++-+-+--+', '++-+++++-+--+-+++++-+'), + 23: ('++---+---+-++-+---+---+', '+-++-++--++++++--++-++-', + '+++---++-+-++-+-++---++', '+++-+++-+------+-+++-++'), + 25: ('++++-+-+-+--++--+-+-+-+++', '++--+--+-++++++++-+--+--+', + '+++--+--++++--++++--+--++', '+-+--+++--++++++--+++--+-'), + 27: ('+--+--+-+++--++--+++-+--+--', '+++-++-+---++--++---+-++-++', + '+---+++++-+-++++-+-+++++---', '+---+++++-+-++++-+-+++++---'), + 29: ('+++---++--+-+----+-+--++---++', '+-+---++--+-++++++-+--++---+-', + '++++-++-+---++++++---+-++-+++', '++--+--+-+++-++++-+++-+--+--+'), + 31: ('++++++-+--+---++++---+--+-+++++', '+--++---+-+-++----++-+-+---++--', + '+--++---+-+-++----++-+-+---++--', '+-----+-++-+++----+++-++-+-----'), + 33: ('++++++-+-+-+++------+++-+-+-+++++', '++-+-++-+----+++--+++----+-++-+-+', + '++--++-+++-+--+-++-+--+-+++-++--+', '+--++--+++++-++----++-+++++--++--'), + 37: ('+--+-+-+-++---+--++++--+---++-+-+-+--', '+---++-++--+-+-++----++-+-+--++-++---', + '+++++-+-----++----++----++-----+-++++', '+--+++-+-----+----++----+-----+-+++--'), + 39: ('+++--+-+-----+--++----++--+-----+-+--++', '+++--++-+---+-+--+----+--+-+---+-++--++', + '++++---+--++----+-+--+-+----++--+---+++', '+---++-+-+-----+++-++-+++-----+-+-++---'), + 41: ('++++--+-++++-++--++----++--++-++++-+--+++', '++++--+-++++-++--++----++--++-++++-+--+++', + '+++-++-+-+-+-----+++--+++-----+-+-+-++-++', '+--+--+-+-+-+++++---++---+++++-+-+-+--+--'), + 43: ('++---++++-+--+--++--------++--+--+-++++---+', '+++-+-++--+-+-++++-+----+-++++-+-+--++-+-++', + '++-++++++----+-+--++-++-++--+-+----++++++-+', '+---++--++++-+-+++-++--++-+++-+-++++--++---'), + 45: ('+++++-++----+-++--++-++-++--++-+----++-++++', '+++---++--+-+-+-++--------++-+-+-+--++---++', + '++-+-++++-+--+--+++--++--+++--+--+-++++-+-+', '+-++-----++++-+-+++-++++-+++-+-++++-----++-'), + 49: ('++++-++-+---++-+++---++-++-++---+++-++---+-++-+++', '++++-++-+---++-+++---++-++-++---+++-++---+-++-+++', + '+----+-++++--+-+++-+-+++--+++-+-+++-+--++++-+----', '+++++-+----++-+---+-+---++---+-+---+-++----+-++++'), + 51: ('+---+++-++-+-+++--+++++--++--+++++--+++-+-++-+++---', '----+++-++-+-+++--+++++--++--+++++--+++-+-++-+++---', + '-+--+----+-+++-+-+++++--+--+--+++++-+-+++-+----+--+', '-+--+----+-+++-+-+++++--+--+--+++++-+-+++-+----+--+'), + 55: ('+-+--+-+-++--+-+++++-+++--++++--+++-+++++-+--++-+-+--+-', '--+--+-+-++--+-+++++-+++--++++--+++-+++++-+--++-+-+--+-', + '+++----++-++--++----+-+-++++++++-+-+----++--++-++----++', '+++----++-++--++----+-+-++++++++-+-+----++--++-++----++'), + 57: ('+---++-+--++++-+++-++---+-++++++-+---++-+++-++++--+-++---', '----++-+--++++-+++-++---+-++++++-+---++-+++-++++--+-++---', + '--+-+-+++--+--+-++---+++++-++++-+++++---++-+--+--+++-+-+-', '--+-+-+++--+--+-++---+++++-++++-+++++---++-+--+--+++-+-+-'), + 61: ('++--+--++--+-+-++++--+-----+------+-----+--++++-+-+--++--+--+', '++--+--++--+-+-++++--+-----+------+-----+--++++-+-+--++--+--+', + '+---+-+-++++---++--+-++-+---++++++---+-++-+--++---++++-+-+---', '++++-+-+----+++--++-+--+-+++------+++-+--+-++--+++----+-+-+++'), + 63: ('++-+++--++-++--+--+-++-+-+++--------+++-+-++-+--+--++-++--+++-+', '-+-+++--++-++--+--+-++-+-+++--------+++-+-++-+--+--++-++--+++-+', + '++++-++-+-++++-+---+---+++---++++++---+++---+---+-++++-+-++-+++', '++++-++-+-++++-+---+---+++---++++++---+++---+---+-++++-+-++-+++'), } + def pmtoZ(s): + return [1 if x == '+' else -1 for x in s] + if existence: return n in db if n not in db: raise ValueError("The Williamson type quadruple of order %s is not yet implemented." % n) - a, b, c, d = map(vector, db[n]) + + a, b, c, d = map(lambda s: vector(pmtoZ(s)), db[n]) return a, b, c, d @@ -406,10 +591,10 @@ def williamson_hadamard_matrix_smallcases(n, existence=False, check=True): 116 x 116 dense matrix over Integer Ring... sage: williamson_hadamard_matrix_smallcases(172) 172 x 172 dense matrix over Integer Ring... - sage: williamson_hadamard_matrix_smallcases(100) + sage: williamson_hadamard_matrix_smallcases(1000) Traceback (most recent call last): ... - ValueError: The Williamson type Hadamard matrix of order 100 is not yet implemented. + ValueError: The Williamson type Hadamard matrix of order 1000 is not yet implemented. """ assert n % 4 == 0 @@ -738,6 +923,92 @@ def _construction_goethals_seidel_matrix(A, B, C, D): [-D*R, -C.T*R, B.T*R, A]]) +def hadamard_matrix_from_sds(n, existence=False, check=True): + r"""Construction of Hadamard matrices from supplementary difference sets. + + Given four SDS with parameters `4-\{n; n_1, n_2, n_3, n_4; \lambda\}` with + `n_1 + n_2 + n_3 + n_4 = n+\lambda` we can construct four (-1,+1) sequences `a_i = (a_{i,0},...,a_{i,n-1})` + where `a_{i,j} = -1` iff `j \in S_i`. These will be the fist rows of four circulant + matrices `A_1, A_2, A_3, A_4` which, when plugged into the Goethals-Seidel array, create an + Hadamard matrix of order `4n` (see [Djo1994b]_). + + The supplementary difference sets are taken from + :func:`sage.combinat.designs.difference_family.supplementary_difference_set`. + + INPUT: + + - ``n`` -- integer, the order of the matrix to be constructed. + + - ``check`` -- boolean: if True (default), check the the matrix is a Hadamard + before returning. + + - ``existence`` -- boolean (default False): if True, only check if the matrix exists. + + OUTPUT: + + If ``existence`` is false, returns the Hadamard matrix of order `n`. It raises + an error if no data is available to construct the matrix of the given order, + or if `n` is not a multiple of `4`. + If ``existence`` is true, returns a boolean representing whether the matrix + can be constructed or not. + + EXAMPLES: + + By default The function returns the Hadamard matrix :: + + sage: from sage.combinat.matrices.hadamard_matrix import hadamard_matrix_from_sds + sage: hadamard_matrix_from_sds(148) + 148 x 148 dense matrix over Integer Ring... + + If ``existence`` is set to True, the function returns a boolean :: + + sage: hadamard_matrix_from_sds(764, existence=True) + True + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import hadamard_matrix_from_sds, is_hadamard_matrix + sage: is_hadamard_matrix(hadamard_matrix_from_sds(172)) + True + sage: hadamard_matrix_from_sds(64, existence=True) + False + sage: hadamard_matrix_from_sds(64) + Traceback (most recent call last): + ... + ValueError: SDS of order 16 not yet implemented. + sage: hadamard_matrix_from_sds(14) + Traceback (most recent call last): + ... + ValueError: n must be a positive multiple of four. + """ + from sage.combinat.designs.difference_family import supplementary_difference_set + + if n <= 0 or n % 4 != 0: + raise ValueError(f'n must be a positive multiple of four.') + t = n // 4 + + if existence: + return supplementary_difference_set(t, existence=True) + + S1, S2, S3, S4 = supplementary_difference_set(t, check=False) + a = [-1 if i in S1 else 1 for i in range(t)] + b = [-1 if i in S2 else 1 for i in range(t)] + c = [-1 if i in S3 else 1 for i in range(t)] + d = [-1 if i in S4 else 1 for i in range(t)] + + if n == 956: + a, b, c, d = [-el for el in d], a, b, c + + 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*t*I(t) + + H = _construction_goethals_seidel_matrix(A, B, C, D) + if check: + assert is_hadamard_matrix(H) + return H + + def hadamard_matrix_cooper_wallis_construction(x1, x2, x3, x4, A, B, C, D, check=True): r""" Create an Hadamard matrix using the contruction detailed in [CW1972]_. @@ -1133,7 +1404,7 @@ def hadamard_matrix_spence_construction(n, existence=False, check=True): r"""Create an Hadamard matrix of order `n` using Spence construction. This construction (detailed in [Spe1975]_), uses supplementary difference sets implemented in - :func:`sage.combinat.designs.difference_family.supplementary_difference_set` to create the + :func:`sage.combinat.designs.difference_family.supplementary_difference_set_from_rel_diff_set` to create the desired matrix. INPUT: @@ -1180,19 +1451,19 @@ def hadamard_matrix_spence_construction(n, existence=False, check=True): ... AssertionError """ - from sage.combinat.designs.difference_family import supplementary_difference_set + from sage.combinat.designs.difference_family import supplementary_difference_set_from_rel_diff_set assert n % 4 == 0 and n > 0 q = n//4 if existence: - return supplementary_difference_set(q, existence=True) + return supplementary_difference_set_from_rel_diff_set(q, existence=True) - if not supplementary_difference_set(q, existence=True): + if not supplementary_difference_set_from_rel_diff_set(q, existence=True): raise ValueError(f'The order {n} is not covered by Spence construction.') - S1, S2, S3, S4 = supplementary_difference_set(q, check=False) + S1, S2, S3, S4 = supplementary_difference_set_from_rel_diff_set(q, check=False) A1 = matrix.circulant([1 if j in S1 else -1 for j in range(q-1)]) A2 = matrix.circulant([1 if j in S4 else -1 for j in range(q-1)]) @@ -1501,6 +1772,14 @@ def hadamard_matrix(n, existence=False, check=True): if existence: return True M = turyn_type_hadamard_matrix_smallcases(n, check=False) + elif hadamard_matrix_miyamoto_construction(n, existence=True): + if existence: + return True + M = hadamard_matrix_miyamoto_construction(n, check=False) + elif hadamard_matrix_from_sds(n, existence=True): + if existence: + return True + M = hadamard_matrix_from_sds(n, check=False) elif hadamard_matrix_spence_construction(n, existence=True): if existence: return True @@ -2261,7 +2540,7 @@ def skew_hadamard_matrix_324(): r""" Construct a skew Hadamard matrix of order 324. - The construction is taken from [Djo1994]_. It uses four supplementary difference sets `S_1, S_2, S_3, S_4`, + The construction is taken from [Djo1994a]_. It uses four supplementary difference sets `S_1, S_2, S_3, S_4`, with `S_1` of skew type. These are then used to generate four matrices of order `81`, which are inserted into the Goethals-Seidel array.