Complex Golay Pairs up to Length 28: A Search via Computer Algebra and Programmatic SAT
CComplex Golay Pairs up to Length 28: A Search viaComputer Algebra and Programmatic SAT
Curtis Bright a , Ilias Kotsireas b , Albert Heinle a , Vijay Ganesh a a University of Waterloo b Wilfrid Laurier University
Abstract
We use techniques from the fields of computer algebra and satisfiability checkingto develop a new algorithm to search for complex Golay pairs. We implement thisalgorithm and use it to perform a complete search for complex Golay pairs of lengthsup to 28. In doing so, we find that complex Golay pairs exist in the lengths 24 and 26but do not exist in the lengths 23, 25, 27, and 28. This independently verifies workdone by F. Fiedler in 2013 and confirms the 2002 conjecture of Craigen, Holzmann, andKharaghani that complex Golay pairs of length 23 don’t exist. Our algorithm is basedon the recently proposed SAT+CAS paradigm of combining SAT solvers with computeralgebra systems to efficiently search large spaces specified by both algebraic and logicalconstraints. The algorithm has two stages: first, a fine-tuned computer program usesfunctionality from computer algebra systems and numerical libraries to construct a listcontaining every sequence that could appear as the first sequence in a complex Golaypair up to equivalence. Second, a programmatic SAT solver constructs every sequence(if any) that pair off with the sequences constructed in the first stage to form a complexGolay pair. This extends work originally presented at the International Symposiumon Symbolic and Algebraic Computation (ISSAC) in 2018; we discuss and implementseveral improvements to our algorithm that enabled us to improve the efficiency of thesearch and increase the maximum length we search from length 25 to 28.
Keywords:
Complex Golay pairs; Boolean satisfiability; SAT solvers; Exhaustivesearch; Autocorrelation
1. Introduction
The sequences that are now referred to
Golay sequences or Golay pairs were firstintroduced by Golay (1949) in a groundbreaking paper on multislit spectrometry. Later,Golay (1961) made a detailed study of their elegant mathematical properties and in thispaper he referred to them as complementary series :Regardless of past or possible future applications, the writer has foundthese complementary series mathematically appealing. . .Since then, Golay pairs and their generalizations have been widely studied and appliedto a huge and surprisingly varied number of problems in engineering. For example, they
To appear in the Journal of Symbolic Computation July 27, 2019 a r X i v : . [ c s . S C ] J u l ave been applied to encoding acoustic surface waves (Tseng, 1971), spread-spectrumsystems (Krone and Sarwate, 1984), optical time domain reflectometry (Nazarathyet al., 1989), CDMA networks (Seberry et al., 2002), medical ultrasounds (Nowickiet al., 2003), train wheel detection (Donato et al., 2004), radar systems (Li et al., 2008),and wireless networks (Lomayev et al., 2017). Golay pairs consist of two sequencesand the property that defines them (roughly speaking) is the fact that one sequence’s“correlation” with itself is the inverse of the other sequence’s “correlation” with itself;see Definition 2 in Section 2 for the formal definition.Although Golay defined his complementary series over an alphabet of {± } , laterauthors have generalized the alphabet to include nonreal roots of unity such as the fourthroot of unity i (cid:66) √−
1. In this paper, we focus on the case where the alphabet is {± , ± i } .In this case the resulting sequence pairs are sometimes referred to as quadriphase , , or quaternary Golay pairs though we will simply refer to them as complex
Golay pairs. If a complex Golay pair of length n exists then we say that n is a complexGolay number .Complex Golay pairs have been extensively studied by many authors in severaldifferent contexts. Initially they were studied in the context of signal processing wheresome researchers ran searches for “polyphase complementary codes” in the process ofstudying signal encoding methods. In particular, Sivaswamy (1978) found a complexGolay pair of length 3 and Frank (1980) found complex Golay pairs of length 5 and 13and stated that 7, 9, 11, 15, and 17 were not complex Golay numbers. Complex Golaypairs were also introduced by Craigen (1994) (independently of the above works) inorder to expand the orders of Hadamard matrices attainable via ordinary Golay pairs—that are only known to exist in the lengths 2 a + b + c · b · c for integers a , b , c ≥
0. Craigenalso proved that if m and n are complex Golay numbers then 2 mn is a complex Golaynumber.An exhaustive search for complex Golay pairs up to length 13 was performed byHolzmann and Kharaghani (1994). This search verified the remark of Frank (1980) that7 and 9 are not complex Golay numbers but found that 11 is in fact a complex Golaynumber (the cause of the mistaken claim is unknown since it was based on computerprograms that were never published). Later, an exhaustive search up to length 19 wasperformed by Craigen, Holzmann, and Kharaghani (2002), showing that 14, 15, 17,and 19 are not complex Golay numbers. In addition, they reported that 21 was not acomplex Golay number and conjectured that 23 was not a complex Golay number basedon a partial search for complex Golay pairs of length 23.Another line of research is to provide explicit constructions for complex Golay pairsand a number of results of this form have been published as well. Fiedler, Jedwab, andParker (2008a,b) provide a construction that explains the existence of all known complexGolay pairs whose lengths are a power of 2, including complex Golay pairs of length 16discovered by Li and Chu (2005) to not fit into a construction given by Davis and Jedwab(1999). Gibson and Jedwab (2011) provide a construction that explains the existence ofall complex Golay pairs up to length 26 and give a table that lists the total number ofcomplex Golay pairs up to length 26. This table was produced by the mathematicianFrank Fiedler, who describes his enumeration method in a subsequent paper (Fiedler,2013) where he also reports that 27 and 28 are not complex Golay numbers.In this paper we give an enumeration method that can be used to verify Fiedler’s2able giving the total number of complex Golay pairs up to length 28. We implementedour method and obtained counts up to length 28 after about 8 . uwaterloo.ca/mathcheck . Toour knowledge, this is the first time that explicit complex Golay pairs (and their countsup to equivalence) have been published for lengths larger than 19. Lastly, we publiclyrelease our code for enumerating complex Golay pairs so that others may verify andreproduce our work; we were not able to find any other code for enumerating complexGolay pairs that was publicly available.Our result is of interest not only because of the verification we provide but alsobecause of the method we use to perform the verification. The method proceeds in twostages. In the first stage, a fine-tuned computer program performs an exhaustive searchamong all sequences that could possibly appear as the first sequence in a complex Golaypair of a given length (up to the equivalence defined in Proposition 8 of Section 2).Several filtering theorems that we describe in Section 2 allow us to discard almost allsequences from consideration. To apply these filtering theorems we use functionalityfrom the computer algebra system M APLE (Monagan et al., 2005) and the mathematicallibrary FFTW (Frigo and Johnson, 2005) as we describe in Section 3. After this filteringis completed we have a list of sequences of a manageable size such that the first sequenceof every complex Golay pair of a given length (up to equivalence) appears in the list.In the second stage, we use the programmatic SAT solver M
APLE
SAT developed byLiang et al. (2017) to determine the sequences from the first stage (if any) that can bepaired up with another sequence to form a complex Golay pair. A programmatic SATinstance is constructed from each sequence found in the first stage such that the instanceis satisfiable if and only if the sequence is part of a complex Golay pair. Furthermore,when the instance is satisfiable the assignment produced by the SAT solver determinesthe second sequence of a complex Golay pair.The “SAT+CAS” method that we use is of interest in its own right because itlinks the two previously separated fields of symbolic computation and satisfiabilitychecking. Recently there has been interest in combining methods from both fields tosolve computational problems as outlined in the invited talk at ISSAC by Ábrahám(2015) and demonstrated by the SC (satisfiability checking + symbolic computation)project initiated by Ábrahám et al. (2016). Our work fits into this paradigm and toour knowledge is the first application of a SAT solver to search for complex Golaypairs, though previous work exists that uses a SAT solver to search for other types ofcomplementary sequences like those defining Williamson matrices (Bright et al., 2016;Zulkoski et al., 2017; Bright, 2017).A preliminary version of our algorithm and results was presented at the conferenceISSAC (Bright et al., 2018b). We have since made several improvements to our algorithmthat allow us to extend our exhaustive search from all lengths n ≤
25 to all lengths n ≤
28. We describe our implementation, give timings for our searches, and compare thetimings with those of our previous algorithm in Section 4. The new timings are aboutan order of magnitude faster on the same hardware, demonstrating the effectiveness ofour improvements. In particular, our new algorithm incorporates the following updates:3 ore efficient filtering.
The preprocessing stage of our algorithm uses a filtering con-dition (Corollary 14 of Section 2) to remove a large number of sequences from con-sideration. In our original algorithm we applied this condition using a large numberof equally-spaced points z along the unit circle. In our new algorithm we apply thiscondition significantly less often but filter approximately the same number of sequences(and sometimes even more) by showing how to only use the condition for the z that arethe most likely to work. More efficient joining.
The first stage of our algorithm requires searching through twolists and finding elements of the first that can be joined with elements of the second. Inour previous algorithm this was done in a totally brute-force manner, i.e., every item inthe first list was paired with every item of the second list to see if they could be joined.In our new algorithm we sort the lists using a custom ordering and show how to find theelements that can be joined via a number of linear scans through the sorted lists (seeSection 3.2), thereby eliminating the need to consider all possible pairs.
Increased filtering.
Because the improved filtering method outlined above is so muchfaster than the previous filtering method we were able to add another round of filteringjust before generating the SAT instances (see Section 3.5). The result is that approxi-mately 40% of the SAT instances generated by our previous algorithm are now quicklyfiltered before even calling a SAT solver.Additionally we include more background details in this version of the paper. Forexample, we add Section 2.1 that gives an alternative definition of complex Golay pairsthat is often useful and we show that this definition is equivalent with the first definition.
2. Background on Complex Golay Pairs
In this section we present the background necessary to describe our method forenumerating complex Golay pairs. First, we require some preliminary definitions todefine what complex Golay pairs are. We let z denote the complex conjugate of z (whichis just the multiplicative inverse of z for z on the unit circle) and if A is a sequence welet A denote the sequence containing the conjugates of A . Definition 1 (cf. Kotsireas (2013)) . The nonperiodic autocorrelation function of asequence A = [ a , . . . , a n − ] of length n isN A ( s ) (cid:66) n − s − (cid:88) k = a k a k + s for s = , . . . , n − . Definition 2.
A pair of sequences ( A , B ) with A and B in {± , ± i } n are called a complexGolay pair if the sum of their nonperiodic autocorrelations is a constant zero for s (cid:44) ,i.e., N A ( s ) + N B ( s ) = for s = , . . . , n − . Note that if A and B are in {± , ± i } n (as we assume throughout this paper) then N A (0) + N B (0) = n by the definition of the nonperiodic autocorrelation function andthe fact that zz = z is ± ± i , explaining why s (cid:44) xample 3. ([1 , , − , [1 , i , is a complex Golay pair of length since the firstsequence A has autocorrelations N A (1) = and N A (2) = − and the second sequence Bhas autocorrelations N B (1) = and N B (2) = .2.1. Alternative definition Instead of viewing complex Golay pairs as pairs of sequences it is also possibleto view them as pairs of polynomials . If A is the sequence [ a , . . . , a n − ] we let A ( z )denote the Hall polynomial a + a z + . . . + a n − z n − which can be viewed as the finitegenerating function of A . This leads to the following alternative definition of complexGolay pairs. Definition 4.
A pair of polynomials ( A ( z ) , B ( z )) that have degrees n − and coefficientsin {± , ± i } are called a complex Golay pair if | A ( z ) | + | B ( z ) | = n for all z on the unitcircle. Example 5. (1 + z − z , + iz + z ) is a complex Golay pair of length since (cid:12)(cid:12)(cid:12) + z − z (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) + iz + z (cid:12)(cid:12)(cid:12) = for all z on the unit circle. These two definitions can be seen to be equivalent using the following lemma thatexpresses the squared absolute values of a polynomial evaluated on the unit circle interms of the autocorrelations of the sequence formed by the polynomial’s coefficients.
Lemma 6.
Let A be a complex sequence of length n. Then | A ( z ) | = N A (0) + (cid:18) n − (cid:88) s = N A ( s ) z − s (cid:19) . for all z on the unit circle.Proof. Since z is on the unit circle we have z = z − , so | A ( z ) | = A ( z ) A ( z ) = A ( z ) A ( z − ).Expanding this, we obtain that | A ( z ) | is equal to n − (cid:88) k = n − (cid:88) l = a k a l z k − l . Collecting terms of z s together for s from − n + n − n − (cid:88) s = (cid:18) n − (cid:88) j = s a j a j − s (cid:19) z s + n − (cid:88) s = (cid:18) n − (cid:88) j = s a j − s a j (cid:19) z − s so the coefficient of z s for positive s is N A ( s ) and for negative s is N A ( s ). Using the factthat N A ( s ) z s = N A ( s ) z − s this becomes N A (0) + n − (cid:88) s = (cid:18) N A ( s ) z − s + N A ( s ) z − s (cid:19) . The desired result now follows from the fact that z + z = z ) for all complex z .5or completeness we now demonstrate the equivalence of Definitions 2 and 4. Theorem 7.
A pair ( A , B ) is a complex Golay pair in the sense of Definition 2 if andonly if ( A ( z ) , B ( z )) is a complex Golay pair in the sense of Definition 4.Proof. If ( A , B ) is a complex Golay pair in the sense of Definition 2 then by Lemma 6 | A ( z ) | + | B ( z ) | = N A (0) + N B (0) + (cid:18) n − (cid:88) s = ( N A ( s ) + N B ( z )) z − s (cid:19) = n for all z on the unit circle, as required.Conversely, if ( A ( z ) , B ( z )) is a complex Golay pair in the sense of Definition 4 thenby Lemma 6 and subtracting off N A (0) + N B (0) = n we have n − (cid:88) s = (cid:16) N A ( s ) z − s + N A ( s ) z s + N B ( z ) z − s + N B ( s ) z s (cid:17) = z on the unit circle. Since any nonzero rational function has a finite number ofzeros in C the function given on the left hand side must be identically zero. In particular,the coefficients of z − s for s = . . . , n − N A ( s ) + N B ( s ) = s = . . . , n − There are certain invertible operations that preserve the property of being a complexGolay pair when applied to a sequence pair ( A , B ). These are summarized in thefollowing proposition. Proposition 8 (cf. Craigen et al. (2002)) . Let ([ a , . . . , a n − ] , [ b , . . . , b n − ]) be a com-plex Golay pair. Then the following are also complex Golay pairs:E1. (Reversal) ([ a n − , . . . , a ] , [ b n − , . . . , b ]) .E2. (Conjugate Reverse A) ([ a n − , . . . , a ] , [ b , . . . , b n − ]) .E3. (Swap) ([ b , . . . , b n − ] , [ a , . . . , a n − ]) .E4. (Scale A) ([ ia , . . . , ia n − ] , [ b , . . . , b n − ]) .E5. (Positional Scaling) ( i (cid:63) A , i (cid:63) B ) where c (cid:63) A denotes the sequence of coefficientsof the polynomial A ( cz ) , i.e., [ a , ca , c a , . . . , c n − a n − ] .Proof. Suppose ( A , B ) is a complex Golay pair and let z be on the unit circle.E1. Note that the reverse of A ( z ) is z n − A ( z − ) and | z n − | =
1. Then | z n − A ( z − ) | + | z n − B ( z − ) | = | A ( z − ) | + | B ( z − ) | = n since ( A , B ) is a complex Golay pair and z − is on the unit circle.E2. The conjugate reverse of A ( z ) is z n − A ( z − ) and | z n − A ( z − ) | = | A ( z ) | = | A ( z ) | .E3. | B ( z ) | + | A ( z ) | = | A ( z ) | + | B ( z ) | = n .E4. A ( z ) scaled by i is iA ( z ) and | iA ( z ) | = | A ( z ) | .65. | A ( iz ) | + | B ( iz ) | = n since ( A , B ) is a complex Golay pair and iz is on the unitcircle. Definition 9.
We call two complex Golay pairs ( A , B ) and ( A (cid:48) , B (cid:48) ) equivalent if ( A (cid:48) , B (cid:48) ) can be obtained from ( A , B ) using the transformations described in Proposition 8. The next lemma provides some normalization conditions that can be used whensearching for complex Golay pairs up to equivalence. Since all complex Golay pairs( A (cid:48) , B (cid:48) ) which are equivalent to a complex Golay pair ( A , B ) can easily be generatedfrom ( A , B ), it suffices to search for complex Golay pairs up to equivalence. Lemma 10 (cf. Fiedler (2013)) . Let ( A (cid:48) , B (cid:48) ) be a complex Golay pair. Then ( A (cid:48) , B (cid:48) ) isequivalent to a complex Golay pair ( A , B ) with a = a = b = and a ∈ {± , i } .Proof. We will transform a given complex Golay sequence pair ( A (cid:48) , B (cid:48) ) into an equiva-lent normalized one using the equivalence operations of Proposition 8. To start with, let A (cid:66) A (cid:48) and B (cid:66) B (cid:48) .First, we ensure that a =
1. To do this, we apply operation E4 (scale A ) enoughtimes until a = a =
1. To do this, we apply operation E5 (positional scaling)enough times until a =
1; note that E5 does not change a .Third, we ensure that a (cid:44) − i . If it is, we apply operation E1 (reversal) and E2(conjugate reverse A ) which has the effect of keeping a = a = a = i .Last, we ensure that b =
1. To do this, we apply operation E3 (swap) and thenoperation E4 (scale A ) enough times so that a = A but setting b = We now prove some useful properties that all complex Golay pairs satisfy and thatwill be exploited by our algorithm for enumerating complex Golay pairs. The propertieswill be used as filtering criteria: if a sequence does not satisfy them then we know itcannot possibly be part of a complex Golay pair and may therefore be filtered away.The next well-known lemma is one of the most powerful results of this form.
Lemma 11 (cf. Popovi´c (1991)) . Let A be a complex sequence of length n. If | A ( z ) | > n for some z on the unit circle then A is not a member of a complex Golay pair.Proof. Suppose the sequence A was a member of a complex Golay pair whose othermember was the sequence B . Since | B ( z ) | ≥
0, we must have | A ( z ) | + | B ( z ) | > n , incontradiction to Definition 4. Example 12.
The sequence A (cid:66) [1 , , cannot be a member of a complex Golay pairsince | A (1) | = is larger than n = . Fiedler (2013) derives a variation of Lemma 11 that is useful because it can beapplied knowing only around half of the entries of A . It is derived using the followingvariation of Lemma 6. Let A even be identical to A with the entries of odd index replacedby zeros and let A odd be identical to A with the entries of even index replaced by zeros.7 emma 13 (cf. Fiedler (2013)) . Let A be a complex sequence of length n. Then | A even ( z ) | + | A odd ( z ) | = N A (0) + (cid:32) n − (cid:88) s = s even N A ( s ) z − s (cid:33) . Proof.
The proof proceeds like in the proof of Lemma 6 and we derive | A even ( z ) | + | A odd ( z ) | = n − (cid:88) k , l = k , l even a k a l z k − l + n − (cid:88) k , l = k , l odd a k a l z k − l = n − (cid:88) k , l = k − l even a k a l z k − l . From this we see that the coefficient on z s for odd s is zero and the coefficient on z s foreven s is the same as in Lemma 6 from which the result follows. Corollary 14.
Let A be a complex sequence of length n. If | A even ( z ) | or | A odd ( z ) | isstrictly larger than n for some z on the unit circle then A is not a member of a complexGolay pair.Proof. If ( A , B ) is a complex Golay pair then by Lemma 13 we have | A even ( z ) | + | A odd ( z ) | + | B even ( z ) | + | B odd ( z ) | = n but if either | A even ( z ) | > n or | A odd ( z ) | > n then this cannot hold. Example 15.
The sequence A (cid:66) [1 , x , , y , , z , cannot be a member of a complexGolay pair (regardless of the values of x, y, and z) since | A even (1) | = is larger than n = .2.4. Sum-of-squares decomposition types The next lemma is useful because it allows us to write 2 n as the sum of four integersquares. It is stated by Holzmann and Kharaghani (1994) using a different notation; weuse the notation resum( A ) and imsum( A ) to represent the real and imaginary parts of thesum of the entries of A . For example, if A (cid:66) [1 , i , − i , i ] then resum( A ) = imsum( A ) = Lemma 16 (cf. Holzmann and Kharaghani (1994)) . If ( A , B ) is a complex Golay pairthen resum( A ) + imsum( A ) + resum( B ) + imsum( B ) = n . Proof.
Using Definition 4 with z = | resum( A ) + imsum( A ) i | + | resum( B ) + imsum( B ) i | = n . Since | resum( X ) + imsum( X ) i | = resum( X ) + imsum( X ) the result follows.A consequence of Lemma 16 is that every complex Golay pair generates a decom-position of 2 n into a sum of four integer squares. In fact, it typically generates severaldecompositions of 2 n into a sum of four squares. Recall that i (cid:63) A denotes positionalscaling by i (operation E5) on the sequence A . If ( A , B ) is a complex Golay pair then8pplying operation E5 to this pair k times shows that ( i k (cid:63) A , i k (cid:63) B ) is also a complexGolay pair. By using Lemma 16 on these complex Golay pairs one obtains the fact that2 n can be decomposed as the sum of four integer squares asresum( i k (cid:63) A ) + imsum( i k (cid:63) A ) + resum( i k (cid:63) B ) + imsum( i k (cid:63) B ) . For k > k =
0, 1, 2, and 3 thisproduces four distinct decompositions of 2 n into a sum of four squares.With the help of a computer algebra system (CAS) one can enumerate every possibleway that 2 n may be written as a sum of four integer squares. For example, when n = + + + = ·
23 and 1 + + + = ·
23 as well as all permutationsof the squares and negations of the integers being squared. During the first stage of ourenumeration method only the first sequence of a complex Golay pair is known, so atthat stage we cannot compute its whole sums-of-squares decomposition. However, it isstill possible to filter some sequences from consideration based on analyzing the twoknown terms in the sums-of-squares decomposition.For example, say that A is the first sequence in a potential complex Golay pair oflength 23 with resum( A ) = A ) =
5. We can immediately discard A fromconsideration because there is no way to chose the resum and imsum of B to completethe sums-of-squares decomposition of 2 n , i.e., there are no integer solutions ( x , y ) of0 + + x + y = n .
3. Enumeration Method
In this section we describe in detail the method we used to perform a completeenumeration of all complex Golay pairs up to length 28. Given a length n our goal is tofind all {± , ± i } sequences A and B of length n such that ( A , B ) is a complex Golay pair. even and A odd The first step of our method uses Fiedler’s trick of considering the entries of A ofeven index separately from the entries of A of odd index. There are approximately n / A even and A odd and there are four possible values for eachnonzero entry. Therefore there are approximately 2 · n / = n + possible sequences tocheck in this step. Additionally, by Lemma 10 we may assume the first nonzero entry ofboth A even and A odd is 1 and that the second nonzero entry of A even is not − i , decreasingthe number of sequences to check in this step by more than a factor of 4. It is quitefeasible to perform a brute-force search through all such sequences when n ≈ A even and A odd . There are an infinitenumber of z on the unit circle so it is not possible to apply Corollary 14 using allsuch z . One simple approach is to try a sufficiently large number of points z so thatwhen a point exists with | A (cid:48) ( z ) | > n (where A (cid:48) is either A even or A odd ) such a pointis usually discovered. For example, Bright et al. (2018b) tested Corollary 14 for 2 equally-spaced points z around the unit circle.Instead, in our implementation we use a nonlinear programming method to estimatethe maximum of | A (cid:48) ( z ) | for z on the unit circle. This can be done with the NLPS OLVE command of the computer algebra system M
APLE though for efficiency we use a custom9 implementation of a variant of the quadratic interpolation method described by Sunand Yuan (2006).We write | A (cid:48) ( z ) | in terms of the real variable θ by using the substitution z = e i θ anddefine f ( θ ) (cid:66) | A (cid:48) ( e i θ ) | . The quadratic interpolation method to estimate the maximumof f ( θ ) over 0 ≤ θ < π proceeds as follows:1. Evaluate f at the 2 points θ , . . . , θ where θ k (cid:66) π k and let f k (cid:66) f ( θ k ). If f k > n for some k we can immediately filter A (cid:48) by Corollary 14.2. For every value of f that is larger than its neighbours (i.e., values of f k that satisfy f k − ≤ f k and f k ≥ f k + ) we use interpolation to find the quadratic polynomial thatpasses through the points ( θ k − , f k − ), ( θ k , f k ), and ( θ k + , f k + ). Let θ ∗ be the valueof θ that maximizes the quadratic polynomial; this can be computed to be exactly12 · f k − ( θ k − θ k + ) + f k ( θ k + − θ k − ) + f k + ( θ k − − θ k ) f k − ( θ k − θ k + ) + f k ( θ k + − θ k − ) + f k + ( θ k − − θ k )and let f ∗ (cid:66) f ( θ ∗ ).3. Note that f ∗ is often a better approximation to a local maximum of f than f k was,and if f ∗ > n then we can filter A (cid:48) by Corollary 14. Otherwise we can use thepoint ( θ ∗ , f ∗ ) to derive a tighter interval in which a local maximum of f must lie.For example, if θ k < θ ∗ < θ k + and f ∗ > f k then we can repeat the previous stepexcept using the points ( θ k , f k ), ( θ ∗ , f ∗ ), and ( θ k + , f k + ).One can use this method to derive more and more accurate approximations to thelocal maxima of f though this is mostly unnecessary for our purposes as we only careabout finding a single value for θ with f ( θ ) > n . A good approximation to a localmaximum was usually found after a single interpolation step so in our implementationwe would move on to looking for another local maximum of f after repeating theinterpolation step three times. If all values of f k were examined and no points θ werefound with f ( θ ) > n then we save A (cid:48) as a sequence that could not be filtered.At the conclusion of this step we have two lists: one list L even of the A even that werenot discarded and one list L odd of the A odd that were not discarded. We now enumerate all possibilities for A by joining all possibilities for A even withall possibilities for A odd . The most straightforward way of doing this would be to simplytry all A ∈ L odd and A ∈ L even ; this was done by Bright et al. (2018b). However,because both L even and L odd can contain millions of sequences it can be inefficient to tryall possible pairings ( A , A ). However, many pairings can be eliminated by using theconditions implied by Lemma 16; we now show how to find all possible pairings thatsatisfy these conditions without needing to try all A ∈ L odd and A ∈ L even .Let ( u , u , u , u ) be an arbitrary quadruple of integers. We will describe how toefficiently find all A whose entries of odd index are in L odd , whose entries of even indexare in L even , and that satisfyresum( A ) = u , imsum( A ) = u , resum( i (cid:63) A ) = u , imsum( i (cid:63) A ) = u . (1)10or each A ∈ L odd we form the vector V (cid:66) (resum( A ) , imsum( A ) , resum( i (cid:63) A ) , imsum( i (cid:63) A ))and for each A ∈ L even we form the vector V (cid:66) ( u − resum( A ) , u − imsum( A ) , u − resum( i (cid:63) A ) , u − imsum( i (cid:63) A )) . We now show that the A that satisfy relationship (1) are exactly those formed by the A and A with V = V . Lemma 17.
Let A ∈ L odd and A ∈ L even and let A be the sequence that satisfies therelationship A ( z ) = A ( z ) + A ( z ) , i.e., A is a sequence of length n with {± , ± i } entries.Then A satisfies (1) if and only if the vectors V and V defined as above satisfy V = V .Proof. Consider the first entry of V and V . If they are equal, we haveresum( A ) = u − resum( A ) . This implies that resum( A ) = u since resum( A ) = resum( A ) + resum( A ). Similarly,we derive imsum( A ) = u , resum( i (cid:63) A ) = u , and imsum( i (cid:63) A ) = u , as required.Conversely, resum( A ) = u implies that the first entries of V and V are equal and theother equations from (1) imply that their other entries are also equal.Therefore we have translated the problem of finding the sequences A that satisfy (1)into the problem of finding A ∈ L odd and A ∈ L even that have matching vectors V and V . We can efficiently solve this problem using a string sorting algorithm asdescribed by (for example) Kotsireas, Koukouvinos, and Seberry (2009). First, we sortthe list L odd so that the vectors V formed above appear in lexicographically increasingorder when iterating through L ∈ L odd . Similarly, we sort L even in the same way so thatthe V appear in lexicographically increasing order.We can then find all V that match with V by a linear scan through the lists L and L and this requires only one pass through each list. For example, if at some pointwe find that V is lexicographically greater than V then we try the next L in the list L even and if we instead find that V is lexicographically greater than V then we try thenext L in the list L odd . If instead we find V = V we iterate through the L and L thatshare the same value of V and V and save the A formed by joining all such A and A in a new list L A .It remains to determine the appropriate values of ( u , u , u , u ) to use in (1). Todo this we use a quadratic Diophantine equation solver and find all solutions of theDiophantine equation u + v + x + y = n in integers u , v , x , y . Let U be the set of all pairs ( u , v ) for which this equation is solvable. By Lemma 16 andoperation E5 we know that if A is a complex Golay pair then (resum( A ) , imsum( A )) and(resum( i (cid:63) A ) , imsum( i (cid:63) A )) must be members of U . Therefore, we apply the abovejoining procedure for all (( u , u ) , ( u , u )) ∈ U .11t the conclusion of this stage we will have a list of sequences L A that couldpotentially be a member of a complex Golay pair. By construction, the first member ofall complex Golay pairs (up to the equivalence described in Lemma 10) of length n willbe in L A . To decrease the size of L A even further we perform additional filtering beforeadding sequences into L A , for example, we ensure that (cid:0) resum( − (cid:63) A ) , imsum( − (cid:63) A ) (cid:1) and (cid:0) resum( − i (cid:63) A ) , imsum( − i (cid:63) A ) (cid:1) are both in U for each A added to L A . We also filter those A with | A ( z ) | > n for some z on the unit circle (see Section 3.5 for optimization details). In the second stage we take as input the list L A generated in the first stage, i.e., a listof the sequences A that were not filtered by any of the filtering theorems we applied. Foreach A ∈ L A we attempt to construct a second sequence B such that ( A , B ) is a complexGolay pair. We do this by generating a SAT instance that encodes the property of ( A , B )being a complex Golay pair where the entries of A are known and the entries of B areunknown and encoded using Boolean variables. Because there are four possible valuesfor each entry of B we use two Boolean variables to encode each entry. Although theexact encoding used is arbitrary, we fixed the following encoding in our implementation,where the variables v k and v k + represent b k , the k th entry of B : v k v k + b k F F F T − T F i T T − i To encode the property that ( A , B ) is a complex Golay pair in our SAT instance weadd the conditions that define ( A , B ) to be a complex Golay pair, i.e., N A ( s ) + N B ( s ) = s = , . . . , n − . These equations could be encoded using clauses in conjunctive normal form (for exampleby constructing logical circuits to perform complex multiplication and addition andthen converting those circuits into CNF clauses). However, we found that a much moreefficient and convenient method was to use a programmatic
SAT solver.The concept of a programmatic SAT solver was first introduced by Ganesh et al.(2012) where a programmatic SAT solver was shown to be more efficient than a standardSAT solver when solving instances derived from RNA folding problems. More recently,a programmatic SAT solver was also shown to be useful when searching for Williamsonmatrices and good matrices by Bright et al. (2018a, 2019). Generally, programmaticSAT solvers perform well when there is additional domain-specific knowledge knownabout the problem being solved. Often this knowledge cannot easily be encoded into aSAT instance directly but can be given to a programmatic SAT solver to help guide thesolver in its search. 12oncretely, a programmatic SAT solver is compiled with a piece of code that encodesa property any solution must satisfy. Periodically the SAT solver will run this codewhile performing its search, and if the current partial assignment violates a property thatis expressed in the provided code then a conflict clause is generated encoding this fact.The conflict clause is added to the SAT solver’s database of learned clauses where itis used to increase the efficiency of the remainder of the search. The reason that theseclauses can be so useful is because they can encode facts that the SAT solver would haveno way of learning otherwise, since the SAT solver has no knowledge of the domain ofthe problem.Not only does this paradigm allow the SAT solver to perform its search moreefficiently, it also allows instances to be much more expressive. Under this frameworkSAT instances do not have to consist solely of Boolean formulas in conjunctive normalform (the typical format of SAT instances) but can consist of clauses in conjunctivenormal form combined with a piece of code that programmatically expresses clauses.Increased expressiveness is also a feature of SMT (SAT modulo theories) solvers, thoughSMT solvers typically require additional overhead and only support a fixed number oftheories such as those specified in the SMT library (Barrett et al., 2016). Additionally,one can compile instance-specific programmatic SAT solvers that are tailored to performsearches for a specific class of problems.For our purposes we use a programmatic SAT solver tailored to search for se-quences B that when paired with a given sequence A form a complex Golay pair. Eachinstance will contain the 2 n variables v , . . . , v n − that encode the entries of B aspreviously specified. In detail, the code given to the SAT solver does the following:1. Compute and store the values N A ( k ) for k = . . . , n − s to n −
1. This will be a variable that controls which autocorrelationcondition we are currently examining.3. Examine the current partial assignment to v , v , v n − , and v n − . If all thesevalues have been assigned then we can determine the values of b and b n − . Fromthese values we compute N B ( s ) = b b n − . If N A ( s ) + N B ( s ) (cid:44) A , B ) cannotbe a complex Golay pair (regardless of the values of b , . . . , b n − ) and thereforewe learn a conflict clause saying that b and b n − cannot both be assigned to theircurrent values. More explicitly, if v cur k represents the literal v k when v k is currentlyassigned to true and the literal ¬ v k when v k is currently assigned to false we learnthe clause ¬ v cur0 ∨ ¬ v cur1 ∨ ¬ v cur2 n − ∨ ¬ v cur2 n − that says that at least one of { v , v , v n − , v n − } must have their value changed.4. Decrement s by 1 and repeat the previous step, computing N B ( s ) if the all the b k that appear in its definition have known values. If N A ( s ) + N B ( s ) (cid:44) b k that appear in the definition of N B ( s ) frombeing assigned the way that they currently are. Continue to repeat this step until s = B are assigned but no clauses have been learned in steps 3–4 thenoutput the complex Golay pair ( A , B ). If an exhaustive search is desired, learn aclause preventing the values of B from being assigned the way they currently are;otherwise learn nothing and return control to the SAT solver.13or each A in the list L A from stage 1 we run a SAT solver with the above programmaticcode. The list of all outputs ( A , B ) in step 5 shown above now form a complete list ofcomplex Golay pairs of length n up to the equivalence given in Lemma 10. In fact, sinceLemma 10 says that we can set b = v and v are alwaysset to false. In other words, we can add the two unit clauses ¬ v and ¬ v into our SATinstance without omitting any complex Golay pairs up to equivalence. At the conclusion of the second stage we have obtained a list of complex Golaypairs of length n such that every complex Golay pair of length n is equivalent to somepair in our list. However, because we have not accounted for all the equivalences inSection 2.2 some pairs in our list may be equivalent to each other. In some sense suchpairs should not actually be considered distinct, so to count how many distinct complexGolay pairs exist in length n we would like to find and remove pairs that are equivalentfrom the list. Additionally, to verify the counts given by Gibson and Jedwab (2011) it isnecessary to produce a list that contains all complex Golay pairs. We now describe analgorithm that does both, i.e., it produces a list of all complex Golay pairs as well as alist of all inequivalent complex Golay pairs.In detail, our algorithm performs the following steps:1. Initialize Ω all to be the empty set. This variable will be a set that will be populatedwith and eventually contain all complex Golay pairs of length n .2. Initialize Ω inequiv to be the empty set. This variable will be a set that will bepopulated with and eventually contain all inequivalent complex Golay pairs oflength n .3. For each ( A , B ) in the set of complex Golay pairs generated in stage 2:(a) If ( A , B ) is already in Ω all then skip this ( A , B ) and proceed to the next pair.(b) Initialize Γ to be the set containing ( A , B ). This variable will be a set that willbe populated with and eventually contain all complex Golay pairs equivalentto ( A , B ).(c) For every γ in Γ add E1( γ ), . . . , E5( γ ) to Γ . Continue to do this until everypair in Γ has been examined and no new pairs are added to Γ .(d) Add ( A , B ) to Ω inequiv and add all pairs in Γ to Ω all .Additionally, pseudocode for this algorithm is given in Algorithm 3.1. Algorithm 3.1
Pseudocode for our postprocessing algorithm. Note that the sets beingiterated through will update during each iteration. Ω all (cid:66) ∅ , Ω inequiv (cid:66) ∅ for ω ∈ { list of generated complex Golay pairs } \ Ω all do Γ (cid:66) { ω } for γ ∈ Γ do Γ (cid:66) Γ ∪ { E1( γ ) , . . . , E5( γ ) } Ω inequiv (cid:66) Ω inequiv ∪ { ω } Ω all (cid:66) Ω all ∪ Γ Ω all gives a list of all complexGolay pairs of length n and listing the members of Ω inequiv gives a list of all inequivalentcomplex Golay pairs of length n . At this point we can also construct the complete list ofsequences that appear in any complex Golay pair of length n . To do this it suffices toadd A and B to a new set Ω seqs for each ( A , B ) ∈ Ω all . Although the method described will correctly enumerate all complex Golay pairs of agiven length n , for the benefit of potential implementors we mention a few optimizationsthat we found helpful.In the preprocessing step and stage 1 it is necessary to evaluate a polynomial atpoints on the unit circle and determine its squared absolute value. The fastest way wefound to do this used the discrete Fourier transform (DFT). For example, let A (cid:48) be thesequence A even , A odd , or A under consideration but padded with trailing zeros so that A (cid:48) is of length N . By definition of the discrete Fourier transform we have that the j th entryof DFT( A (cid:48) ) is exactly A (cid:48) ( z j ) where z j (cid:66) exp(2 π i j / N ). Thus, we determine the values of | A (cid:48) ( z j ) | by taking the squared absolute values of the entries of DFT( A (cid:48) ). If | A (cid:48) ( z ) | > n for some z then by Lemma 11 or Corollary 14 we can discard A (cid:48) from consideration.To guard against potential inaccuracies introduced by the algorithms used to computethe DFT we actually ensure that | A (cid:48) ( z ) | > n + (cid:15) for some tolerance (cid:15) that is small butlarger than the accuracy of the DFT (e.g., (cid:15) = − ).In stage 1 we solve the quadratic Diophantine equation u + v + x + y = n inintegers u , v , x , y . In fact, we can also add the constraints u + v ≡ n (mod 2) and x + y ≡ n (mod 2)(the second is unnecessary, being implied by the first) because of the following lemma. Lemma 18.
Suppose R and I are the resum and imsum of a sequence A ∈ {± , ± i } n .Then R + I ≡ n (mod 2) .Proof. Let c denote the number of entries in A with value c . Then R + I = ( − − ) + ( i − − i ) ≡ + − + i + − i (mod 2)since − ≡ n since there are n entries in A .In the final filtering step of stage 1 we ensure that Diophantine equations of the form R + I + x + y = n are solvable in integers ( x , y ) where R and I are given. This canbe efficiently done by precomputing a Boolean two dimensional array D such that D | R | , | I | is true if and only if this equation has a solution, making the check for solvability a fastlookup. At this point we also check if we can find a z on the unit circle with | A ( z ) | > n .Because the joining process is the bottleneck of our algorithm it is important to makethis filtering step as efficient as possible; we did this by computing the values of | A ( z ) | via the expression | A ( z ) + A ( z ) | and precomputing the values of A ( z ) and A ( z ) for A ∈ L odd and A ∈ L even on points of the form z = exp(2 π i j /
32) with j = , . . . , j in ascendingorder (i.e., for each z in ascending complex argument) but to first perform the check15n points z with larger spacing between them. In our implementation we checked thecondition for z of the form exp(2 π i j / N ) with j odd and N =
8, 16, and 32 in thatorder. (This ignores checking the condition when z ∈ {± , ± i } but that is desirable since | A ( i k ) | = resum( i k (cid:63) A ) + imsum( i k (cid:63) A ) and the sums-of-squares condition is a strictlystronger filtering method.)The downside of precomputing the values of A ( z ) and A ( z ) for A ∈ L odd and A ∈ L even is that when the lists L odd and L even are large this requires a significantamount of memory. In our implementation storing the values of A ( z ) for all A ∈ L even in the lengths 27 and 28 each required about 8GB of RAM, the searches in lengths25 and 26 each required about 2GB of RAM, and the searches in lengths 23 and 24each required about 0.5GB of RAM. However, the amount of precomputation could beincreased or decreased depending on the amount of memory available.In stage 1 we make a pass through the lists L odd and L even for each valid possibilityof the values ( u , u , u , u ). However, it is possible to reuse some work between passes.For example, it is only necessary to sort the lists L odd and L even once. Even thoughthe values of V depend on the values ( u , u , u , u ) the relative ordering of L even isunaffected since all V are shifted by the same amount in each pass. It is even possibleto do some parts of the passes simultaneously: for example, if the first coordinate of V is greater than the first coordinate of V for some u then V is lexicographically greaterthan V in all passes with that value of u and so we can iterate to the next L in all suchpasses.At the conclusion of stage 1 we added an additional round of filtering to the se-quences in L A that was not in our previous work. Since this step was not a bottleneckwe used a more involved filtering process here; in our implementation we ensured that | A ( z ) | ≤ n for the 2 points of the form z j = exp(2 π i j / ). Additionally, we keeptrack of the local maxima found by this process and use the quadratic interpolationmethod described in Section 3.1 to find more accurate approximations to the localmaxima of | A ( z ) | .In stage 2 one can also include properties that complex Golay sequences must satisfyin the SAT instances. As an example of this, we state the following proposition that waspublished by Bright et al. (2018b) at ISSAC 2018. Proposition 19.
Let ( A , B ) be a complex Golay pair. Thena k a n − k − b k b n − k − = ± for k = , . . . , n − . To prove this, we use the following simple lemma.
Lemma 20.
Let c k ∈ Z for k = , . . . , n − . Then n − (cid:88) k = i c k = implies n − (cid:88) k = c k ≡ . Proof.
Let c denote the number of c k with value c . Note that the sum on the left impliesthat = and = because the 1s must cancel with the −
1s and the i s must cancelwith the − i s. Then (cid:80) n − k = c k = + + = + ≡ roof. Let c k , d k ∈ Z be such that a k = i c k and b k = i d k . Using this notation themultiplicative equation from Proposition 19 becomes the additive congruence c k + c n − k − + d k + d n − k − ≡ . (2)Since ( A , B ) is a complex Golay pair, the autocorrelation equations give us n − s − (cid:88) k = (cid:16) i c k − c k + s + i d k − d k + s (cid:17) = s = . . . , n −
1. Using Lemma 20 and the fact that − ≡ n − s − (cid:88) k = (cid:0) c k + c k + s + d k + d k + s (cid:1) ≡ s = . . . , n −
1. With s = n − s =
1) one immediately derives (2) for k = s = n − s = n − k = s = n − s = n − k = k .In short, Proposition 19 tells us that an even number of a k , a n − k − , b k , and b n − k − arereal for each k = . . . , n −
1. For example, if exactly one of a k and a n − k − is real thenexactly one of b k and b n − k − must also be real. In this case, using our encoding fromSection 3.3 we add the two binary clauses v k ∨ v n − k − and ¬ v k ∨ ¬ v n − k − to our SAT instance. These clauses say that exactly one of v k and v n − k − is true.Conversely, if an even number of a k and a n − k − are real then an even number of b k and b n − k − must also be real. In this case we add the two binary clauses v k ∨ ¬ v n − k − and ¬ v k ∨ v n − k − to our SAT instance.
4. Results
In order to provide a verification of the counts given by Fiedler (2013), Gibson andJedwab (2011), and Craigen, Holzmann, and Kharaghani (2002) we implemented theenumeration method described in Section 3. The preprocessing step was performed bya C program and used the mathematical library FFTW by Frigo and Johnson (2005)for computing the values of f , . . . , f as described in Section 3.5. Stage 1 wasperformed by a C++ program, used FFTW for computing the values of A ( z ) and A ( z ),and a M APLE script by Riel (2006) for determining the solvability of the Diophantineequations given in Section 3.3. Stage 2 was performed by the programmatic SAT solverM
APLE
SAT by Liang et al. (2017). The postprocessing step was performed by a Pythonscript. 17otal CPU Time in hours n Preproc. Stage 1 Stage 217 0.00 0.00 0.0018 0.00 0.01 0.0119 0.00 0.01 0.0120 0.00 0.08 0.0421 0.00 0.71 0.1322 0.00 0.88 0.1523 0.01 4.67 0.1924 0.01 7.09 1.2425 0.03 97.86 1.3726 0.06 234.43 4.7227 0.12 3255.56 49.5628 0.22 2543.34 25.73
Table 1: The time used to run the various stages of our algorithm in lengths 17 ≤ n ≤ We ran our implementation on a cluster of machines running CentOS 7 and usingIntel Xeon E5-2683V4 processors running at 2 . N cores by partitioning the list L odd into N sublists of approximately equal size. Everything in the stages proceededexactly as before except that in stage 1 the list L odd was about N times shorter than itwould have been otherwise. In practice, this allowed us to complete the search nearly N times faster; for example, our search in length 28 required about 3.5 months of CPUtime but completed in about 3 hours when distributed across 1000 cores. The timingsfor the preprocessing step and the two stages of our algorithm are given in Table 1;the timings for the postprocessing step were negligible. The times are given as thetotal amount of CPU time used across all cores. Our code is available online as apart of the M ATH C HECK project (see uwaterloo.ca/mathcheck ) and the resultingenumeration of complex Golay pairs is available from https://doi.org/10.5281/zenodo.1246337 .The sizes of the lists L even and L odd computed in the preprocessing step and the sizeof the list L A computed in stage 1 are given in Table 2 for all lengths n ≤
28. Withoutapplying any filtering L A would have size 4 n so Table 2 demonstrates the power of thecriteria we used to perform filtering; for example in length 28 less than 10 − % of thepossible sequences A were added to L A . A result of Gersho et al. (1979) implies thatfor z on the unit circle the maximum value of | A ( z ) | is n log n + O ( n log log n ) for almostall {± , ± i } -polynomials A of degree n . In particular, the number of polynomials that donot satisfy this is at most 4 n / (log n ) , implying that the size of L A will be o (4 n ) thoughin practice we even have | L A | ≤ n .The generated SAT instances had 2 n variables (encoding the entries b , . . . , b n − ),two unit clauses (encoding b = (cid:98) n / (cid:99) binary clauses (encoding Proposition 19),18nd n − b k are known for k = . . . , (cid:98) n / (cid:99) , the programmatic constraints uniquely determine theremaining values of b k . Thus 4 n / is a crude upper bound on the number of possibleassignments to the values of b , . . . , b n − .We now give an upper bound on the complexity of the parts of our algorithm that donot run in polynomial time. The preprocessing and postprocessing costs are dominatedby the costs of stage 1 and 2. The cost of stage 1 is at most 4 n , an upper bound on thenumber of generated SAT instances. The cost of stage 2 is difficult to specify rigorouslysince it may depend on the specific SAT solver that is used. However, there are atmost 4 n / possible assignments for each instance, so using a solver that checks eachassignment once gives an overall cost of 4 n / . The running time is better in practice,though seemingly still exponential.Finally, we provide counts of the total number of complex Golay pairs of length n ≤
28 in Table 3. The sizes of Ω seqs match those given by Fiedler (2013) in all cases,the sizes of Ω all match those given by Gibson and Jedwab (2011) for all n ≤
26 (thelargest length they included) and the sizes of Ω inequiv match those given by Craigen et al.(2002) for all n ≤
19 (the largest length they exhaustively solved).Because Fiedler (2013), Gibson and Jedwab (2011), and Craigen et al. (2002)do not provide implementations or timings for the enumerations they completed it isnot possible for us to compare the efficiency of our algorithm to previous algorithms.However, in Table 4 we compare our implementation’s timings to the timings wepreviously presented (Bright et al., 2018b). Table 4 shows that the improved version ofour algorithm performs about an order of magnitude faster in general.
5. Future Work
Besides increasing the length to which complex Golay pairs have been enumeratedthere are a number of avenues for improvements that could be made in future work. Asone example, we remark that we have not exploited the algebraic structure of complexGolay pairs revealed by Craigen, Holzmann, and Kharaghani (2002). In particular,their work contains a theorem implying that if p ≡ n and A is a member of a complex Golay pair of length n then A ( z ) is not irreducibleover F p ( i ). Ensuring that this property holds could be added to the filtering conditionsused in stage 1. In fact, the authors relate the factorization of A ( z ) over F p ( i ) to thefactorization of B ( z ) over F p ( i ) for any complex Golay pair ( A , B ). This factorizationcould potentially be used to perform stage 2 more efficiently, possibly supplementing orreplacing the SAT solver entirely, though it is unclear if such a method would performbetter than our method in practice. In any case, it would not be possible to apply theirtheorem in all lengths since the length n might only be divisible by primes p = p ≡ f ( θ ) defined inSection 3.1 to help find the maximum of f ( θ ). For example, a consequence of Lemma 6and Euler’s identity is that f ( θ ) = N A (cid:48) (0) + n − (cid:88) s = Re( N A (cid:48) ( s )) cos( s θ ) + n − (cid:88) s = Im( N A (cid:48) ( s )) sin( s θ )19 | L even | | L odd | | L A | −
12 3 1 33 3 1 14 3 4 35 12 4 56 12 16 147 39 16 128 48 64 369 153 64 4410 153 204 11811 561 252 9912 645 860 44513 2121 884 27914 2463 3284 29415 8340 3572 165016 9087 12116 82917 31275 12824 323318 34560 46080 1115919 117597 50944 1091820 130215 173620 2687621 446052 194004 8194122 500478 667304 9016323 1694871 732232 11874724 1886562 2515416 20013825 6447250 2727452 70958426 7183879 9578506 73789127 24426370 10591928 761847428 27265578 36354113 3687209
Table 2: The number of sequences A even , A odd , and A that passed the filtering conditions of our algorithm inlengths up to 28. | Ω seqs | | Ω all | | Ω inequiv | Table 3: The number complex Golay pairs in lengths up to 28. The table counts the number of individualsequences, the number of pairs, and the number of pairs up to equivalence. n ISSAC’18 This paper17 0.07 0.0018 0.27 0.0219 0.26 0.0220 0.80 0.1221 3.86 0.8422 10.77 1.0323 45.18 4.8724 86.97 8.3425 702.39 99.2626 − − − Table 4: A comparison between the method presented at ISSAC 2018 (Bright et al., 2018b) and the methodused in this paper. The times measure the total amount of computation time that was used to run a completesearch in the lengths 17 ≤ n ≤
28 when run on the same hardware. showing that f has a simple form as a sum of sines and cosines. Finding the real rootsof the derivative of f would reveal the locations of the local maxima and minima of f .However, again it is unclear if this method would perform better than our approximationmethod in practice. Note that our method is not guaranteed to find the global maximumof f since we only apply the quadratic interpolation optimization method to ranges thatwe know contain a local maximum by examining the initial sampling of f . With a morecareful branch-and-bound algorithm it would be possible to guarantee that the globalmaximum was found but the approximations that we found were effective enough thatwe did not pursue this.Another possible improvement could be obtained by deriving further propertieslike Proposition 19 that complex Golay pairs must satisfy. For example, consider thefollowing property that could be viewed as a strengthening of Proposition 19: a k a n − k − = ( − n + b k b n − k − for k = . . . , n − . An examination of all complex Golay pairs up to length 28 reveals that they all satisfythis property except for a single complex Golay pair up to equivalence. The only pairsthat don’t satisfy this property are equivalent to (cid:0) [1 , , , − , , , − , , [1 , i , i , − , , − i , − i , − (cid:1) and were already singled out by Fiedler, Jedwab, and Parker (2008b) for being theonly known examples of what they call “cross-over” Golay sequence pairs. Since acounterexample exists to this property there is no hope of proving it in general, butperhaps a suitable generalization could be proven.Lastly, the running time analysis of our algorithm could be improved. We havegiven a crude upper bound of 4 n / = n for the complexity of the parts of our algorithm22hat do not run in polynomial time. However, this running time would be improved bytighter bounds on either the size of L A or the number of possible assignments that theSAT solver will try. Acknowledgements
This work was made possible by the facilities of the Shared Hierarchical AcademicResearch Computing Network (SHARCNET) and Compute/Calcul Canada. The authorswould also like to thank the anonymous reviewers for their comments on this paper andon our ISSAC 2018 submission. Their careful reviews improved our work’s clarity andpresentation.
References
Ábrahám, E., 2015. Building bridges between symbolic computation and satisfiabilitychecking. In: Proceedings of the 2015 ACM on International Symposium on Symbolicand Algebraic Computation. ACM, New York, pp. 1–6.URL https://doi.org/10.1145/2755996.2756636
Ábrahám, E., Abbott, J., Becker, B., Bigatti, A. M., Brain, M., Buchberger, B., Cimatti,A., Davenport, J. H., England, M., Fontaine, P., Forrest, S., Griggio, A., Kroening,D., Seiler, W. M., Sturm, T., 2016. SC : Satisfiability Checking meets SymbolicComputation. In: Intelligent Computer Mathematics: 9th International Conference,CICM 2016, Bialystok, Poland, July 25–29, 2016, Proceedings. Springer InternationalPublishing, Cham, pp. 28–43, .URL https://doi.org/10.1007/978-3-319-42547-4_3 Barrett, C., Fontaine, P., Tinelli, C., 2016. The Satisfiability Modulo Theories Library(SMT-LIB).URL
Bright, C., 2017. Computational methods for combinatorial and number theoreticproblems. Ph.D. thesis, University of Waterloo.URL http://hdl.handle.net/10012/11761
Bright, C., Ganesh, V., Heinle, A., Kotsireas, I. S., Nejati, S., Czarnecki, K., 2016.M
ATH C HECK
2: A SAT+CAS Verifier for Combinatorial Conjectures. In: Com-puter Algebra in Scientific Computing - 18th International Workshop, CASC 2016,Bucharest, Romania, September 19–23, 2016, Proceedings. pp. 117–133.URL https://doi.org/10.1007/978-3-319-45641-6_9
Bright, C., Kotsireas, I., Ganesh, V., 2018a. A SAT+CAS method for enumeratingWilliamson matrices of even order. In: Proceedings of the Thirty-Second AAAIConference on Artificial Intelligence, New Orleans, Louisiana, USA, February 2–7,2018. pp. 6573–6580.URL https://doi.org/10.1145/3208976.3209006
Bright, C., Ðokovi´c, D. Ž., Kotsireas, I., Ganesh, V., 2019. A SAT+CAS approach tofinding good matrices: New examples and counterexamples. In: Proceedings of theThirty-Third AAAI Conference on Artificial Intelligence, Honolulu, Hawaii, USA,January 27 – February 1, 2019. pp. 1435–1442.URL https://doi.org/10.1609/aaai.v33i01.33011435
Craigen, R., 1994. Complex Golay sequences. Journal of Combinatorial Mathematicsand Combinatorial Computing 15, 161–169.Craigen, R., Holzmann, W., Kharaghani, H., 2002. Complex Golay sequences: structureand applications. Discrete Math. 252 (1–3), 73–89.URL https://doi.org/10.1016/S0012-365X(01)00162-5
Davis, J. A., Jedwab, J., 1999. Peak-to-mean power control in OFDM, Golay com-plementary sequences, and Reed–Muller codes. IEEE Transactions on InformationTheory 45 (7), 2397–2417.URL https://doi.org/10.1109/18.796380
Donato, P. G., Urena, J., Mazo, M., Alvarez, F., June 2004. Train wheel detection withoutelectronic equipment near the rail line. In: IEEE Intelligent Vehicles Symposium,2004. pp. 876–880.URL https://doi.org/10.1109/IVS.2004.1336500
Fiedler, F., 2013. Small Golay sequences. Advances in Mathematics of Communications7 (4).URL http://doi.org/10.3934/amc.2013.7.379
Fiedler, F., Jedwab, J., Parker, M. G., 2008a. A framework for the construction of Golaysequences. IEEE Transactions on Information Theory 54 (7), 3114–3129.URL https://doi.org/10.1109/TIT.2008.924667
Fiedler, F., Jedwab, J., Parker, M. G., 2008b. A multi-dimensional approach to theconstruction and enumeration of Golay complementary sequences. Journal of Combi-natorial Theory, Series A 115 (5), 753–776.URL https://doi.org/10.1016/j.jcta.2007.10.001
Frank, R., 1980. Polyphase complementary codes. IEEE Transactions on Informationtheory 26 (6), 641–647.URL https://doi.org/10.1109/TIT.1980.1056272
Frigo, M., Johnson, S. G., 2005. The design and implementation of FFTW3. Proceedingsof the IEEE 93 (2), 216–231.URL https://doi.org/10.1109/JPROC.2004.840301 https://doi.org/10.1007/978-3-642-31612-8_12
Gersho, A., Gopinath, B., Odlyzko, A., 1979. Coefficient inaccuracy in transversalfiltering. Bell System Technical Journal 58 (10), 2301–2316.URL https://doi.org/10.1002/j.1538-7305.1979.tb02968.x
Gibson, R. G., Jedwab, J., 2011. Quaternary Golay sequence pairs I: Even length.Designs, Codes and Cryptography 59 (1-3), 131–146.URL https://doi.org/10.1007/s10623-010-9471-z
Golay, M., 1961. Complementary series. IRE Transactions on Information Theory 7 (2),82–87.URL https://doi.org/10.1109/TIT.1961.1057620
Golay, M. J., 1949. Multi-slit spectrometry. JOSA 39 (6), 437–444.URL https://doi.org/10.1364/JOSA.39.000437
Holzmann, W. H., Kharaghani, H., 1994. A computer search for complex Golay se-quences. Australasian Journal of Combinatorics 10, 251–258.Kotsireas, I. S., 2013. Algorithms and metaheuristics for combinatorial matrices. In:Handbook of Combinatorial Optimization. Springer, pp. 283–309.URL https://doi.org/10.1007/978-1-4419-7997-1_13
Kotsireas, I. S., Koukouvinos, C., Seberry, J., 2009. Weighing matrices and stringsorting. Annals of Combinatorics 13 (3), 305–313.URL https://doi.org/10.1007/s00026-009-0027-8
Krone, S., Sarwate, D., 1984. Quadriphase sequences for spread-spectrum multiple-access communication. IEEE Transactions on Information Theory 30 (3), 520–529.URL https://doi.org/10.1109/TIT.1984.1056913
Li, S., Chen, J., Zhang, L., Zhou, Y., 2008. Construction of quadri-phase complete com-plementary pairs applied in MIMO radar systems. In: ICSP 2008: 9th InternationalConference on Signal Processing. IEEE, pp. 2298–2301.URL https://doi.org/10.1109/ICOSP.2008.4697608
Li, Y., Chu, W. B., 2005. More golay sequences. IEEE Transactions on InformationTheory 51 (3), 1141–1145.URL https://doi.org/10.1109/TIT.2004.842775
Liang, J. H., Poupart, P., Czarnecki, K., Ganesh, V., 2017. An empirical study of branch-ing heuristics through the lens of global learning rate. In: International Conferenceon Theory and Applications of Satisfiability Testing. Springer, pp. 119–135.URL https://doi.org/10.1007/978-3-319-66263-3_8
Monagan, M. B., Geddes, K. O., Heal, K. M., Labahn, G., Vorkoetter, S. M., McCarron,J., DeMarco, P., 2005. Maple 10 Programming Guide. Maplesoft, Waterloo ON,Canada.Nazarathy, M., Newton, S. A., Giffard, R., Moberly, D., Sischka, F., Trutna, W., Foster,S., 1989. Real-time long range complementary correlation optical time domainreflectometer. Journal of Lightwave Technology 7 (1), 24–38.URL https://doi.org/10.1109/50.17729
Nowicki, A., Secomski, W., Litniewski, J., Trots, I., Lewin, P., 2003. On the applica-tion of signal compression using Golay’s codes sequences in ultrasound diagnostic.Archives of Acoustics 28 (4).URL http://acoustics.ippt.gov.pl/index.php/aa/article/view/460
Popovi´c, B. M., 1991. Synthesis of power efficient multitone signals with flat amplitudespectrum. IEEE transactions on Communications 39 (7), 1031–1033.URL https://doi.org/10.1109/26.87205
Riel, J., 2006.
NSOKS : A M
APLE script for writing n as a sum of k squares.URL Seberry, J. R., Wysocki, B. J., Wysocki, T. A., 2002. On a use of Golay sequencesfor asynchronous DS CDMA applications. In: Advanced Signal Processing forCommunication Systems. Springer, pp. 183–196.URL https://doi.org/10.1007/0-306-47791-2_14
Sivaswamy, R., 1978. Multiphase complementary codes. IEEE Transactions on Infor-mation Theory 24 (5), 546–552.URL https://doi.org/10.1109/TIT.1978.1055936
Sun, W., Yuan, Y., 2006. Optimization theory and methods: nonlinear programming.Springer Science & Business Media.URL https://doi.org/10.1007/b106451
Tseng, C. C., 1971. Signal multiplexing in surface-wave delay lines using orthogonalpairs of Golay’s complementary sequences. IEEE Transactions on Sonics and Ultra-sonics 18 (2), 103–107.URL https://doi.org/10.1109/T-SU.1971.29599
Zulkoski, E., Bright, C., Heinle, A., Kotsireas, I. S., Czarnecki, K., Ganesh, V., 2017.Combining SAT solvers with computer algebra systems to verify combinatorialconjectures. Journal of Automated Reasoning 58 (3), 313–339.URL https://doi.org/10.1007/s10817-016-9396-yhttps://doi.org/10.1007/s10817-016-9396-y