diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index 7bef2a9a4e5..9fbd670ef14 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -624,7 +624,7 @@ REFERENCES: and variants*. :arxiv:`1810.00789` -.. [BDKR2013] \D. Best, D.Z. Dokovic, H. Kharaghani and H. Ramp. +.. [BDKR2013] \D. Best, D.Ž. Đoković, H. Kharaghani and H. Ramp. *Turyn-Type Sequences: Classification, Enumeration, and Construction*, Journal of Combinatorial Designs 21(1) (2013): 24-35. :doi:`10.1002/jcd.21318` @@ -1216,6 +1216,11 @@ REFERENCES: Proc. London Math. Soc. (3) **77** (1998), 551–581. :mathscinet:`MR1643413` +.. [BS1969] \D. Blatt, G. Szekeres. + *A Skew Hadamard Matrix of Order 52*, + Canadian Journal of Mathematics 21 (1969): 1319-1322. + :doi:`10.4153/CJM-1969-144-2` + .. [BS1996] Eric Bach, Jeffrey Shallit. *Algorithmic Number Theory, Vol. 1: Efficient Algorithms*. MIT Press, 1996. ISBN 978-0262024051. @@ -2041,15 +2046,40 @@ REFERENCES: and some constructions of de Luca and Rauzy*, Theoret. Comput. Sci. 255 (2001) 539--553. -.. [Djo1992] \D. Đoković. +.. [Djo1992a] \D. Đoković. *Construction of some new Hadamard matrices*, Bulletin of the Australian Mathematical Society 45(2) (1992): 327-332. :doi:`10.1017/S0004972700030185` +.. [Djo1992b] \D. Đoković. + *Skew Hadamard matrices of order 4 x 37 and 4 x 43*, + Journal of Combinatorial Theory, Series A 61(2) (1992): 319-321. + :doi:`10.1016/0097-3165(92)90029-T` + +.. [Djo1992c] \D. Đoković. + *Ten New Orders for Hadamard Matrices of Skew Type*, + Publikacije Elektrotehničkog fakulteta. Serija Matematika 2 (1992): 47-59. + .. [Djo1994] \D. Đoković. *Five New Orders for Hadamard Matrices of Skew Type*, Australasian Journal of Combinatorics 10 (1994): 259-264. +.. [Djo2008a] \D. Đoković. + *Skew-Hadamard matrices of orders 188 and 388 exist*, + International Mathematical Forum 3 no.22 (2008): 1063-1068. + :arxiv:`0704.0640` + +.. [Djo2008b] \D. Đoković. + *Skew-Hadamard matrices of orders 436, 580 and 988 exist*, + Journal of Combinatorial Designs 16 (2008): 493-498. + :arxiv:`0706.1973` + +.. [Djo2023a] \D. Đoković. + *Skew-Hadamard matrices of order 276*. + :arxiv:`10.48550/ARXIV.2301.02751` + +.. [Djo2023b] \D. Đoković, Email Communication. 26 January 2023. + .. [DK2013] John R. Doyle and David Krumm, *Computing algebraic numbers of bounded height*, :arxiv:`1111.4963v4` (2013). @@ -2397,6 +2427,11 @@ REFERENCES: code, by I.A.Faradjev (with contributions by A.E.Brouwer, D.V.Pasechnik). https://github.com/dimpase/coco +.. [FKS2004] \R. J. Fletcher, C. Koukouvinos and J. Seberry. + *New skew-Hadamard matrices of order 4·59 and new D-optimal designs of order 2·59*, + Discrete Mathematics 286(3) (2004): 251-253. + :doi:`10.1016/j.disc.2004.05.009` + .. [FL2001] David Forge and Michel Las Vergnas. *Orlik-Solomon type algebras*. European J. Combin. **22** (5), (2001). pp. 699-704. @@ -5493,6 +5528,11 @@ REFERENCES: :doi:`10.1007/978-1-4684-9322-1`, ISBN 978-1-4684-9322-1. +.. [Spe1977] \E. Spence. + *Skew-Hadamard matrices of order 2(q + 1)*, + Discrete Mathematics 18(1) (1977): 79-85. + :doi:`10.1016/0012-365X(77)90009-7` + .. [Spe2013] \D. Speyer, *An infinitely generated upper cluster algebra*, :arxiv:`1305.6867`. @@ -5842,6 +5882,10 @@ REFERENCES: Hall-Littlewood vertex operators and generalized Kostka polynomials. Adv. Math. 158 (2001), no. 1, 66-85. +.. [Sze1988] \G. Szekeres. + *A note on skew type orthogonal ±1 matrices*, + Combinatorics, Colloquia Mathematica Societatis, Janos Bolyai, 52 (1988): 489-498. + .. _ref-T: **T** @@ -6122,6 +6166,11 @@ REFERENCES: pages 150--168, 1932, `available on JSTOR `_ +.. [Whi1971] \A. Whiteman. + *An infinite family of skew Hadamard matrices*, + Pacific Journal of Mathematics 38(3) (1971): 817-822. + :doi:`10.2140/pjm.1971.38.817` + .. [White2015] Noah White. *The monodromy of real Bethe vectors for the Gaudin model*. J. Combin. Algebra, **2** no. 3 (2018). pp. 259-300. diff --git a/src/sage/combinat/designs/difference_family.py b/src/sage/combinat/designs/difference_family.py index 410bb446fee..609cdf6542e 100644 --- a/src/sage/combinat/designs/difference_family.py +++ b/src/sage/combinat/designs/difference_family.py @@ -1294,24 +1294,24 @@ def _is_periodic_sequence(seq, period): first = seq[:per] periodic = True for j in range(1, len(seq)//per): - if seq[j*per:(j+1)*per] != first: + if seq[j*per : (j+1)*per] != first: periodic = False break if periodic: return False - if seq[:period] != seq[period:2*period]: + if seq[:period] != seq[period : 2*period]: return False return True def _create_m_sequence(q, n, check=True): - r"""Create an m-sequence over GF(q) with period `q^n-1`. + r"""Create an m-sequence over GF(q) with period `q^n - 1`. Given a prime power `q`, the m-sequence is created as described by [Zie1959]_ from a primitive function over the finite field `GF(q)`. - Given a primitive function `f=c_0+c_1x+...+c_nx^n` over `K=GF(q)` of degree `n`, - the recurrence is given by: `a_i = -c_0^{-1}(c_1a_{i-1}+...+c_na{i-n})`. - The first `n` elements will be `0,0,...,0,1` and these will give a maximal length recurrence sequence + Given a primitive function `f=c_0+c_1x+...+c_nx^n` over `K = GF(q)` of degree `n`, + the recurrence is given by: `a_i = -c_0^{-1}(c_1a_{i-1} + ... + c_na{i-n})`. + The first `n` elements will be `0, 0, ..., 0, 1` and these will give a maximal length recurrence sequence as shown in [Mit2008]_. INPUT: @@ -1320,7 +1320,7 @@ def _create_m_sequence(q, n, check=True): - ``n`` -- a nonnegative number. - - ``check`` -- boolean (dafault True): if true, check that the result is a seqauence with correct period. + - ``check`` -- boolean (default True): if true, check that the result is a seqauence with correct period. Setting it to false may speed up considerably the computation. EXAMPLES:: @@ -1354,13 +1354,13 @@ def _create_m_sequence(q, n, check=True): period = q**n - 1 seq_len = period*2 if check else period - seq = [1]+[0]*(n-1) + seq = [1] + [0]*(n-1) while len(seq) < seq_len: nxt = 0 for i, coeff in zip(exps[1:], coeffs[1:]): - nxt += coeff*seq[-i] - seq.append(-coeffs[0].inverse()*nxt) + nxt += coeff * seq[-i] + seq.append(-coeffs[0].inverse() * nxt) if check: assert _is_periodic_sequence(seq, period) @@ -1444,14 +1444,14 @@ def relative_difference_set_from_m_sequence(q, N, check=True): raise ValueError('N must be at least 2') m_seq = _create_m_sequence(q, N, check=False) - period = q**N-1 + period = q**N - 1 G = AdditiveAbelianGroup([period]) set1 = [i for i in G if m_seq[i[0]] == 1] if check: H = _get_submodule_of_order(G, q-1) - assert is_relative_difference_set(set1, G, H, (period//(q-1), q-1, q**(N-1), q**(N-2))) + assert is_relative_difference_set(set1, G, H, (period // (q-1), q - 1, q**(N-1), q**(N-2))) return set1 def relative_difference_set_from_homomorphism(q, N, d, check=True): @@ -1506,7 +1506,7 @@ def relative_difference_set_from_homomorphism(q, N, d, check=True): if (q-1)%d != 0: raise ValueError('q-1 must be a multiple of d') - G = AdditiveAbelianGroup([q**N-1]) + G = AdditiveAbelianGroup([q**N - 1]) K = _get_submodule_of_order(G, d) assert K is not None, 'Could not find kernel' @@ -1517,11 +1517,11 @@ def relative_difference_set_from_homomorphism(q, N, d, check=True): second_diff_set = [theta(x) for x in diff_set] if check: - H = _get_submodule_of_order(G2, (q-1)//d) - assert is_relative_difference_set(second_diff_set, G2, H, ((q**N-1)//(q-1), (q-1)//d, q**(N-1), q**(N-2)*d)) + H = _get_submodule_of_order(G2, (q-1) // d) + assert is_relative_difference_set(second_diff_set, G2, H, ((q**N-1) // (q-1), (q-1) // d, q**(N-1), q**(N-2) * d)) return second_diff_set -def is_relative_difference_set(R, G, H, params, verbose =False): +def is_relative_difference_set(R, G, H, params, verbose=False): r"""Check if `R` is a difference set of `G` relative to `H`, with the given parameters. This function checks that `G`, `H` and `R` have the orders specified in the parameters, and @@ -1546,9 +1546,9 @@ def is_relative_difference_set(R, G, H, params, verbose =False): sage: from sage.combinat.designs.difference_family import _get_submodule_of_order, relative_difference_set_from_m_sequence, is_relative_difference_set sage: q, N = 5, 2 - sage: G = AdditiveAbelianGroup([q^N-1]) - sage: H = _get_submodule_of_order(G, q-1) - sage: params = ((q^N-1)//(q-1), q-1, q^(N-1), q^(N-2)) + sage: G = AdditiveAbelianGroup([q^N - 1]) + sage: H = _get_submodule_of_order(G, q - 1) + sage: params = ((q^N-1) // (q-1), q - 1, q^(N-1), q^(N-2)) sage: R = relative_difference_set_from_m_sequence(q, N) sage: is_relative_difference_set(R, G, H, params) True @@ -1561,7 +1561,7 @@ def is_relative_difference_set(R, G, H, params, verbose =False): False """ m, n, k, d = params - if G.order() != m*n: + if G.order() != m * n: if verbose: print('Incorrect order of G:', G.order()) return False @@ -1579,7 +1579,7 @@ def is_relative_difference_set(R, G, H, params, verbose =False): for el1 in R: for el2 in R: if el1 != el2: - idx = el1-el2 + idx = el1 - el2 if idx not in diff_set: diff_set[idx] = 0 diff_set[idx] += 1 @@ -1602,11 +1602,11 @@ def is_relative_difference_set(R, G, H, params, verbose =False): return True def is_supplementary_difference_set(Ks, v, lmbda): - r"""Check that the sets in ``Ks`` are `n-\{v; k_1,...,k_n; \lambda \}` supplementary difference sets. - + r"""Check that the sets in ``Ks`` are `n-\{v; k_1, ..., k_n; \lambda \}` supplementary difference sets. + From the definition in [Spe1975]_: let `S_1, S_2, ..., S_n` be `n` subsets of an additive abelian group `G` of order `v` - such that `|S_i|= k_i`. If, for each `g\in G`, `g \neq 0`, the total number of solutions of `a_i-a'_i = g`, with - `a_i,a'_i \in S_i` is `\lambda`, then `S_1, S_2, ..., S_n` are `n-\{v; k_1,...,k_n;\lambda\}` supplementary difference sets. + such that `|S_i| = k_i`. If, for each `g \in G`, `g \neq 0`, the total number of solutions of `a_i - a'_i = g`, with + `a_i, a'_i \in S_i` is `\lambda`, then `S_1, S_2, ..., S_n` are `n-\{v; k_1, ..., k_n; \lambda\}` supplementary difference sets. INPUT: @@ -1638,8 +1638,7 @@ def is_supplementary_difference_set(Ks, v, lmbda): for K in Ks: for el1 in K: for el2 in K: - diff = G[el1]-G[el2] - + diff = G[el1] - G[el2] if diff not in differences_counter: differences_counter[diff] = 0 differences_counter[diff] += 1 @@ -1722,13 +1721,12 @@ def supplementary_difference_set(q, existence=False, check=True): .. SEEALSO:: :func:`is_supplementary_difference_set` - """ s = 0 m = -1 while q > 2**(s+1) and (q-1) % 2**(s+1) == 0: - prime_pow = (q-1)//2**(s+1)-1 + prime_pow = (q-1)//2**(s+1) - 1 if is_prime_power(prime_pow) and prime_pow % 2 == 1: m = (q - (2**(s+1) + 1)) // 2**(s+1) + 1 break @@ -1742,7 +1740,7 @@ def supplementary_difference_set(q, existence=False, check=True): if m == -1: raise ValueError('There is no s for which m-1 is an odd prime power') - set1 = relative_difference_set_from_homomorphism(m-1, 2, (m-2)//2, check=False) + set1 = relative_difference_set_from_homomorphism(m - 1, 2, (m-2) // 2, check=False) from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing P = PolynomialRing(ZZ, 'x') @@ -1753,15 +1751,15 @@ def supplementary_difference_set(q, existence=False, check=True): hall += P.monomial(d[0]) T_2m = 0 - for i in range(2*m): + for i in range(2 * m): T_2m += P.monomial(i) - modulo = P.monomial(2*m)-1 + modulo = P.monomial(2*m) - 1 diff = T_2m - (1+P.monomial(m))*hall diff = diff.mod(modulo) exp1, exp2 = diff.exponents() - a = (exp1+exp2-m)//2 + a = (exp1+exp2-m) // 2 alfa3 = (P.monomial(a) + hall).mod(modulo) alfa4 = (P.monomial(a+m) + hall).mod(modulo) @@ -1769,22 +1767,22 @@ def supplementary_difference_set(q, existence=False, check=True): psi3 = alfa3 psi4 = alfa4 for i in range(s): - psi3 = (alfa3(P.monomial(2))+P.monomial(1)*alfa4(P.monomial(2))).mod(P.monomial(4*m)-1) + psi3 = (alfa3(P.monomial(2)) + P.monomial(1)*alfa4(P.monomial(2))).mod(P.monomial(4*m)-1) psi4 = (alfa3(P.monomial(2)) + P.monomial(1)*(T_2m(P.monomial(2)) - alfa4(P.monomial(2)))).mod(P.monomial(4*m)-1) # Construction of psi1, psi2 set2 = relative_difference_set_from_m_sequence(q, 2, check=False) - s3 = _get_fixed_relative_difference_set(set2) + s3 = get_fixed_relative_difference_set(set2) phi_exps = [] for i in range(len(s3)): for j in range(i+1, len(s3)): - diff = (s3[i]-s3[j]) - if diff%(q-1) == 0 and diff%(q**2-1) != 0: + diff = s3[i] - s3[j] + if diff % (q-1) == 0 and diff % (q**2-1) != 0: phi_exps.append(s3[i]) - exps1 = [(x+1)//2 for x in phi_exps if x%2 == 1] - exps2 = [x//2 for x in phi_exps if x%2 == 0] + exps1 = [(x+1)//2 for x in phi_exps if x % 2 == 1] + exps2 = [x//2 for x in phi_exps if x % 2 == 0] theta1 = 0 for exp in exps1: @@ -1794,12 +1792,12 @@ def supplementary_difference_set(q, existence=False, check=True): theta2 = 0 for exp in exps2: theta2 += P.monomial(exp) - theta2 = theta2.mod(P.monomial(q-1)-1) + theta2 = theta2.mod(P.monomial(q-1) - 1) phi1 = 0 phi2 = 0 for exp in phi_exps: - if exp%2 == 0: + if exp % 2 == 0: phi2 += P.monomial(exp) else: phi1 += P.monomial(exp) @@ -1817,8 +1815,8 @@ def supplementary_difference_set(q, existence=False, check=True): return K1, K2, K3, K4 -def _get_fixed_relative_difference_set(rel_diff_set, as_elements=False): - r"""Contruct an equivalent relative difference set fixed by the size of the set. +def get_fixed_relative_difference_set(rel_diff_set, as_elements=False): + r"""Construct an equivalent relative difference set fixed by the size of the set. Given a relative difference set `R(q+1, q-1, q, 1)`, it is possible to find a translation of this set fixed by `q` (see Section 3 of [Spe1975]_). We say that a set is fixed by `t` if @@ -1842,30 +1840,30 @@ def _get_fixed_relative_difference_set(rel_diff_set, as_elements=False): EXAMPLES:: - sage: from sage.combinat.designs.difference_family import relative_difference_set_from_m_sequence, _get_fixed_relative_difference_set + sage: from sage.combinat.designs.difference_family import relative_difference_set_from_m_sequence, get_fixed_relative_difference_set sage: s1 = relative_difference_set_from_m_sequence(5, 2) - sage: _get_fixed_relative_difference_set(s1) #random + sage: get_fixed_relative_difference_set(s1) #random [2, 10, 19, 23, 0] If ``as_elements`` is true, the reuslt will contain elements of the group:: - sage: _get_fixed_relative_difference_set(s1, as_elements=True) #random + sage: get_fixed_relative_difference_set(s1, as_elements=True) #random [(2), (10), (19), (23), (0)] TESTS:: - sage: from sage.combinat.designs.difference_family import _is_fixed_relative_difference_set + sage: from sage.combinat.designs.difference_family import is_fixed_relative_difference_set sage: s1 = relative_difference_set_from_m_sequence(5, 2) - sage: s2 = _get_fixed_relative_difference_set(s1, as_elements=True) - sage: _is_fixed_relative_difference_set(s2, len(s2)) + sage: s2 = get_fixed_relative_difference_set(s1, as_elements=True) + sage: is_fixed_relative_difference_set(s2, len(s2)) True sage: s1 = relative_difference_set_from_m_sequence(9, 2) - sage: s2 = _get_fixed_relative_difference_set(s1, as_elements=True) - sage: _is_fixed_relative_difference_set(s2, len(s2)) + sage: s2 = get_fixed_relative_difference_set(s1, as_elements=True) + sage: is_fixed_relative_difference_set(s2, len(s2)) True sage: type(s2[0]) - sage: s2 = _get_fixed_relative_difference_set(s1) + sage: s2 = get_fixed_relative_difference_set(s1) sage: type(s2[0]) """ @@ -1875,14 +1873,14 @@ def _get_fixed_relative_difference_set(rel_diff_set, as_elements=False): s2 = None for el in G: fixed_set = [el+x for x in rel_diff_set] - if _is_fixed_relative_difference_set(fixed_set, q): + if is_fixed_relative_difference_set(fixed_set, q): s2 = fixed_set break assert s2 is not None, 'Cannot find fixed translation of the set' s3 = None for i in range(G.order()): - temp = [((q+1)*i+x[0])%G.order() for x in s2] + temp = [((q+1)*i+x[0]) % G.order() for x in s2] if 0 in temp: s3 = temp break @@ -1892,10 +1890,10 @@ def _get_fixed_relative_difference_set(rel_diff_set, as_elements=False): return [G[x] for x in s3] return s3 -def _is_fixed_relative_difference_set(R, q): +def is_fixed_relative_difference_set(R, q): r"""Check if the relative difference set `R` is fixed by `q`. - A relative difference set `R` is fixed by `q` if `\{qd | d\in R\}= R` (see Section 3 of [Spe1975]_). + A relative difference set `R` is fixed by `q` if `\{qd | d \in R\}= R` (see Section 3 of [Spe1975]_). INPUT: @@ -1905,40 +1903,45 @@ def _is_fixed_relative_difference_set(R, q): EXAMPLES:: - sage: from sage.combinat.designs.difference_family import relative_difference_set_from_m_sequence, _get_fixed_relative_difference_set, _is_fixed_relative_difference_set + sage: from sage.combinat.designs.difference_family import relative_difference_set_from_m_sequence, get_fixed_relative_difference_set, is_fixed_relative_difference_set sage: s1 = relative_difference_set_from_m_sequence(7, 2) - sage: s2 = _get_fixed_relative_difference_set(s1, as_elements=True) - sage: _is_fixed_relative_difference_set(s2, len(s2)) + sage: s2 = get_fixed_relative_difference_set(s1, as_elements=True) + sage: is_fixed_relative_difference_set(s2, len(s2)) True sage: G = AdditiveAbelianGroup([15]) sage: s3 = [G[1], G[2], G[3], G[4]] - sage: _is_fixed_relative_difference_set(s3, len(s3)) + sage: is_fixed_relative_difference_set(s3, len(s3)) False If the relative difference set does not contain elements of the group, the method returns false:: sage: s1 = relative_difference_set_from_m_sequence(7, 2) - sage: s2 = _get_fixed_relative_difference_set(s1, as_elements=False) - sage: _is_fixed_relative_difference_set(s2, len(s2)) + sage: s2 = get_fixed_relative_difference_set(s1, as_elements=False) + sage: is_fixed_relative_difference_set(s2, len(s2)) False - - """ for el in R: - if q*el not in R: + if q * el not in R: return False return True 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`. + 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 - always skew, i.e. `S_1 \cap (-S_1) = \emptyset` and `S_1 \cup (-S_1) = G\setminus\{0\}`. - - The data for `n = 103, 151` is taken from [Djo1994]_ and the data for `n = 67, 113, 127, 157, 163, 181, 241` - is taken from [Djo1992]_. - + 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 = 67, 113, 127, 157, 163, 181, 241`: [Djo1992a]_ + * `n = 37, 43`: [Djo1992b]_ + * `n = 39, 49, 65, 93, 121, 129, 133, 217, 219, 267`: [Djo1992c]_ + * `n = 97`: [Djo2008a]_ + * `n = 109, 145, 247`: [Djo2008b]_ + * `n = 73`: [Djo2023b]_ + INPUT: - ``n`` -- integer, the parameter of the supplementary difference set. @@ -1989,79 +1992,178 @@ def skew_supplementary_difference_set(n, existence=False, check=True): False sage: skew_supplementary_difference_set(127, existence=True) True - """ + .. NOTE:: + + The data for `n=247` in [Djo2008b]_ contains a typo: the set `\alpha_2` should contain `223` instead of `233`. + This can be verified by checking the resulting sets, which are given explicitly in the paper. + """ + # If -1 is present in an index set, it means that {0} should be added to that set indices = { - 67: [[0,3,5,6,9,10,13,14,17,18,20], - [0,2,4,9,11,12,13,16,19,21], - [1,3,6,10,11,13,14,16,20,21], - [2,4,6,8,9,11,14,17,19]], - 103: [[1,3,4,6,8,11,12,14,17,18,20,22,25,27,28,30,32], - [2,9,10,12,13,14,15,16,20,21,22,23,24,26,28,29,30], - [0,1,2,3,4,11,12,13,16,17,19,20,21,24,25,26,28,30,31], - [0,1,2,3,4,5,6,13,15,18,19,20,23,24,25,26,27,28,29,31]], - 113: [[0,3,4,6,8,10,13,14], - [1,3,8,9,10,11,12,13], - [0,2,3,5,6,7,12], - [1,2,3,5,8,9,15]], - 127: [[0,3,5,7,8,10,12,14,16], - [0,1,3,6,7,9,10,12,14,15], - [0,1,3,4,5,7,8,9,15,16], - [1,4,5,6,9,10,13,14,15,16]], - 151: [[0,3,5,6,8,11,13,14,16,19,21,23,25,27,28], - [2,3,6,13,16,17,20,23,25,26,27,28,29], - [0,1,2,3,4,6,7,8,9,10,11,12,23,24,27,28], - [1,4,5,10,11,12,13,14,16,18,19,22,25,26,27,28]], - 157:[[0,2,5,7,8,11], - [0,4,5,6,9,11], - [6,7,8,9,10,11], - [0,5,6,7,8,10,11]], - 163: [[0,2,5,6,9,10,13,14,17], - [0,1,7,10,12,15,16,17], - [0,1,3,5,8,13,15,16,17], - [3,6,7,8,11,12,13,14,16,17]], - 181: [[0,3,5,6,8,10,13,15,16,19], - [4,5,7,8,11,14,15,16,18,19], - [0,4,10,11,13,15,16,18,19], - [2,4,5,7,11,13,15,17,19]], - 241: [[0,2,4,6,8,11,12,14], - [1,3,4,6,7,13,14,15], - [6,8,9,10,12,13,14,15], - [3,4,5,9,10,13,14]], + 37: [[0, 3, 5, 7, 9, 10], [0, 5, 6, 7, 8], + [1, 2, 6, 7, 9], [2, 6, 8, 9, 10]], + 39: [[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]], + 43: [[1, 2, 4], [1, 2, 4], [0, 2, 3], [3, 4, -1]], + 49: [[1, 2, 5, 7, 8, 10, 13, 14], [4, 5, 6, 7, 10, 11], + [0, 1, 2, 4, 6, 7, 12, 14], [1, 2, 3, 5, 6, 10, 12, 13, 14]], + 65: [[1, 3, 5, 6, 8, 10, 13, 14, 17, 18, 20, 22], + [0, 3, 7, 10, 16, 17, 18, 20, 21], + [2, 4, 6, 8, 9, 10, 14, 15, 16, 17, 18, 20], + [5, 7, 8, 9, 11, 12, 13, 14, 16, 18, 19, 20, 21]], + 67: [[0, 3, 5, 6, 9, 10, 13, 14, 17, 18, 20], + [0, 2, 4, 9, 11, 12, 13, 16, 19, 21], + [1, 3, 6, 10, 11, 13, 14, 16, 20, 21], + [2, 4, 6, 8, 9, 11, 14, 17, 19]], + 73: [[4, 6, 8, 14], [8, 10, 12, 14], [4, 6, 10, 12], [-1, 0, 2, 10]], + 93: [[0, 3, 4, 6, 9, 10, 12, 14, 17, 18], + [2, 3, 4, 5, 9, 13, 15, 18, 19], + [1, 2, 3, 4, 5, 6, 7, 8, 16], + [1, 4, 6, 11, 12, 13, 15, 16, 17, 18]], + 97: [[1, 2, 4, 6, 9, 11, 13, 14, 17, 18, 21, 23, 25, 27, 29, 30], + [1, 2, 6, 7, 8, 9, 10, 11, 12, 13, 23, 27, 29], + [0, 1, 2, 5, 6, 12, 13, 15, 16, 20, 24, 25, 26, 29, 30, 31], + [0, 2, 3, 4, 7, 8, 9, 11, 12, 13, 15, 16, 17, 18, 23, 28, 29]], + 103: [[1, 3, 4, 6, 8, 11, 12, 14, 17, 18, 20, 22, 25, 27, 28, 30, 32], + [2, 9, 10, 12, 13, 14, 15, 16, 20, 21, 22, 23, 24, 26, 28, 29, 30], + [0, 1, 2, 3, 4, 11, 12, 13, 16, 17, 19, 20, 21, 24, 25, 26, 28, 30, 31], + [0, 1, 2, 3, 4, 5, 6, 13, 15, 18, 19, 20, 23, 24, 25, 26, 27, 28, 29, 31]], + 109: [[0, 2, 5, 7, 8, 10, 12, 15, 16, 19, 20, 23, 24, 26, 29, 30, 33, 34], + [4, 5, 6, 7, 11, 15, 18, 19, 20, 22, 25, 30, 32, 33, 35], + [0, 1, 5, 6, 9, 10, 11, 14, 17, 20, 24, 26, 27, 28, 29, 31, 32], + [0, 3, 4, 6, 7, 9, 10, 12, 13, 22, 24, 25, 26, 27, 28, 29, 31, 33, 35]], + 113: [[0, 3, 4, 6, 8, 10, 13, 14], + [1, 3, 8, 9, 10, 11, 12, 13], + [0, 2, 3, 5, 6, 7, 12], + [1, 2, 3, 5, 8, 9, 15]], + 121: [[0, 2, 4, 7, 8, 11, 13, 14, 16, 19, 20, 22], + [0, 1, 4, 5, 8, 9, 10, 15, 17, 20, 23], + [1, 2, 3, 7, 9, 16, 18, 19, 20, 21, 22, 23], + [0, 2, 9, 10, 11, 12, 13, 14, 15, 17, 18, 21, 22, 23]], + 127: [[0, 3, 5, 7, 8, 10, 12, 14, 16], + [0, 1, 3, 6, 7, 9, 10, 12, 14, 15], + [0, 1, 3, 4, 5, 7, 8, 9, 15, 16], + [1, 4, 5, 6, 9, 10, 13, 14, 15, 16]], + 129: [[1, 2, 4, 7, 9, 11, 12, 14,16,18], + [0, 1, 2, 3, 9, 11, 14, 15, 19], + [0, 1, 3, 6, 8, 10, 12, 16, 18, 19], + [0, 3, 7, 8, 9, 10, 12, 14, 15, 17]], + 133: [[1, 2, 5, 6, 9, 11, 12, 14], [1, 4, 7, 9, 10, 12, 13, 15], + [0, 5, 6, 8, 11, 12, 13, 15], [0, 1, 2, 5, 7, 8, 9, 13, 14, 15]], + 145: [[1, 2, 4, 7, 9, 10, 13, 14, 16, 19, 20, 22], [0, 2, 4, 7, 10, 11, 14, 18, 19, 20, 21, 22], + [1, 3, 6, 9, 12, 13, 14, 17, 19, 20, 21, 22, 23], [2, 3, 5, 6, 7, 9, 12, 13, 15, 16, 19, 20, 21, 22, 23]], + 151: [[0, 3, 5, 6, 8, 11, 13, 14, 16, 19, 21, 23, 25, 27, 28], + [2, 3, 6, 13, 16, 17, 20, 23, 25, 26, 27, 28, 29], + [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 23, 24, 27, 28], + [1, 4, 5, 10, 11, 12, 13, 14, 16, 18, 19, 22, 25, 26, 27, 28]], + 157: [[0, 2, 5, 7, 8, 11], + [0, 4, 5, 6, 9, 11], + [6, 7, 8, 9, 10, 11], + [0, 5, 6, 7, 8, 10, 11]], + 163: [[0, 2, 5, 6, 9, 10, 13, 14, 17], + [0, 1, 7, 10, 12, 15, 16, 17], + [0, 1, 3, 5, 8, 13, 15, 16, 17], + [3, 6, 7, 8, 11, 12, 13, 14, 16, 17]], + 181: [[0, 3, 5, 6, 8, 10, 13, 15, 16, 19], + [4, 5, 7, 8, 11, 14, 15, 16, 18, 19], + [0, 4, 10, 11, 13, 15, 16, 18, 19], + [2, 4, 5, 7, 11, 13, 15, 17, 19]], + 217: [[0, 3, 5, 7, 8, 11, 12, 14], [1, 3, 4, 7, 9, 11, 12, 15], + [3, 4, 5, 6, 7, 9, 10, 14, 15], [1, 3, 4, 5, 7, 8, 11, 13, 14]], + 219: [[1, 3, 5, 6, 8, 11, 12, 15, 17, 18, 21, 22, 24], + [2, 6, 8, 10, 11, 12, 13, 16, 19, 22, 23, 24], + [0, 1, 5, 6, 10, 11, 13, 14, 17, 20, 21, 24, 25], + [0, 2, 3, 4, 5, 6, 7, 11, 12, 13, 16, 20, 23]], + 241: [[0, 2, 4, 6, 8, 11, 12, 14], + [1, 3, 4, 6, 7, 13, 14, 15], + [6, 8, 9, 10, 12, 13, 14, 15], + [3, 4, 5, 9, 10, 13, 14]], + 247: [[0, 2, 4, 7, 8, 10, 12, 15, 16, 18, 20, 23, 25, 27, 29], + [0, 2, 7, 9, 11, 12, 14, 15, 16, 18, 20, 22, 26], + [2, 3, 4, 12, 13, 14, 15, 16, 18, 20, 23, 24, 26, 27, 29], + [0, 3, 4, 6, 10, 11, 12, 14, 18, 19, 20, 22, 25, 29]], + 267: [[0, 3, 4, 7, 8, 11, 13, 15, 16, 19, 21, 22, 25], + [0, 1, 4, 5, 6, 8, 14, 15, 18, 21, 23], + [0, 2, 4, 5, 7, 9, 10, 11, 14, 15, 16, 17, 25], + [0, 1, 3, 4, 6, 14, 15, 16, 17, 18, 20, 22, 23, 25]], } + # If the element is a list, that is the coset. cosets_gens = { - 67: [1,2,3,4,5,6,8,10,12,15,17], + 37: [1, 2, 3, 5, 6, 9], + 39: [1, 2, 3, 4, 6, 8, [13]], + 43: [1, 3, 7], + 49: [1, 2, 3, 4, 6, 7, 9, 12], + 65: [1, 2, 3, 5, 6, 7, 9, 10, 11, [13], 22, [26]], + 67: [1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 17], + 73: [1, 3, 5, 9, 11, 13, 17, 25], + 93: [1, 2, 3, 5, 7, 9, 10, 14, 15, [31]], + 97: [1, 2, 3, 4, 5, 6, 7, 9, 10, 12, 13, 15, 18, 20, 23, 26], 103: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 15, 17, 19, 21, 23, 30], + 109: [1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 13, 15, 16, 18, 20, 23, 25, 30], 113: [1, 2, 3, 5, 6, 9, 10, 13], + 121: [1, 2, 4, 5, 7, 8, 10, 11, 16, 17, 19, 20], 127: [1, 3, 5, 7, 9, 11, 13, 19, 21], - 151: [1, 2, 3 ,4, 5, 6, 9, 10, 11, 12, 15, 22, 27, 29, 30], + 129: [1, 3, 5, 7, 9, 11, 13, 19, 21, [43]], + 133: [1, 2, 3, 6, 7, 9, 18, [19, 38, 76]], + 145: [1, [2, 17, 32, 72, 77, 127, 137], [3, 43, 48, 98, 108, 118, 133], [6, 51, 71, 86, 91, 96, 121], + [7, 52, 82, 107, 112, 117, 132], [11, 21, 31, 46, 61, 101, 106], [14, 19, 69, 79, 89, 104, 119], + [22, 42, 57, 62, 67, 92, 122], [5, 35, 80, 100, 115, 120, 125], [10, 15, 55, 70, 85, 95, 105], + [29], [58]], + 151: [1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 15, 22, 27, 29, 30], 157: [1, 2, 3, 5, 9, 15], 163: [1, 2, 3, 5, 6, 9, 10, 15, 18], 181: [1, 2, 3, 4, 6, 7, 8, 12, 13, 24], + 217: [1, 2, 4, 5, 7, 10, 19, [31, 62, 124]], + 219: [1, 2, 3, 5, 7, 9, 11, 15, 19, 22, 23, 33, [73]], 241: [1, 2, 4, 5, 7, 13, 19, 35], + 247: [1, [2, 18, 31, 32, 41, 110, 122, 162, 223], [3, 27, 48, 165, 170, 183, 185, 211, 243], + [5, 28, 45, 58, 80, 158, 187, 201, 226], [6, 54, 83, 93, 96, 119, 123, 175, 239], [7, 20, 63, 73, 112, 138, 163, 180, 232], + [10, 56, 69, 90, 116, 127, 155, 160, 205], [11, 47, 99, 102, 111, 115, 150, 176, 177], [13, 52, 65, 78, 91, 117, 143, 208, 221], + [14, 29, 40, 79, 113, 126, 146, 217, 224], [17, 25, 43, 49, 140, 142, 153, 194, 225], [19, 57, 171], + [33, 34, 37, 50, 59, 86, 98, 141, 203], [35, 66, 68, 74, 100, 118, 159, 172, 196], [38, 95, 114]], + 267: [1, 2, 3, 5, 7, 9, 10, 13, 14, 15, 19, 39, [89]], } H_db = { + 37: [1, 10, -11], + 39: [1, 16, 22], + 43: [1, 4, 11, 16, 21, -2, -8], + 49: [1, 18, 30], + 65: [1, 16, 61], 67: [1, 29, 37], + 73: [1, 2, 4, 8, 16, 32, 64, 55, 37], + 93: [1, 4, 16, 64, 70], + 97: [1, 35, 61], 103: [1, 46, 56], - 113: [1,16,28,30,49,106,109], - 127: [1,2,4,8,16,32,64], - 151: [1, 8,19,59, 64], - 157: [1,14,16,39,46,67,75,93,99,101,108,130,153], - 163: [1,38,40,53,58,85,104,133,140], - 181: [1,39,43,48,62,65,73,80,132], - 241: [1,15,24,54,87,91,94,98,100,119,160,183,205,225,231], + 109: [1, 45, 63], + 113: [1, 16, 28, 30, 49, 106, 109], + 121: [1, 3, 9, 27, 81], + 127: [1, 2, 4, 8, 16, 32, 64], + 129: [1, 4, 16, 64, 97, 121, 127], + 133: [1, 4, 16, 25, 64, 93, 100, 106, 123], + 145: [1, 16, 36, 81, 111, 136, 141], + 151: [1, 8, 19, 59, 64], + 157: [1, 14, 16, 39, 46, 67, 75, 93, 99, 101, 108, 130, 153], + 163: [1, 38, 40, 53, 58, 85, 104, 133, 140], + 181: [1, 39, 43, 48, 62, 65, 73, 80, 132], + 217: [1, 8, 9, 25, 51, 64, 72, 78, 81, 142, 190, 191, 193, 200, 214], + 219: [1, 4, 16, 37, 55, 64, 148, 154, 178], + 241: [1, 15, 24, 54, 87, 91, 94, 98, 100, 119, 160, 183, 205, 225, 231], + 247: [1, 9, 16, 55, 61, 81, 139, 144, 235], + 267: [1, 4, 16, 64, 67, 91, 97, 121, 217, 223, 256], } def generate_set(index_set, cosets): S = [] for idx in index_set: - S += cosets[idx] + if idx == -1: + S.append(Z(0)) + else: + S += cosets[idx] return S - if existence: return n in indices @@ -2073,11 +2175,11 @@ def generate_set(index_set, cosets): cosets = [] for el in cosets_gens[n]: - even_coset = [] - odd_coset = [] - for x in H: - even_coset.append(x*el) - odd_coset.append(-x*el) + if isinstance(el, list): + even_coset = [Z(x) for x in el] + else: + even_coset = [x*el for x in H] + odd_coset = [-x for x in even_coset] cosets.append(even_coset) cosets.append(odd_coset) @@ -2087,7 +2189,7 @@ def generate_set(index_set, cosets): S4 = generate_set(indices[n][3], cosets) if check: - lmbda = len(S1)+len(S2)+len(S3)+len(S4) - n + 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) diff --git a/src/sage/combinat/matrices/hadamard_matrix.py b/src/sage/combinat/matrices/hadamard_matrix.py index 91f73d29101..bf915399d20 100644 --- a/src/sage/combinat/matrices/hadamard_matrix.py +++ b/src/sage/combinat/matrices/hadamard_matrix.py @@ -46,16 +46,16 @@ - [Hora]_ """ -#***************************************************************************** +# ***************************************************************************** # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # https://www.gnu.org/licenses/ -#***************************************************************************** +# ***************************************************************************** from urllib.request import urlopen -from sage.combinat.designs.difference_family import skew_supplementary_difference_set +from sage.combinat.designs.difference_family import get_fixed_relative_difference_set, relative_difference_set_from_homomorphism, skew_supplementary_difference_set from sage.rings.integer_ring import ZZ from sage.matrix.constructor import matrix, block_matrix, block_diagonal_matrix, diagonal_matrix @@ -63,34 +63,53 @@ from math import sqrt from sage.matrix.constructor import identity_matrix as I from sage.matrix.constructor import ones_matrix as J -from sage.matrix.constructor import zero_matrix +from sage.matrix.constructor import zero_matrix, matrix_method from sage.misc.unknown import Unknown from sage.cpython.string import bytes_to_str from sage.modules.free_module_element import vector from sage.combinat.t_sequences import T_sequences_smallcases +from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing -def normalise_hadamard(H): +def normalise_hadamard(H, skew=False): r""" Return the normalised Hadamard matrix corresponding to ``H``. The normalised Hadamard matrix corresponding to a Hadamard matrix `H` is a matrix whose every entry in the first row and column is +1. + If ``skew`` is True, the matrix returned will be skew-normal: a skew Hadamard + matrix with first row of all `+1`. + EXAMPLES:: - sage: from sage.combinat.matrices.hadamard_matrix import normalise_hadamard + sage: from sage.combinat.matrices.hadamard_matrix import normalise_hadamard, is_hadamard_matrix, skew_hadamard_matrix sage: H = normalise_hadamard(hadamard_matrix(4)) sage: H == hadamard_matrix(4) True + sage: H = normalise_hadamard(skew_hadamard_matrix(20, skew_normalize=False), skew=True) + sage: is_hadamard_matrix(H, skew=True, normalized=True) + True + + If ``skew`` is True but the Hadamard matrix is not skew, the matrix returned + will not be normalized:: + + sage: H = normalise_hadamard(hadamard_matrix(92), skew=True) + sage: is_hadamard_matrix(H, normalized=True) + False """ - for i in range(H.ncols()): - if H[0, i] < 0: - H.rescale_col(i, -1) - for i in range(H.nrows()): - if H[i, 0] < 0: - H.rescale_row(i, -1) - return H + + if skew: + dd = diagonal_matrix(H[0]) + return dd*H*dd + else: + for i in range(H.ncols()): + if H[0, i] < 0: + H.rescale_col(i, -1) + for i in range(H.nrows()): + if H[i, 0] < 0: + H.rescale_row(i, -1) + return H def hadamard_matrix_paleyI(n, normalize=True): @@ -140,22 +159,22 @@ def hadamard_matrix_paleyI(n, normalize=True): True """ p = n - 1 - if not(is_prime_power(p) and (p % 4 == 3)): + 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 = FiniteField(p, 'x') K_list = list(K) - K_list.insert(0,K.zero()) + 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 + H[i, 0] = -1 + H[0, i] = 1 if normalize: for i in range(n): - H[i,i] = -1 + H[i, i] = -1 H = normalise_hadamard(H) return H @@ -204,29 +223,30 @@ def hadamard_matrix_paleyII(n): True """ q = n//2 - 1 - if not(n%2==0 and is_prime_power(q) and (q % 4 == 1)): + 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 = FiniteField(q, 'x') K_list = list(K) - K_list.insert(0,K.zero()) + 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[0, i] = 1 + H[i, 0] = 1 + H[i, i] = 0 - tr = { 0: matrix(2,2,[ 1,-1,-1,-1]), - 1: matrix(2,2,[ 1, 1, 1,-1]), - -1: matrix(2,2,[-1,-1,-1, 1])} + tr = { 0: matrix(2, 2, [ 1, -1, -1, -1]), + 1: matrix(2, 2, [ 1, 1, 1, -1]), + -1: matrix(2, 2, [-1, -1, -1, 1])} - H = block_matrix(q+1,q+1,[tr[v] for r in H for v in r]) + H = block_matrix(q+1, q+1, [tr[v] for r in H for v in r]) return normalise_hadamard(H) + def hadamard_matrix_williamson_type(a, b, c, d, check=True): r""" Construction of Williamson type Hadamard matrix. @@ -280,16 +300,17 @@ def hadamard_matrix_williamson_type(a, b, c, d, check=True): n = len(a) assert len(a) == len(b) == len(c) == len(d) - assert A*A.T+B*B.T+C*C.T+D*D.T==4*n*I(n) + assert A*A.T+B*B.T+C*C.T+D*D.T == 4*n*I(n) M = block_matrix([[ A, B, C, D], [-B, A, -D, C], [-C, D, A, -B], - [-D, -C, B, A]]) + [-D, -C, B, A]]) if check: assert is_hadamard_matrix(M, normalized=False, skew=False) return M + def williamson_type_quadruples_smallcases(n, existence=False): r""" Quadruples of matrices that can be used to construct Williamson type Hadamard matrices. @@ -342,14 +363,14 @@ def williamson_type_quadruples_smallcases(n, existence=False): [1, -1, -1, 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]), + 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, 1, 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, 1, 1], + [1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, 1], + [1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1]), } if existence: @@ -360,6 +381,7 @@ def williamson_type_quadruples_smallcases(n, existence=False): a, b, c, d = map(vector, db[n]) return a, b, c, d + def williamson_hadamard_matrix_smallcases(n, existence=False, check=True): r""" Construct Williamson type Hadamard matrices for some small values of n. @@ -389,9 +411,9 @@ def williamson_hadamard_matrix_smallcases(n, existence=False, check=True): ... ValueError: The Williamson type Hadamard matrix of order 100 is not yet implemented. """ - assert n%4 == 0 + assert n % 4 == 0 - if not williamson_type_quadruples_smallcases(n//4, existence=True): + if not williamson_type_quadruples_smallcases(n//4, existence=True): if existence: return False raise ValueError("The Williamson type Hadamard matrix of order %s is not yet implemented." % n) @@ -419,25 +441,26 @@ def hadamard_matrix_156(): sage: hadamard_matrix_156() 156 x 156 dense matrix over Integer Ring... """ - a = [1, 1,-1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1] - b = [1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,-1] - c = [1, 1, 1,-1, 1, 1,-1,-1, 1, 1,-1, 1, 1] - d = [1, 1,-1, 1,-1, 1, 1, 1, 1,-1, 1,-1, 1] + a = [1, 1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1] + b = [1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1] + c = [1, 1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1] + d = [1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1] A, B, C, D = map(matrix.circulant, [a, b, c, d]) - return block_matrix([[ A, A, A, B,-B, C,-C,-D, B, C,-D,-D], - [ A,-A, B,-A,-B,-D, D,-C,-B,-D,-C,-C], - [ A,-B,-A, A,-D, D,-B, B,-C,-D, C,-C], - [ B, A,-A,-A, D, D, D, C, C,-B,-B,-C], - [ B,-D, D, D, A, A, A, C,-C, B,-C, B], - [ B, C,-D, D, A,-A, C,-A,-D, C, B,-B], - [ D,-C, B,-B, A,-C,-A, A, B, C, D,-D], - [-C,-D,-C,-D, C, A,-A,-A,-D, B,-B,-B], - [ D,-C,-B,-B,-B, C, C,-D, A, A, A, D], - [-D,-B, C, C, C, B, B,-D, A,-A, D,-A], - [ C,-B,-C, C, D,-B,-D,-B, A,-D,-A, A], - [-C,-D,-D, C,-C,-B, B, B, D, A,-A,-A]]) + return block_matrix([[ A, A, A, B, -B, C, -C, -D, B, C, -D, -D], + [ A, -A, B, -A, -B, -D, D, -C, -B, -D, -C, -C], + [ A, -B, -A, A, -D, D, -B, B, -C, -D, C, -C], + [ B, A, -A, -A, D, D, D, C, C, -B, -B, -C], + [ B, -D, D, D, A, A, A, C, -C, B, -C, B], + [ B, C, -D, D, A, -A, C, -A, -D, C, B, -B], + [ D, -C, B, -B, A, -C, -A, A, B, C, D, -D], + [-C, -D, -C, -D, C, A, -A, -A, -D, B, -B, -B], + [ D, -C, -B, -B, -B, C, C, -D, A, A, A, D], + [-D, -B, C, C, C, B, B, -D, A, -A, D, -A], + [ C, -B, -C, C, D, -B, -D, -B, A, -D, -A, A], + [-C, -D, -D, C, -C, -B, B, B, D, A, -A, -A]]) + def construction_four_symbol_delta_code_I(X, Y, Z, W): r""" @@ -500,7 +523,8 @@ def construction_four_symbol_delta_code_I(X, Y, Z, W): n = len(X) assert len(Y) == n and len(Z) == n-1 and len(W) == n-1 - autocorrelation = lambda seq, j: sum([seq[i]*seq[i+j] for i in range(len(seq)-j)]) + def autocorrelation(seq, j): + return sum([seq[i]*seq[i+j] for i in range(len(seq)-j)]) for j in range(1, n): assert sum(map(lambda seq: autocorrelation(seq, j), [X, Y, Z, W])) == 0 @@ -577,12 +601,14 @@ def construction_four_symbol_delta_code_II(X, Y, Z, W): n = len(Z) assert len(X) == n+1 and len(Y) == n+1 and len(W) == n - autocorrelation = lambda seq, j: sum([seq[i]*seq[i+j] for i in range(len(seq)-j)]) + def autocorrelation(seq, j): + return sum([seq[i]*seq[i+j] for i in range(len(seq)-j)]) + for j in range(1, n): assert sum(map(lambda seq: autocorrelation(seq, j), [X, Y, Z, W])) == 0 def alternate(seq1, seq2): - return [seq1[i//2] if i%2 == 0 else seq2[(i-1)//2] for i in range(len(seq1)+len(seq2))] + return [seq1[i//2] if i % 2 == 0 else seq2[(i-1)//2] for i in range(len(seq1)+len(seq2))] XaltZ = alternate(X, Z) Wneg = [-w for w in W] @@ -594,6 +620,7 @@ def alternate(seq1, seq2): T4 = XaltZ + alternate(Yneg, W) + [-1] return T1, T2, T3, T4 + def four_symbol_delta_code_smallcases(n, existence=False): r""" Return the 4-symobl `\delta` code of length `n` if available. @@ -628,18 +655,18 @@ def four_symbol_delta_code_smallcases(n, existence=False): db = { 1: ([1, -1], [1, 1], [1], [1]), 14: ([1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1], - [1, 1, 1,-1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1], - [1, 1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1], - [1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1]) + [1, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1], + [1, 1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1], + [1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1]) } T1, T2, T3, T4 = None, None, None, None - if n%2 == 1 and (n-1)//2 in db: + if n % 2 == 1 and (n-1)//2 in db: if existence: return True X, Y, Z, W = db[(n-1)//2] T1, T2, T3, T4 = construction_four_symbol_delta_code_I(X, Y, Z, W) - elif n%4 == 3 and (n-3)//4 in db: + elif n % 4 == 3 and (n-3) // 4 in db: if existence: return True X, Y, Z, W = db[(n-3)//4] @@ -653,7 +680,8 @@ def four_symbol_delta_code_smallcases(n, existence=False): return T1, T2, T3, T4 -def _construction_goethals_seidel_matrix(A, B ,C, D): + +def _construction_goethals_seidel_matrix(A, B, C, D): r""" Construct the Goethals Seidel matrix. @@ -703,12 +731,13 @@ def _construction_goethals_seidel_matrix(A, B ,C, D): [-1 1| 1 -1| 1 -1| 1 1] """ n = len(A[0]) - R = matrix(ZZ, n, n, lambda i,j: 1 if i+j==n-1 else 0) + R = matrix(ZZ, n, n, lambda i, j: 1 if i+j == n-1 else 0) return block_matrix([[ A, B*R, C*R, D*R], [-B*R, A, -D.T*R, C.T*R], [-C*R, D.T*R, A, -B.T*R], [-D*R, -C.T*R, B.T*R, A]]) + 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]_. @@ -797,6 +826,7 @@ def hadamard_matrix_cooper_wallis_construction(x1, x2, x3, x4, A, B, C, D, check assert is_hadamard_matrix(H) return H + def hadamard_matrix_cooper_wallis_smallcases(n, check=True, existence=False): r""" Construct Hadamard matrices using the Cooper-Wallis construction for some small values of `n`. @@ -856,7 +886,7 @@ def hadamard_matrix_cooper_wallis_smallcases(n, check=True, existence=False): ... AssertionError """ - assert n%4 == 0 and n > 0 + assert n % 4 == 0 and n > 0 db = { 67: ( @@ -868,7 +898,7 @@ def hadamard_matrix_cooper_wallis_smallcases(n, check=True, existence=False): } for T_seq_len in divisors(n//4): - will_size = n//(4* T_seq_len) + will_size = n // (4*T_seq_len) if (T_seq_len in db or T_sequences_smallcases(T_seq_len, existence=True)) and williamson_type_quadruples_smallcases(will_size, existence=True): if existence: return True @@ -891,6 +921,7 @@ def hadamard_matrix_cooper_wallis_smallcases(n, check=True, existence=False): return False raise ValueError("The Cooper-Wallis construction for Hadamard matrices of order %s is not yet implemented." % n) + def _get_baumert_hall_units(n, existence=False): r""" Construct Baumert-Hall units of size `n` from available 4-symbol `\delta` codes. @@ -935,8 +966,7 @@ def _get_baumert_hall_units(n, existence=False): ... ValueError: The Baumert-Hall units of size 200 have not yet been implemented """ - assert n%4 == 0 and n > 0 - + assert n % 4 == 0 and n > 0 delta_codes_len = n//4 if not four_symbol_delta_code_smallcases(delta_codes_len, existence=True): @@ -964,6 +994,7 @@ def _get_baumert_hall_units(n, existence=False): e4 = _construction_goethals_seidel_matrix(M4hat, M3hat, -M2hat, M1hat) return e1, e2, e3, e4 + def hadamard_matrix_turyn_type(a, b, c, d, e1, e2, e3, e4, check=True): r""" Construction of Turyn type Hadamard matrix. @@ -1019,11 +1050,11 @@ def hadamard_matrix_turyn_type(a, b, c, d, e1, e2, e3, e4, check=True): n = len(a) assert len(a) == len(b) == len(c) == len(d) - assert A*A.T+B*B.T+C*C.T+D*D.T==4*n*I(n) + assert A*A.T+B*B.T+C*C.T+D*D.T == 4*n*I(n) t4 = len(e1[0]) - assert t4 %4 == 0 - t = t4//4 + assert t4 % 4 == 0 + t = t4 // 4 # Check that e1, e2, e3, e4 are valid Baumert-Hall units for i in range(t4): @@ -1037,12 +1068,12 @@ def hadamard_matrix_turyn_type(a, b, c, d, e1, e2, e3, e4, check=True): for j in range(i+1, len(units)): assert units[i]*units[j].T + units[j]*units[i].T == 0*I(t4) - H = e1.tensor_product(A) + e2.tensor_product(B) + e3.tensor_product(C) + e4.tensor_product(D) if check: assert is_hadamard_matrix(H) return H + def turyn_type_hadamard_matrix_smallcases(n, existence=False, check=True): r""" Construct an Hadamard matrix of order `n` from available 4-symbol `\delta` codes and Williamson quadruples. @@ -1058,7 +1089,7 @@ def turyn_type_hadamard_matrix_smallcases(n, existence=False, check=True): - ``existence`` -- boolean (default False): if True, only check if matrix exists. - - ``check`` -- bolean: if True (default), check the the matrix is an Hadamard matrix before returning. + - ``check`` -- boolean: if True (default), check the the matrix is an Hadamard matrix before returning. EXAMPLES:: @@ -1080,7 +1111,7 @@ def turyn_type_hadamard_matrix_smallcases(n, existence=False, check=True): ... ValueError: The Turyn type construction for Hadamard matrices of order 64 is not yet implemented. """ - assert n%4 == 0 and n > 0 + assert n % 4 == 0 and n > 0 for delta_code_len in divisors(n//4): units_size = delta_code_len*4 @@ -1097,6 +1128,7 @@ def turyn_type_hadamard_matrix_smallcases(n, existence=False, check=True): return False raise ValueError("The Turyn type construction for Hadamard matrices of order %s is not yet implemented." % n) + def hadamard_matrix_spence_construction(n, existence=False, check=True): r"""Create an Hadamard matrix of order `n` using Spence construction. @@ -1150,7 +1182,7 @@ def hadamard_matrix_spence_construction(n, existence=False, check=True): """ from sage.combinat.designs.difference_family import supplementary_difference_set - assert n%4 == 0 and n > 0 + assert n % 4 == 0 and n > 0 q = n//4 @@ -1167,7 +1199,7 @@ def hadamard_matrix_spence_construction(n, existence=False, check=True): A3 = matrix.circulant([1 if j in S3 else -1 for j in range(q-1)]) A4 = matrix.circulant([1 if j in S2 else -1 for j in range(q-1)]) - P = matrix(ZZ, [[1 if (i + j)%(q-1) == 0 else 0 for i in range(1, q)] for j in range(1, q)]) + P = matrix(ZZ, [[1 if (i + j) % (q-1) == 0 else 0 for i in range(1, q)] for j in range(1, q)]) e = matrix([1]*(q-1)) m1 = matrix([-1]) @@ -1185,6 +1217,7 @@ def hadamard_matrix_spence_construction(n, existence=False, check=True): return H + def is_hadamard_matrix(M, normalized=False, skew=False, verbose=False): r""" Test if `M` is a Hadamard matrix. @@ -1243,6 +1276,14 @@ def is_hadamard_matrix(M, normalized=False, skew=False, verbose=False): False sage: is_hadamard_matrix(h, skew=False, verbose=True) True + sage: h = skew_hadamard_matrix(20, skew_normalize=False) + sage: is_hadamard_matrix(h, skew=True, normalized=True, verbose=True) + The matrix is not skew-normalized + False + sage: from sage.combinat.matrices.hadamard_matrix import normalise_hadamard + sage: h = normalise_hadamard(h, skew=True) + sage: is_hadamard_matrix(h, skew=True, normalized=True, verbose=True) + True """ n = M.ncols() if n != M.nrows(): @@ -1261,16 +1302,17 @@ def is_hadamard_matrix(M, normalized=False, skew=False, verbose=False): return False prod = (M*M.transpose()).dict() - if (len(prod) != n or - set(prod.values()) != {n} or - any((i, i) not in prod for i in range(n))): + if (len(prod) != n or set(prod.values()) != {n} or any((i, i) not in prod for i in range(n))): if verbose: print("The product M*M.transpose() is not equal to nI") return False if normalized: - if (set(M.row(0) ) != {1} or - set(M.column(0)) != {1}): + if skew and (set(M.row(0)) != {1}): + if verbose: + print("The matrix is not skew-normalized") + return False + elif not skew and (set(M.row(0)) != {1} or set(M.column(0)) != {1}): if verbose: print("The matrix is not normalized") return False @@ -1278,20 +1320,54 @@ def is_hadamard_matrix(M, normalized=False, skew=False, verbose=False): if skew: for i in range(n-1): for j in range(i+1, n): - if M[i,j] != -M[j,i]: + if M[i, j] != -M[j, i]: if verbose: print("The matrix is not skew") return False for i in range(n): - if M[i,i] != 1: + if M[i, i] != 1: if verbose: print("The matrix is not skew - diagonal entries must be all 1") return False return True -from sage.matrix.constructor import matrix_method + +def is_skew_hadamard_matrix(M, normalized=False, verbose=False): + r""" + Test if `M` is a skew Hadamard matrix. + + this is a wrapper around the function :func:`is_hadamard_matrix` + + INPUT: + + - ``M`` -- a matrix + + - ``normalized`` (boolean) -- whether to test if ``M`` is a skew-normalized + Hadamard matrix, i.e. has its first row filled with +1. + + - ``verbose`` (boolean) -- whether to be verbose when the matrix is not + skew Hadamard. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_skew_hadamard_matrix, skew_hadamard_matrix + sage: h = matrix.hadamard(12) + sage: is_skew_hadamard_matrix(h, verbose=True) + The matrix is not skew + False + sage: h = skew_hadamard_matrix(12) + sage: is_skew_hadamard_matrix(h) + True + sage: from sage.combinat.matrices.hadamard_matrix import normalise_hadamard + sage: h = normalise_hadamard(skew_hadamard_matrix(12), skew=True) + sage: is_skew_hadamard_matrix(h, verbose=True, normalized=True) + True + """ + return is_hadamard_matrix(M, skew=True, normalized=normalized, verbose=verbose) + + @matrix_method -def hadamard_matrix(n,existence=False, check=True): +def hadamard_matrix(n, existence=False, check=True): r""" Tries to construct a Hadamard matrix using the available methods. @@ -1379,7 +1455,7 @@ def hadamard_matrix(n,existence=False, check=True): sage: matrix.hadamard(324, existence=True) True """ - if not(n % 4 == 0) and (n > 2): + if not (n % 4 == 0) and (n > 2): if existence: return False raise ValueError("The Hadamard matrix of order %s does not exist" % n) @@ -1395,10 +1471,10 @@ def hadamard_matrix(n,existence=False, check=True): if existence: return True M = hadamard_matrix_paleyII(n) - elif n == 4 or n % 8 == 0 and hadamard_matrix(n//2,existence=True) is True: + elif n == 4 or n % 8 == 0 and hadamard_matrix(n//2, existence=True) is True: if existence: return True - had = hadamard_matrix(n//2,check=False) + had = hadamard_matrix(n//2, check=False) chad1 = matrix([list(r) + list(r) for r in had.rows()]) mhad = (-1) * had R = len(had.rows()) @@ -1425,7 +1501,7 @@ def hadamard_matrix(n,existence=False, check=True): if existence: return True M = turyn_type_hadamard_matrix_smallcases(n, check=False) - elif hadamard_matrix_spence_construction(n ,existence=True): + elif hadamard_matrix_spence_construction(n, existence=True): if existence: return True M = hadamard_matrix_spence_construction(n, check=False) @@ -1504,7 +1580,7 @@ def hadamard_matrix_www(url_file, comments=False): _rshcd_cache = {} -def regular_symmetric_hadamard_matrix_with_constant_diagonal(n,e,existence=False): +def regular_symmetric_hadamard_matrix_with_constant_diagonal(n, e, existence=False): r""" Return a Regular Symmetric Hadamard Matrix with Constant Diagonal. @@ -1573,13 +1649,13 @@ def regular_symmetric_hadamard_matrix_with_constant_diagonal(n,e,existence=False - [HX2010]_ """ - if existence and (n,e) in _rshcd_cache: - return _rshcd_cache[n,e] + 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 + _rshcd_cache[n, e] = True return True M = None @@ -1588,7 +1664,7 @@ def true(): sqn = None if is_square(n): sqn = int(sqrt(n)) - if n<0: + if n < 0: if existence: return False raise ValueError @@ -1596,31 +1672,31 @@ def true(): 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)]) + 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: + 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) + 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 = 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 = 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 = strongly_regular_graph(196, 91, 42, 42).adjacency_matrix() M = J(196) - 2*M elif n == 324: if existence: @@ -1638,7 +1714,7 @@ def true(): elif (e == 1 and sqn is not None and sqn % 4 == 2 and - strongly_regular_graph(sqn-1,(sqn-2)//2,(sqn-6)//4, + strongly_regular_graph(sqn-1, (sqn-2)//2, (sqn-6)//4, existence=True) is True and is_prime_power(ZZ(sqn + 1))): if existence: @@ -1648,30 +1724,31 @@ def true(): # 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]): + 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) is True and - regular_symmetric_hadamard_matrix_with_constant_diagonal(n2,e2,existence=True)) is True: + if (regular_symmetric_hadamard_matrix_with_constant_diagonal(n1, e1, existence=True) is True and + regular_symmetric_hadamard_matrix_with_constant_diagonal(n2, e2, existence=True)) is 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) + 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 + _rshcd_cache[n, e] = Unknown if existence: return Unknown - raise ValueError("I do not know how to build a {}-RSHCD".format((n,e))) + 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)) == {ZZ(e*sqn)} + assert set(map(sum, M)) == {ZZ(e*sqn)} return M + def RSHCD_324(e): r""" Return a size 324x324 Regular Symmetric Hadamard Matrix with Constant Diagonal. @@ -1717,15 +1794,16 @@ def RSHCD_324(e): from sage.graphs.generators.smallgraphs import JankoKharaghaniTonchevGraph as JKTG M = JKTG().adjacency_matrix() M = J(324) - 2*M - if e==-1: - M1=M[:162].T - M2=M[162:].T - M11=M1[:162] - M12=M1[162:].T - M21=M2[:162].T - M=block_matrix([[M12,-M11],[-M11.T,M21]]) + if e == -1: + M1 = M[:162].T + M2 = M[162:].T + M11 = M1[:162] + M12 = M1[162:].T + M21 = M2[:162].T + M = block_matrix([[M12, -M11], [-M11.T, M21]]) return M + def _helper_payley_matrix(n, zero_position=True): r""" Return the matrix constructed in Lemma 1.19 page 291 of [SWW1972]_. @@ -1852,23 +1930,23 @@ def rshcd_from_close_prime_powers(n): - [SWW1972]_ """ - if n%4: + if n % 4: raise ValueError("n(={}) must be congruent to 0 mod 4") - a,b = sorted([n-1,n+1],key=lambda x:-x%4) - Sa = _helper_payley_matrix(a) - Sb = _helper_payley_matrix(b) - U = matrix(a,[[int(i+j == a-1) for i in range(a)] for j in range(a)]) + a, b = sorted([n-1, n+1], key=lambda x: -x % 4) + Sa = _helper_payley_matrix(a) + Sb = _helper_payley_matrix(b) + U = matrix(a, [[int(i+j == a-1) for i in range(a)] for j in range(a)]) K = (U*Sa).tensor_product(Sb) + U.tensor_product(J(b)-I(b)) - J(a).tensor_product(I(b)) - F = lambda x:diagonal_matrix([-(-1)**i for i in range(x)]) - G = block_diagonal_matrix([J(1),I(a).tensor_product(F(b))]) - e = matrix(a*b,[1]*(a*b)) - H = block_matrix(2,[-J(1),e.transpose(),e,K]) + F = lambda x: diagonal_matrix([-(-1)**i for i in range(x)]) + G = block_diagonal_matrix([J(1), I(a).tensor_product(F(b))]) + e = matrix(a*b, [1]*(a*b)) + H = block_matrix(2, [-J(1), e.transpose(), e, K]) HH = G*H*G - assert len(set(map(sum,HH))) == 1 + assert len(set(map(sum, HH))) == 1 assert HH**2 == n**2*I(n**2) return HH @@ -1909,16 +1987,103 @@ def williamson_goethals_seidel_skew_hadamard_matrix(a, b, c, d, check=True): - [KoSt08]_ """ n = len(a) - A,B,C,D=map(matrix.circulant, [a,b,c,d]) + A, B, C, D = map(matrix.circulant, [a, b, c, d]) if check: - assert A*A.T+B*B.T+C*C.T+D*D.T==4*n*I(n) - assert A+A.T==2*I(n) + assert A*A.T+B*B.T+C*C.T+D*D.T == 4*n*I(n) + assert A+A.T == 2*I(n) M = _construction_goethals_seidel_matrix(A, B, C, D) if check: assert is_hadamard_matrix(M, normalized=False, skew=True) return M + +def skew_hadamard_matrix_spence_construction(n, check=True): + r""" + Construct skew Hadamard matrix of order `n` using Spence constrution. + + This function will construct skew Hadamard matrix of order `n=2(q+1)` where `q` is + a prime power with `q = 5` (mod 8). The construction is taken from [Spe1977]_, and the + relative difference sets are constructed using :func:`sage.combinat.designs.difference_family.relative_difference_set_from_homomorphism`. + + INPUT: + + - ``n`` -- A positive integer. + + - ``check`` -- boolean. If true (default), check that the resulting matrix is Hadamard + before returning it. + + OUTPUT: + + If `n` satisfies the requirements described above, the function returns a `n\times n` Hadamard matrix. + Otherwise, an exception is raised. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import skew_hadamard_matrix_spence_construction + sage: skew_hadamard_matrix_spence_construction(28) + 28 x 28 dense matrix over Integer Ring... + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix + sage: is_hadamard_matrix(skew_hadamard_matrix_spence_construction(12, check=False), skew=True) + True + sage: is_hadamard_matrix(skew_hadamard_matrix_spence_construction(60, check=False), skew=True) + True + sage: skew_hadamard_matrix_spence_construction(31) + Traceback (most recent call last): + ... + ValueError: The order 31 is not covered by the Spence construction. + sage: skew_hadamard_matrix_spence_construction(16) + Traceback (most recent call last): + ... + ValueError: The order 16 is not covered by the Spence construction. + """ + q = n//2 - 1 + m = (q+1)//2 + if n % 4 != 0 or not is_prime_power(q) or q % 8 != 5: + raise ValueError(f'The order {n} is not covered by the Spence construction.') + + D = relative_difference_set_from_homomorphism(q, 2, (q-1)//4, check=False) + D_fixed = get_fixed_relative_difference_set(D) + D_union = D_fixed + [q+1+el for el in D_fixed] + D_union = list(set([el % (4*(q+1)) for el in D_union])) + + def find_a(i): + for a in range(8): + if (a*(q+1)//2+i) % 8 == 0: + return a + + ai = [find_a(0), find_a(1), find_a(2), find_a(3)] + P = PolynomialRing(ZZ, 'x') + + Tm = 0 + for i in range(m): + Tm += P.monomial(i) + + Ds = [[] for i in range(4)] + + for el in D_union: + i = el % 8 + if i > 3: + continue + Ds[i].append(((el + ai[i] * m) // 8) % m) + + psis = [0 for i in range(4)] + for i in range(4): + for el in Ds[i]: + psis[i] += P.monomial(el) + + diffs = [(2*psis[i] - Tm).mod(P.monomial(m)-1) for i in range(4)] + a = [-el for el in diffs[1].coefficients()] + b = diffs[0].coefficients() + c = diffs[2].coefficients() + d = diffs[3].coefficients() + + return williamson_goethals_seidel_skew_hadamard_matrix(a, b, c, d, check=check) + + def GS_skew_hadamard_smallcases(n, existence=False, check=True): r""" Data for Williamson-Goethals-Seidel construction of skew Hadamard matrices @@ -1926,12 +2091,17 @@ def GS_skew_hadamard_smallcases(n, existence=False, check=True): Here we keep the data for this construction. Namely, it needs 4 circulant matrices with extra properties, as described in :func:`sage.combinat.matrices.hadamard_matrix.williamson_goethals_seidel_skew_hadamard_matrix` - Matrices for `n=36` and `52` are given in [GS70s]_. Matrices for `n=92` are given - in [Wall71]_. + Matrices are taken from: + + * `n=36, 52`: [GS70s]_ + * `n=92`: [Wall71]_ + * `n=188`: [Djo2008a]_ + * `n=236`: [FKS2004]_ + * `n=276`: [Djo2023a]_ Additional data is obtained from skew supplementary difference sets contained in :func:`sage.combinat.designs.difference_family.skew_supplementary_difference_set`, using the - construction described in [Djo1992]_. + construction described in [Djo1992a]_. INPUT: @@ -1957,27 +2127,29 @@ def GS_skew_hadamard_smallcases(n, existence=False, check=True): def pmtoZ(s): return [1 if x == '+' else -1 for x in s] - if existence: - return n in [36, 52, 92] or skew_supplementary_difference_set(n//4, existence=True) + db = { + 36: ['+++-+-+--', '+-++--++-', '--++++++-', '+++-++-++'], + 52: ['++++-++--+---', '-+-++----++-+', '--+-+++++-+++', '--+-+++++-+++'], + 92: ['+-------++-+-+--+++++++', '++--+--++++--++++--+--+', '++---+-+-+-++-+-+-+---+', '+----+--+--++--+--+----'], + 188: ['+----+----++-+-+---++-++--+--+++-+-+--++++-++++', + '++--+---+------++------++-+-++--+-+-+----+---++', + '+-+-++---++-+---+++---++-++-++-++-+++++-+-+----', + '+++-++-+-+---+-+++++--+-----++---+--+++++--++-+'], + 236: ['+-+---+-+-++-++---+----++-----+++++--++++-+++--+--+-+-+++-+', + '+-+---+-+-++-++---+----++-----+++++--++++-+++--+--+-+-+++-+', + '+++-++----+++-+-+++--+--++------+---+-----+--+-+--+---+----', + '++++++--+++--+---++-+-+-+---+-+----++++-++-+--++-+--+------'], + 276: ['+--+++--+-+++--+---++-+++++-+++-++-+--+---+-----+--+++-++---+-++---++', + '+-++--+-+----++-+---++++-+---+-++++++++-+---+-++++---+-++----+-+--++-', + '--+--+-++---+--++--+-+-+++-+--++---++++-+-+-+--+-++-+++++++--+--+++++', + '-+---+++-----++---+++-+++--+++++--+---+-+-++++-++++-++-++-+-+++++++++'] + } - if n == 36: - a = [ 1, 1, 1, -1, 1, -1, 1, -1, -1] - b = [ 1, -1, 1, 1, -1, -1, 1, 1, -1] - c = [-1, -1]+[1]*6+[-1] - d = [ 1, 1, 1, -1, 1, 1, -1, 1, 1] - return WGS(a, b, c, d, check=check) + if existence: + return n in db or skew_supplementary_difference_set(n//4, existence=True) - if n == 52: - a = pmtoZ('++++-++--+---') - b = pmtoZ('-+-++----++-+') - c = pmtoZ('--+-+++++-+++') - return WGS(a, b, c, c, check=check) - - if n == 92: - a = [1,-1,-1,-1,-1,-1,-1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1, 1, 1, 1, 1, 1, 1] - b = [1, 1,-1,-1, 1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1, 1,-1,-1, 1] - c = [1, 1,-1,-1,-1, 1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1, 1,-1,-1,-1, 1] - d = [1,-1,-1,-1,-1, 1,-1,-1, 1,-1,-1, 1, 1,-1,-1, 1,-1,-1, 1,-1,-1,-1,-1] + if n in db: + a, b, c, d = map(pmtoZ, db[n]) return WGS(a, b, c, d, check=check) if skew_supplementary_difference_set(n//4, existence=True): @@ -1991,9 +2163,348 @@ def pmtoZ(s): return None -_skew_had_cache={} -def skew_hadamard_matrix(n,existence=False, skew_normalize=True, check=True): +def skew_hadamard_matrix_whiteman_construction(n, existence=False, check=True): + r""" + Construct a skew Hadamard matrix of order `n=2(q+1)` where `q=p^t` is a prime power with `p \cong 5 \mod 8` and `t \cong 2 \mod 4`. + + Assuming `n` satisfies the conditions above, it is possible to construct two supplementary difference sets + `A, B` (see [Whi1971]_), and these can be used to construct a skew Hadamard matrix, as described in [BS1969]_. + + INPUT: + + - ``n`` -- A positive integer, the order of the matrix to be constructed. + + - ``existence`` -- boolean (default False). If True, only return whether the Hadamard matrix can be constructed. + + - ``check`` -- boolean: if True (default), check the the result is a skew Hadamard matrix + before returning it. + + OUTPUT: + + If ``existence`` is false, returns the skew Hadamard matrix of order `n`. It raises an error if `n` does + not satisfy the required conditions. + If ``existence`` is true, returns a boolean representing whether the matrix can be constructed or not. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import skew_hadamard_matrix_whiteman_construction + sage: skew_hadamard_matrix_whiteman_construction(52) + 52 x 52 dense matrix over Integer Ring... + sage: skew_hadamard_matrix_whiteman_construction(52, existence=True) + True + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix + sage: is_hadamard_matrix(skew_hadamard_matrix_whiteman_construction(52, check=False), skew=True) + True + sage: is_hadamard_matrix(skew_hadamard_matrix_whiteman_construction(340, check=False), skew=True) + True + sage: skew_hadamard_matrix_whiteman_construction(31) + Traceback (most recent call last): + ... + ValueError: The order 31 is not covered by the Whiteman construction. + sage: skew_hadamard_matrix_whiteman_construction(100) + Traceback (most recent call last): + ... + ValueError: The order 100 is not covered by the Whiteman construction. + sage: skew_hadamard_matrix_whiteman_construction(100, existence=True) + False + """ + + q = n // 2 - 1 + p, t = is_prime_power(q, get_data=True) + + is_order_valid = n % 4 == 0 and t > 0 and p % 8 == 5 and t % 4 == 2 + if existence: + return is_order_valid + + if not is_order_valid: + raise ValueError(f'The order {n} is not covered by the Whiteman construction.') + + from sage.rings.finite_rings.finite_field_constructor import GF + G = GF(q) + f = (q-1) // 8 + Cs = {i: [G.gen()**(8*s+i) for s in range(f)] for i in [0, 1, 2, 3, 6, 7]} + A = Cs[0] + Cs[1] + Cs[2] + Cs[3] + B = Cs[0] + Cs[1] + Cs[6] + Cs[7] + + m = n//4 - 1 + Glist = list(G) + + S = [[0 for i in range(n)] for j in range(n)] + for i in range(2*m + 1): + for j in range(2*m + 1): + S[2*m + 1 + i][2*m + 1 + j] = -1 if Glist[j] - Glist[i] in A else 1 + S[i][j] = -S[2*m + 1 + i][2*m + 1 + j] + S[2*m + 1 + j][i] = -1 if Glist[j] - Glist[i] in B else 1 + S[i][2*m + 1 + j] = -S[2*m + 1 + j][i] + S[4*m + 2][i] = -1 + S[4*m + 2][2*m + 1 + i] = 1 + S[i][4*m + 2] = 1 + S[i + 2*m + 1][4*m + 2] = -1 + for i in range(4*m + 3): + S[4*m + 3][i] = 1 + S[i][4*m + 3] = -1 + for i in range(4*m + 4): + S[i][i] = 1 + + H = matrix(S) + + if check: + assert is_hadamard_matrix(H, skew=True) + return H + + +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`, + 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. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import skew_hadamard_matrix_324 + sage: skew_hadamard_matrix_324() + 324 x 324 dense matrix over Integer Ring... + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix + sage: is_hadamard_matrix(skew_hadamard_matrix_324(), skew=True) + True + """ + from sage.symbolic.ring import SymbolicRing + from sage.rings.finite_rings.integer_mod_ring import Zmod + + R = SymbolicRing() + x = R.var('x') + Z3 = Zmod(3) + F = Z3.extension(x**4 - x**3 - 1) + H = [(F.gen()**16)**i for i in range(10)] + + cosets = [] + for i in range(8): + cosets.append([F.gen()**i * el for el in H]) + cosets.append([-F.gen()**i * el for el in H]) + + def generate_set(index_set, cosets): + S = [] + for idx in index_set: + S += cosets[idx] + return S + + S1 = generate_set([1, 2, 4, 6, 8, 10, 12, 14], cosets) + S2 = generate_set([1, 2, 3, 4, 10, 11, 13], cosets) + S3 = generate_set([4, 5, 6, 8, 12, 13, 14], cosets) + S4 = generate_set([2, 4, 5, 6, 7, 11, 12, 13, 15], cosets) + + A = matrix([[-1 if y-x in S1 else +1 for y in F] for x in F]) + B = matrix([[-1 if y-x in S2 else +1 for y in F] for x in F]) + C = matrix([[-1 if y-x in S3 else +1 for y in F] for x in F]) + D = matrix([[-1 if y-x in S4 else +1 for y in F] for x in F]) + + return _construction_goethals_seidel_matrix(A, B, C, D) + + +def skew_hadamard_matrix_from_good_matrices(a, b, c, d, check=True): + r""" + Construct skew Hadamard matrix from good matrices. + + Given good matrices `A`, `B`, `C`, `D` (`A` circulant, `B, C, D` back-circulant) they can be used to construct + a skew Hadamard matrix using the following block matrix (as described in [Sze1988]_): + + .. MATH:: + + \left(\begin{array}{rrrr} + A & B & C & D \\ + -B & A & D & -C \\ + -C & -D & A & B \\ + -D & C & -B & A + \end{array}\right) + + INPUT: + + - ``a`` -- (1,-1) list specifying the 1st row of `A`. + + - ``b`` -- (1,-1) list specifying the 1st row of `B`. + + - ``d`` -- (1,-1) list specifying the 1st row of `C`. + + - ``c`` -- (1,-1) list specifying the 1st row of `D`. + + - ``check`` -- boolean: if True (default), check that the matrix is a skew Hadamard matrix before returning it. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import skew_hadamard_matrix_from_good_matrices + sage: a, b, c, d = ([1, 1, 1, -1, -1], [1, -1, 1, 1, -1], [1, -1, -1, -1, -1], [1, -1, -1, -1, -1]) + sage: skew_hadamard_matrix_from_good_matrices(a, b, c, d) + 20 x 20 dense matrix over Integer Ring... + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix + sage: a, b, c, d = ([1, 1, 1, -1, -1], [1, -1, 1, 1, -1], [1, -1, -1, -1, -1], [1, -1, -1, -1, -1]) + sage: is_hadamard_matrix(skew_hadamard_matrix_from_good_matrices(a, b, c, d, check=False), skew=True) + True + sage: a, b, c, d = ([1, 1, 1, -1, -1], [1, -1, 1, 1, -1], [1, -1, -1, -1, -1], [1, -1, -1, -1, 1]) + sage: skew_hadamard_matrix_from_good_matrices(a, b, c, d) + Traceback (most recent call last): + ... + AssertionError + sage: a, b, c, d = ([1, 1, 1], [1, -1, 1, 1, -1], [1, -1, -1, -1, -1], [1, -1, -1, -1, -1]) + sage: skew_hadamard_matrix_from_good_matrices(a, b, c, d) + Traceback (most recent call last): + ... + AssertionError + """ + n = len(a) + m = (n-1) // 2 + + assert len(a) == len(b) == len(c) == len(d) + assert a[0] == 1 and b[0] == 1 and c[0] == 1 and d[0] == 1 + for i in range(1, m+1): + assert a[i] == -a[n-i] and b[i] == b[n-i] and c[i] == c[n-i] and d[i] == d[n-i] + + def back_circulant(row): + length = len(row) + return matrix([[row[(j+i) % length] for j in range(length)] for i in range(length)]) + + A = matrix.circulant(a) + B = back_circulant(b) + C = back_circulant(c) + D = back_circulant(d) + + H = block_matrix([[ A, B, C, D], + [-B, A, D, -C], + [-C, -D, A, B], + [-D, C, -B, A]]) + + if check: + assert is_hadamard_matrix(H, skew=True) + return H + + +def skew_hadamard_matrix_from_good_matrices_smallcases(n, existence=False, check=True): + r""" + Construct skew Hadamard matrices from good matrices for some small values of `n=4m`, with `m` odd. + + The function stores good matrices of odd orders `\le 31`, taken from [Sze1988]_. + These are used to create skew Hadamard matrices of order `4m`, `1 \le m \le 31` (`m` odd), using the function + :func:`skew_hadamard_matrix_from_good_matrices`. + + ALGORITHM: + + Given four sequences (stored in ``E_sequences``) of length `m`, they can be used to construct four `E-sequences` + of length `n=2m+1`, as follows: + + .. MATH:: + + \begin{aligned} + a &= 1, a_0, a_1, ..., a_m, -a_m, -a_{m-1}, ..., -a_0 \\ + b &= 1, b_0, b_1, ..., b_m, b_m, b_{m-1}, ..., b_0 \\ + c &= 1, c_0, c_1, ..., c_m, c_m, c_{m-1}, ..., c_0 \\ + d &= 1, d_0, d_1, ..., d_m, d_m, d_{m-1}, ..., d_0 \\ + \end{aligned} + + These E-sequences will be the first rows of the four good matrices needed to construct a skew Hadamard matrix + of order `4n`. + + INPUT: + + - ``n`` -- the order of the skew Hadamard matrix to be constructed. + + - ``existence`` -- boolean (default False). If True, only return whether the Hadamard matrix can be constructed. + + - ``check`` -- boolean: if True (default), check the the matrix is an Hadamard matrix before returning it. + + OUTPUT: + + If ``existence`` is false, returns the skew Hadamard matrix of order `n`. It raises an error if no data + is available to construct the matrix of the given order. + If ``existence`` is true, returns a boolean representing whether the matrix can be constructed or not. + + EXAMPLES:: + + sage: from sage.combinat.matrices.hadamard_matrix import skew_hadamard_matrix_from_good_matrices_smallcases + sage: skew_hadamard_matrix_from_good_matrices_smallcases(20) + 20 x 20 dense matrix over Integer Ring... + sage: skew_hadamard_matrix_from_good_matrices_smallcases(20, existence=True) + True + + TESTS:: + + sage: from sage.combinat.matrices.hadamard_matrix import is_hadamard_matrix + sage: is_hadamard_matrix(skew_hadamard_matrix_from_good_matrices_smallcases(28, check=False), skew=True) + True + sage: skew_hadamard_matrix_from_good_matrices_smallcases(140) + Traceback (most recent call last): + ... + ValueError: The Good matrices of order 35 are not yet implemented. + sage: skew_hadamard_matrix_from_good_matrices_smallcases(14) + Traceback (most recent call last): + ... + ValueError: The skew Hadamard matrix of order 14 from good matrices does not exist. + sage: skew_hadamard_matrix_from_good_matrices_smallcases(140, existence=True) + False + sage: skew_hadamard_matrix_from_good_matrices_smallcases(14, existence=True) + False + """ + E_sequences = { + 0: ['', '', '', ''], + 1: ['+', '-', '-', '+'], + 2: ['++', '-+', '--', '--'], + 3: ['++-', '++-', '+-+', '-++'], + 4: ['+++-', '+-+-', '---+', '++-+'], + 5: ['+-+--', '+++--', '-+++-', '---+-'], + 6: ['+-+---', '---+++', '+-+--+', '----+-'], + 7: ['+++++--', '-++--++', '----+-+', '-+---+-'], + 8: ['+--++-+-', '+--+----', '++----+-', '+---+-+-'], + 9: ['-+-----++', '+-+++++--', '-+----++-', '--+-+-++-'], + 10: ['+--+++++++', '++--++++-+', '--++-+-+-+', '---+++-+-+'], + 11: ['++-+-------', '+----+--+--', '+-+--++---+', '--++-+-+-++'], + 12: ['+-----+-+---', '+-++++-+-++-', '---+--++++--', '--+-+++--+--'], + 13: ['+---+-+--++-+', '+++---++-++-+', '+++-+++-++---', '+---++++-+-+-'], + 14: ['+--+----+-+-++', '+---++++-++--+', '+-+----++-+--+', '++++++---+-+-+'], + 15: ['+--++----+---+-', '-++-+---+-+++--', '++---+--+--+++-', '-++++++++--+-+-'] + } + + def pm_to_good_matrix(s, sign=1): + e1 = [1 if x == '+' else -1 for x in s] + e2 = [sign * (1 if x == '+' else -1) for x in s] + e2.reverse() + return [1] + e1 + e2 + + if not (n % 4 == 0 and (n//4) % 2 == 1): + if existence: + return False + raise ValueError("The skew Hadamard matrix of order %s from good matrices does not exist." % n) + + m = n//4 + l = (m-1) // 2 + + if existence: + return l in E_sequences + + if l not in E_sequences: + raise ValueError("The Good matrices of order %s are not yet implemented." % m) + + e1, e2, e3, e4 = E_sequences[l] + a = pm_to_good_matrix(e1, sign=-1) + b = pm_to_good_matrix(e2) + c = pm_to_good_matrix(e3) + d = pm_to_good_matrix(e4) + return skew_hadamard_matrix_from_good_matrices(a, b, c, d, check=check) + + +_skew_had_cache = {} + + +def skew_hadamard_matrix(n, existence=False, skew_normalize=True, check=True): r""" Tries to construct a skew Hadamard matrix @@ -2059,11 +2570,11 @@ def skew_hadamard_matrix(n,existence=False, skew_normalize=True, check=True): 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) + sage: skew_hadamard_matrix(356) Traceback (most recent call last): ... - ValueError: A skew Hadamard matrix of order 100 is not yet implemented. - sage: skew_hadamard_matrix(100,existence=True) + ValueError: A skew Hadamard matrix of order 356 is not yet implemented. + sage: skew_hadamard_matrix(356,existence=True) Unknown Check that :trac:`28526` is fixed:: @@ -2086,7 +2597,7 @@ def true(): M = None if existence and n in _skew_had_cache: return True - if not(n % 4 == 0) and (n > 2): + 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) @@ -2098,35 +2609,50 @@ def true(): if existence: return true() M = matrix([1]) + elif n == 324: + if existence: + return true() + M = skew_hadamard_matrix_324() + elif skew_hadamard_matrix_from_good_matrices_smallcases(n, existence=True): + if existence: + return true() + M = skew_hadamard_matrix_from_good_matrices_smallcases(n, check=False) elif is_prime_power(n - 1) and ((n - 1) % 4 == 3): if existence: return true() M = hadamard_matrix_paleyI(n, normalize=False) - + elif is_prime_power(n//2 - 1) and (n//2 - 1) % 8 == 5: + if existence: + return true() + M = skew_hadamard_matrix_spence_construction(n, check=False) + elif skew_hadamard_matrix_whiteman_construction(n, existence=True): + if existence: + return true() + M = skew_hadamard_matrix_whiteman_construction(n, check=False) elif n % 8 == 0: - if skew_hadamard_matrix(n//2,existence=True) is True: # (Lemma 14.1.6 in [Ha83]_) + if skew_hadamard_matrix(n//2, existence=True) is True: # (Lemma 14.1.6 in [Ha83]_) if existence: return true() - H = skew_hadamard_matrix(n//2,check=False) - M = block_matrix([[H,H], [-H.T,H.T]]) + H = skew_hadamard_matrix(n//2, check=False) + M = block_matrix([[H, H], [-H.T, H.T]]) - else: # try Williamson construction (Lemma 14.1.5 in [Ha83]_) - for d in divisors(n)[2:-2]: # skip 1, 2, n/2, and n + else: # try Williamson construction (Lemma 14.1.5 in [Ha83]_) + for d in divisors(n)[2:-2]: # skip 1, 2, n/2, and n n1 = n//d if is_prime_power(d - 1) and (d % 4 == 0) and (n1 % 4 == 0)\ - and skew_hadamard_matrix(n1,existence=True) is True: + and skew_hadamard_matrix(n1, existence=True) is True: if existence: return true() H = skew_hadamard_matrix(n1, check=False)-I(n1) - U = matrix(ZZ, d, lambda i, j: -1 if i==j==0 else\ - 1 if i==j==1 or (i>1 and j-1==d-i)\ + U = matrix(ZZ, d, lambda i, j: -1 if i == j == 0 else + 1 if i == j == 1 or (i > 1 and j-1 == d-i) else 0) - A = block_matrix([[matrix([0]), matrix(ZZ,1,d-1,[1]*(d-1))], - [ matrix(ZZ,d-1,1,[-1]*(d-1)), - _helper_payley_matrix(d-1,zero_position=0)]])+I(d) + A = block_matrix([[matrix([0]), matrix(ZZ, 1, d-1, [1]*(d-1))], + [matrix(ZZ, d-1, 1, [-1]*(d-1)), + _helper_payley_matrix(d-1, zero_position=0)]])+I(d) M = A.tensor_product(I(n1))+(U*A).tensor_product(H) break - if M is None: # try Williamson-Goethals-Seidel construction + if M is None: # try Williamson-Goethals-Seidel construction if GS_skew_hadamard_smallcases(n, existence=True) is True: if existence: return true() @@ -2137,15 +2663,13 @@ def true(): return Unknown raise ValueError("A skew Hadamard matrix of order %s is not yet implemented." % n) if skew_normalize: - dd = diagonal_matrix(M[0]) - M = dd*M*dd + M = normalise_hadamard(M, skew=True) if check: - assert is_hadamard_matrix(M, normalized=False, skew=True) - if skew_normalize: - assert M[0]==vector([1]*n) - _skew_had_cache[n]=True + assert is_hadamard_matrix(M, normalized=skew_normalize, skew=True) + _skew_had_cache[n] = True return M + def symmetric_conference_matrix(n, check=True): r""" Tries to construct a symmetric conference matrix @@ -2185,12 +2709,12 @@ def symmetric_conference_matrix(n, check=True): """ from sage.graphs.strongly_regular_db import strongly_regular_graph as srg try: - m = srg(n-1,(n-2)/2,(n-6)/4,(n-2)/4) + m = srg(n-1, (n-2)/2, (n-6)/4, (n-2)/4) except ValueError: raise C = matrix([0]+[1]*(n-1)).stack(matrix([1]*(n-1)).stack(m.seidel_adjacency_matrix()).T) if check: - assert (C==C.T and C**2==(n-1)*I(n)) + assert (C == C.T and C**2 == (n-1)*I(n)) return C @@ -2235,16 +2759,16 @@ def szekeres_difference_set_pair(m, check=True): B = [b for b in G if b + F.one() in sG] if check: from itertools import product, chain - assert(len(A) == len(B) == m) + assert (len(A) == len(B) == m) if m > 1: - assert(sG == set([xy[0] / xy[1] + assert (sG == set([xy[0] / xy[1] for xy in chain(product(A, A), product(B, B))])) - assert(all(F.one() / b + F.one() in sG for b in B)) - assert(not any(F.one() / a - F.one() in sG for a in A)) + assert (all(F.one() / b + F.one() in sG for b in B)) + assert (not any(F.one() / a - F.one() in sG for a in A)) return G, A, B -def typeI_matrix_difference_set(G,A): +def typeI_matrix_difference_set(G, A): r""" (1,-1)-incidence type I matrix of a difference set `A` in `G` @@ -2264,7 +2788,8 @@ def typeI_matrix_difference_set(G,A): [-1 -1 1 1 -1] """ n = len(G) - return matrix(n,n, lambda i,j: 1 if G[i]/G[j] in A else -1) + return matrix(n, n, lambda i, j: 1 if G[i]/G[j] in A else -1) + def rshcd_from_prime_power_and_conference_matrix(n): r""" @@ -2317,27 +2842,27 @@ def rshcd_from_prime_power_and_conference_matrix(n): - [WW1972]_ """ from sage.graphs.strongly_regular_db import strongly_regular_graph as srg - if is_prime_power(n) and 2==(n-1)%4: + if is_prime_power(n) and 2 == (n-1) % 4: try: - M = srg(n-2,(n-3)//2,(n-7)//4) + 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 + 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) + 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) + 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], + [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)]])