Complex Golay Pairs up to Length 28: A Search via Computer Algebra and Programmatic SAT
Curtis Bright Ilias Kotsireas Albert Heinle Vijay Ganesh
cs.SC cs.LO math.CO






uw]Curtis Bright
wlu]Ilias Kotsireas
uw]Albert Heinle
uw]Vijay Ganesh
[uw]University of Waterloo
[wlu]Wilfrid Laurier University


We use techniques from the fields of computer algebra and satisfiability checking
to develop a new algorithm to search for complex Golay pairs.
We implement this algorithm and use it to
perform a complete search for complex Golay pairs of lengths up to 28.
In doing so, we find that complex Golay pairs exist in the lengths 24 and 26 but
do not exist in the lengths 23, 25, 27, and 28.
This independently verifies work done by F. Fiedler in 2013 
and confirms the 2002 conjecture of Craigen, Holzmann, and Kharaghani 
that complex Golay pairs of length 23 don't exist.
Our algorithm is based on the recently proposed SAT+CAS paradigm 
of combining SAT solvers with computer algebra systems to efficiently search large
spaces specified by both algebraic and logical constraints.
The algorithm has two stages:
first, a fine-tuned computer program uses functionality
from computer algebra systems and numerical libraries to construct a list containing every sequence
that could appear as the first sequence in a complex Golay pair 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 complex Golay pair.
This extends work originally presented at the
International Symposium on Symbolic and Algebraic Computation (ISSAC) in 2018; 
we discuss and implement several improvements to our algorithm that enabled us to improve the
efficiency of the search and increase the maximum length
we search from length 25 to 28.



Complex Golay pairs; Boolean satisfiability; SAT solvers; Exhaustive search; Autocorrelation






§ INTRODUCTION


The sequences that are now referred to Golay sequences or Golay pairs were first 
introduced by <cit.> in a groundbreaking paper on multislit spectrometry.
Later, <cit.> made a detailed study of their elegant mathematical
properties and in this paper he referred to them as complementary series:

Regardless of past or possible future applications, the writer has found these complementary
series mathematically appealing…
Since then, Golay pairs and their generalizations have been widely studied
and applied to a huge and surprisingly varied number of problems in engineering.
For example, they have been applied to
encoding acoustic surface waves <cit.>, 
spread-spectrum systems <cit.>, 
optical time domain reflectometry <cit.>, 
CDMA networks <cit.>, 
medical ultrasounds <cit.>, 
train wheel detection <cit.>, 
radar systems <cit.>, 
and wireless networks <cit.>. 
Golay pairs consist of two sequences and 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 <ref> in Section <ref> for
the formal definition.

Although Golay defined his complementary series over an alphabet of ±1, later authors have generalized
the alphabet to include nonreal roots of unity such as 
the fourth root of unity i√(-1).
In this paper, we focus on the case where the alphabet
is ±1,± i. 
In this case the resulting sequence pairs
are sometimes referred to as quadriphase, 4-phase, 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 complex Golay number.


Complex Golay pairs have been extensively studied by many authors
in several different contexts.
Initially they were studied in the context of signal processing
where some researchers ran searches for “polyphase complementary codes”
in the process of studying signal encoding methods.
In particular,
<cit.> found
a complex Golay pair of length 3 and <cit.>
found complex Golay pairs of length 5 and 13
and stated that 7, 9, 11, 15, and 17 were not complex Golay numbers.
Complex Golay pairs were also introduced by <cit.>
(independently of the above works)
in order 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·5^b·13^c
for integers a,b,c≥0.
Craigen also proved that if m and n are complex Golay numbers
then 2mn is a complex Golay number.

An exhaustive search for complex Golay pairs up to length 13 was performed by <cit.>.
This search verified the remark of <cit.>
that 7 and 9 are not complex Golay numbers but
found that 11 is in fact a complex Golay number
(the cause of the mistaken claim is unknown
since it was based on computer programs that were never published).
Later, an exhaustive search up to length 19 was performed by <cit.>,
showing that 14, 15, 17, and 19 are not complex Golay numbers.
In addition,
they reported
that 21 was not a complex Golay number
and conjectured that 23 was not a complex Golay number
based on a partial search for complex Golay pairs of length 23.

Another line of research is to provide explicit
constructions for complex Golay pairs and a number of results of this
form have been published as well.
<cit.> provide a construction that explains
the existence of all known complex Golay pairs whose lengths are a power
of 2,
including complex Golay pairs of length 16 discovered by <cit.> to
not fit into a construction given by <cit.>.
<cit.> provide
a construction that explains the existence of all complex Golay pairs
up to length 26 and give a table that lists the total number of complex Golay pairs up to length 26.
This table was produced by the mathematician Frank Fiedler,
who describes his enumeration method in a subsequent paper <cit.>
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's table giving the total number of complex Golay pairs
up to length 28.
We implemented our method 
and obtained counts up to length 28 after
about 8.5 months of CPU time.
The counts we obtain match those in Fiedler's table in each case, increasing the
confidence that the enumeration was performed without error.
In addition, we also provide counts
for the total number of complex Golay pairs up to well-known equivalence operations
and explicitly publish the sequences
on our website
https://uwaterloo.ca/mathcheckuwaterloo.ca/mathcheck. 
To our knowledge, this is the first time 
that explicit complex Golay pairs (and their counts up to equivalence) have been
published for lengths larger than 19.
Lastly, we publicly release our code
for enumerating complex Golay pairs
so that others may verify and reproduce our work; we were not able to find
any other code for enumerating complex Golay pairs that was publicly available.

Our result is of interest not only because of the verification we provide
but also because of the method we use to perform the verification.
The method proceeds in two stages.  In the first stage,
a fine-tuned computer program performs an exhaustive search
among all sequences that could possibly appear as the first
sequence in a complex Golay pair of a given length
(up to the equivalence defined in Proposition <ref> of Section <ref>).
Several filtering theorems that we
describe in Section <ref> allow us to discard almost all
sequences from consideration.  To apply these filtering theorems we use
functionality from the computer algebra system Maple <cit.> and the
mathematical library FFTW <cit.> as we describe in Section <ref>.
After this filtering is completed we have
a list of sequences of a manageable size such that 
the first sequence of 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 MapleSAT
developed by <cit.>
to determine the sequences from the first stage (if any) that can be
paired up with another sequence to form a complex Golay pair.  A
programmatic SAT
instance is constructed from each sequence found in the first stage
such that the instance is 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 determines the second sequence of a complex Golay pair.

The “SAT+CAS” method that 
we use is of interest in its own right because it links the two previously
separated fields of symbolic computation and satisfiability checking.
Recently there has been interest in combining methods from both fields
to solve computational problems as outlined in the
invited talk at ISSAC by <cit.> and demonstrated by the
SC^2 (satisfiability checking + symbolic computation) project initiated by <cit.>.  Our work fits into
this paradigm and to our knowledge is the first application of a SAT solver
to search for complex Golay pairs,
though previous work exists that uses
a SAT solver to search for other types of
complementary sequences
like those defining Williamson matrices
<cit.>.

A preliminary version of our algorithm and results
was presented at the conference ISSAC <cit.>.
We have since made several improvements to our algorithm
that 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 the timings with those of our previous algorithm
in Section <ref>.
The new timings are about an order of magnitude faster
on the same hardware,
demonstrating the effectiveness of our improvements.
In particular, our new algorithm incorporates the following
updates:



  
More efficient filtering
The preprocessing stage of our algorithm
uses a filtering condition (Corollary <ref> of Section <ref>)
to remove a large number of sequences from consideration.
In our original algorithm we applied this condition using
a large number of equally-spaced points z along the unit circle.
In our new algorithm we apply this condition 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 are the most likely to work.



  
More efficient joining
The first stage of our algorithm requires searching
through two lists and finding elements of the first that
can be joined with elements of the second.  In our
previous algorithm this was done in a totally
brute-force manner, i.e., every item in the 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 the
elements that can be joined via a number of linear
scans through the sorted lists (see Section <ref>),
thereby eliminating
the need to consider all possible pairs.



  
Increased filtering
Because the improved filtering method outlined above
is so much faster than the previous filtering method
we were able to add another round of filtering just
before generating the SAT instances (see Section <ref>).
The result is that approximately 40% of the SAT instances
generated by our previous algorithm are now quickly filtered
before even calling a SAT solver. 


Additionally we include more background details
in this version of the paper.  For example,
we add Section <ref> that gives an alternative
definition of complex Golay pairs that is often useful
and we show that this definition is equivalent with the first definition.



§ BACKGROUND ON COMPLEX GOLAY PAIRS


In this section we present the background necessary to describe our
method for enumerating complex Golay pairs.
First, we require some preliminary definitions
to define what complex Golay pairs are.
We let z denote the complex conjugate of z
(which is just the multiplicative inverse of z
for z on the unit circle) and if A is a sequence we let
A denote the sequence containing the conjugates
of A.


The nonperiodic autocorrelation function
of a sequence A = [a_0, …, a_n-1] of length n is 

    N_A(s) ∑_k=0^n-s-1 a_ka_k+s    for     s=0,…, n-1 .




  A pair of sequences (A,B) with A and B in {±1,± i}^n
  are called a complex Golay pair
  if the sum of their nonperiodic autocorrelations is a constant zero for s≠0, i.e.,
  
    N_A(s) + N_B(s) = 0     for     s = 1, …, n-1 .



Note that if A and B are in {±1,± i}^n (as we assume throughout this paper)
then N_A(0)+N_B(0)=2n by the definition of
the nonperiodic autocorrelation function
and the fact that
zz=1 if z is ±1 or ± i, explaining why s≠0 in Definition <ref>.

([1,1,-1], [1, i, 1]) is a complex Golay pair of length 3 since the
first sequence A has autocorrelations N_A(1)=0 and N_A(2)=-1
and the second sequence B has autocorrelations N_B(1)=0 and N_B(2)=1.




 §.§ Alternative definition


Instead of viewing complex Golay pairs as pairs of sequences
it is also possible to view them as pairs of polynomials.
If A is the sequence [a_0,…,a_n-1] we let A(z) denote the
Hall polynomial a_0+a_1z+…+a_n-1z^n-1 which can be
viewed as the finite generating function of A.
This leads to the following alternative definition of complex Golay pairs.


A pair of polynomials (A(z),B(z)) that have degrees n-1 and coefficients in {±1,± i}
are called a complex Golay pair if A(z)^2 + B(z)^2 = 2n
for all z on the unit circle.



(1+z-z^2,1+iz+z^2) is a complex Golay pair of length 3 since

    []1+z-z^2^2 + []1+iz+z^2^2 = 6

for all z on the unit circle.


These two definitions can be seen to be equivalent using the following lemma
that expresses the squared absolute values of a polynomial
evaluated on the unit circle in terms of the autocorrelations of the sequence
formed by the polynomial's coefficients.


Let A be a complex sequence of length n.  Then

    A(z)^2 = N_A(0) + 2[]∑_s=1^n-1N_A(s)z^-s.

for all z on the unit circle.



Since z is on the unit circle we have z=z^-1, so A(z)^2=A(z)A(z)=A(z)A(z^-1).
Expanding this, we obtain that A(z)^2 is equal to

    ∑_k=0^n-1∑_l=0^n-1 a_k a_l z^k-l .

Collecting terms of z^s together for s from -n+1 to n-1 this becomes

    ∑_s=0^n-1[]∑_j=s^n-1a_ja_j-sz^s + ∑_s=1^n-1[]∑_j=s^n-1a_j-sa_jz^-s

so the coefficient of z^s for positive s is N_A(s) and for negative s is N_A(s).
Using the fact that N_A(s)z^s=N_A(s)z^-s 
this becomes

    N_A(0) + ∑_s=1^n-1[]N_A(s)z^-s + N_A(s)z^-s .

The desired result now follows from the fact that z+z=2(z) for all complex z.


For completeness we now demonstrate the equivalence of Definitions <ref> and <ref>.


A pair (A,B) is a complex Golay pair in the
sense of Definition <ref> if and only if
(A(z),B(z)) is a complex Golay pair in the sense
of Definition <ref>.



If (A,B) is a complex Golay pair in the
sense of Definition <ref> then
by Lemma <ref>

    A(z)^2 + B(z)^2 = N_A(0) + N_B(0) + 2[]∑_s=1^n-1(N_A(s)+N_B(z))z^-s = 2n

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 <ref> then
by Lemma <ref> and subtracting off N_A(0) + N_B(0) = 2n we have

    ∑_s=1^n-1[]N_A(s)z^-s+N_A(s)z^s+N_B(z)z^-s+N_B(s)z^s = 0

for all z on the unit circle.  Since any nonzero rational function has a finite number of zeros in 
the function given on the left hand side must be identically zero.  
In particular, the coefficients of z^-s for s=1, …, n-1 must be zero
and we derive N_A(s)+N_B(s)=0 for s=1, …, n-1 as required.




 §.§ Equivalence operations


There are certain invertible operations that preserve the property of
being a complex Golay pair when applied to a sequence pair (A,B).
These are summarized in the following proposition.



Let ([a_0,…,a_n-1], [b_0,…,b_n-1]) be a complex
Golay pair. Then the following are also complex Golay pairs:


  E1. (Reversal) ([a_n-1,…,a_0], [b_n-1,…,b_0]).

  E2. (Conjugate Reverse A) ([a_n-1,…,a_0], [b_0,…,b_n-1]).

  E3. (Swap) ([b_0,…,b_n-1], [a_0,…,a_n-1]). 

  E4. (Scale A) ([ia_0,…,ia_n-1], [b_0,…,b_n-1]). 

  E5. (Positional Scaling) (i⋆ A, i⋆ B)
where c⋆ A denotes the sequence of coefficients of the polynomial A(cz), i.e.,
[a_0,ca_1,c^2a_2,…,c^n-1a_n-1].



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-1A(z^-1) and z^n-1=1.  Then

    z^n-1A(z^-1)^2 + z^n-1B(z^-1)^2 = A(z^-1)^2 + B(z^-1)^2 = 2n

since (A,B) is a complex Golay pair and z^-1 is on the unit circle.

  E2. The conjugate reverse of A(z) is z^n-1A(z^-1) and z^n-1A(z^-1)=A(z)=A(z).

  E3. B(z)^2+A(z)^2=A(z)^2+B(z)^2=2n.

  E4. A(z) scaled by i is iA(z) and iA(z)=A(z).

  E5. A(iz)^2+B(iz)^2=2n since (A,B) is a complex Golay pair and iz is on the unit circle.




We call two complex Golay pairs (A, B) and (A', B') equivalent
if (A', B') can be obtained from
(A, B) using the transformations described in Proposition <ref>.


The next lemma provides some normalization conditions that can be used when searching
for complex Golay pairs up to equivalence.
Since all complex Golay pairs (A',B') which are equivalent to
a complex Golay pair (A,B) can easily
be generated from (A,B), it suffices to search for complex Golay pairs up to equivalence.


  Let (A', B') be a complex Golay pair.
  Then (A', B') is equivalent to a complex Golay pair
  (A, B) with a_0=a_1=b_0=1 and a_2∈±1, i.


We will transform a given complex Golay sequence pair (A', B') into 
an equivalent normalized one
using the equivalence operations of
Proposition <ref>.
To start with, let A A' and B B'.

First, we ensure that a_0=1.  To do this, we apply operation E4 (scale A) enough times
until a_0=1.

Second, we ensure that a_1=1.  To do this, we apply operation E5 (positional scaling) enough times
until a_1=1; note that E5 does not change a_0.

Third, we ensure that a_2≠ -i.  If it is, we apply operation E1 (reversal) and E2 (conjugate reverse A)
which has the effect of keeping a_0=a_1=1 and setting a_2=i.

Last, we ensure that b_0=1.  To do this, we apply operation E3 (swap) and then operation E4 (scale A) enough times
so that a_0=1 and then operation E3 (swap) again.  This has the effect of not changing A but setting b_0=1.




 §.§ Filtering properties


We now prove some useful properties that all complex
Golay pairs satisfy and that will be exploited by our algorithm
for enumerating complex Golay pairs.
The properties will be used as filtering criteria:
if a sequence does not satisfy them then we know it cannot
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.



Let A be a complex sequence of length n.
If A(z)^2 > 2n for some z on the unit circle
then A is not a member of a complex Golay pair.


Suppose the sequence A was a member of a complex Golay pair whose other
member was the sequence B.
Since B(z)^2≥0, we must have A(z)^2+B(z)^2>2n,
in contradiction to Definition <ref>.



The sequence A[1,1,1] cannot be a member of a complex Golay pair
since A(1)^2=9 is larger than 2n=6.


<cit.> derives a variation of Lemma <ref>
that is useful because it can be applied knowing only around half of the entries of A.
It is derived using the following variation of Lemma <ref>.
Let  be identical to A with the entries of odd index replaced by zeros and let  be
identical to A with the entries of even index replaced by zeros.



Let A be a complex sequence of length n.  Then

    (z)^2 + (z)^2 = N_A(0) + 2[]∑_s=1
    s even^n-1N_A(s)z^-s.




The proof proceeds like in the proof of Lemma <ref> and we derive

    (z)^2 + (z)^2 = ∑_k,l=0
    k,l even^n-1 a_k a_l z^k-l + ∑_k,l=0
    k,l odd^n-1 a_k a_l z^k-l = ∑_k,l=0
    k-l even^n-1 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 for even s
is the same as in Lemma <ref> from which the result follows.



Let A be a complex sequence of length n.
If (z)^2 or (z)^2 is strictly larger than 2n
for some z on the unit circle
then A is not a member of a complex Golay pair.



If (A,B) is a complex Golay pair then by Lemma <ref>
we have

    (z)^2 + (z)^2 + (z)^2 + (z)^2 = 2n

but if either (z)^2>2n
or (z)^2>2n then this cannot hold.



The sequence A[1,x,1,y,1,z,1] cannot be a member of a complex Golay pair
(regardless of the values of x, y, and z)
since (1)^2=16 is larger than 2n=14.




 §.§ Sum-of-squares decomposition types


The next lemma is useful because it allows us to write 2n
as the sum of four integer squares.  It is stated
by <cit.> using a different notation; we use the notation
(A) and (A) to represent the real and imaginary parts of
the sum of the entries of A.
For example, if A[1,i,-i,i] then (A)=(A)=1. 


If (A, B) is a complex Golay pair then

    (A)^2 + (A)^2 + (B)^2 + (B)^2 = 2n .



Using Definition <ref> with z=1 we have

    (A)+(A)i^2 + (B)+(B)i^2 = 2n .

Since (X)+(X)i^2 = (X)^2 + (X)^2 the result follows.


A consequence of Lemma <ref> is
that every complex Golay pair generates a decomposition of 2n into a sum of four integer squares.
In fact, it typically generates several decompositions of 2n into a sum of four squares.
Recall that i⋆ A denotes positional scaling by i (operation E5) on the sequence A.
If (A,B) is a complex Golay pair then applying operation E5 to this pair k times 
shows that (i^k⋆ A,i^k⋆ B) 
is also a complex Golay pair.  
By using Lemma <ref> on these complex Golay pairs
one obtains the fact that 2n can be decomposed as the sum of 
four integer squares as

    i^k⋆ A^2 + i^k⋆ A^2 + i^k⋆ B^2 + i^k⋆ B^2 .

For k>3 this produces no new decompositions but in general
for k=0, 1, 2, and 3 this produces four distinct
decompositions of 2n into a sum of four squares.

With the help of a computer algebra system (CAS) one can
enumerate every possible way that 2n may be written as a sum of four integer squares.  For example,
when n=23 one has 0^2+1^2+3^2+6^2=2·23 and 1^2+2^2+4^2+5^2=2·23 as well as all
permutations of the squares and negations of the integers being squared.
During the first stage of our enumeration method only the first sequence of a complex Golay pair
is known, so at that stage we cannot compute its whole
sums-of-squares decomposition.
However, it is still possible to filter some sequences from consideration based on analyzing
the two known terms in the sums-of-squares decomposition.

For example, say that A is the
first sequence in a potential complex Golay pair of length 23 with (A)=0 and (A)=5.
We can immediately discard A from consideration because 
there is no way to chose the  and  of B to complete the sums-of-squares
decomposition of 2n, i.e., there are no integer solutions (x,y) of 0^2+5^2+x^2+y^2=2n.



§ ENUMERATION METHOD


In this section we describe in detail the method we used to
perform a complete enumeration of all complex Golay pairs up to length 28.
Given a length n our goal is to find all {±1,± i} sequences A and B
of length n such that (A,B) is a complex Golay pair.



 §.§ Preprocessing: Enumerate possibilities for  and 


The first step of our method uses Fiedler's trick of considering
the entries of A of even index separately from the entries of A of odd index.
There are approximately n/2 nonzero entries in each of  and 
and there are four possible values for each nonzero entry.  Therefore there
are approximately 2·4^n/2=2^n+1 possible sequences to check in this
step.  Additionally, by Lemma <ref> we may assume the first nonzero
entry of both  and  is 1 and that the second nonzero
entry of  is not -i, decreasing the number of sequences
to check in this step by more than a factor of 4.
It is quite feasible to perform a brute-force search
through all such sequences when n≈30.

We apply Corollary <ref> to every possibility for  and .
There are an infinite number of z on the unit circle so
it is not possible to apply Corollary <ref> using all such z.
One simple approach 
is to try a sufficiently large number of points z so that 
when a point exists with A'(z)^2>2n
(where A' is either  or ) such a point is usually discovered.
For example, <cit.> tested Corollary <ref>
for 2^14 equally-spaced points z around the unit circle.

Instead, in our implementation we use a nonlinear programming method
to estimate the maximum of A'(z)^2 for z on the unit circle.
This can be done with the NLPSolve command of the
computer algebra system Maple
though for efficiency we use a custom C implementation
of a variant of the quadratic interpolation
method described by <cit.>.

We write A'(z)^2 in terms of the real variable
θ by using the substitution z=e^iθ and define f(θ)A'(e^iθ)^2.
The quadratic interpolation method to estimate the maximum of f(θ) over 0≤θ<2π proceeds as follows:


  * Evaluate f at the 2^7 points θ_0, …, θ_127 where θ_k2π k/128
and let f_k f(θ_k).
If f_k>2n for some k we can immediately filter A' by Corollary <ref>.

  * For every value of f that is larger than its neighbours (i.e.,
values of f_k that satisfy f_k-1≤ f_k and f_k≥ f_k+1)
we use interpolation to find the quadratic polynomial that passes through the points
(θ_k-1,f_k-1), (θ_k,f_k), and (θ_k+1,f_k+1).
Let θ^* be the value of θ that maximizes the quadratic polynomial;
this can be computed to be exactly

    1/2·f_k-1(θ_k^2-θ_k+1^2)+f_k(θ_k+1^2-θ_k-1^2)+f_k+1(θ_k-1^2-θ_k^2)/f_k-1(θ_k-θ_k+1)+f_k(θ_k+1-θ_k-1)+f_k+1(θ_k-1-θ_k)

and let f^* f(θ^*).

  * Note that f^* is often a better approximation to a local maximum of f than
f_k was, and if f^* > 2n then we can filter A' by Corollary <ref>.
Otherwise we can use the point (θ^*,f^*) to derive a tighter interval in
which a local maximum of f must lie.
For example, if θ_k<θ^*<θ_k+1 and f^*>f_k then we
can repeat the previous step except using the points (θ_k,f_k), (θ^*,f^*),
and (θ_k+1,f_k+1).


One can use this method to derive more and more accurate approximations
to the local maxima of f though this is mostly unnecessary for our purposes
as we only care about finding a single value for θ with f(θ)>2n.
A good approximation to a local maximum was usually
found after a single interpolation step so in our implementation
we would move on to looking for another local maximum of f
after repeating the interpolation step three times.
If all values of f_k were examined and no
points θ were found with f(θ)>2n then
we save A' as a sequence that could not be filtered.


At the conclusion of this step we have two lists: one list 
of the  that were not discarded and one list  of the 
that were not discarded.



 §.§ Stage 1: Enumerate possibilities for A


We now enumerate all possibilities for A by joining all
possibilities for  with all possibilities for .
The most straightforward way of doing this would be to simply try
all A_1∈ and A_2∈; this was done by <cit.>.
However, because both  and  can contain millions
of sequences it can be inefficient to try all possible pairings
(A_1,A_2).
However, many pairings can be eliminated
by using the conditions implied by Lemma <ref>;
we now show how to find all possible
pairings that satisfy these conditions without
needing to try all A_1∈ and A_2∈.

Let (u_0,u_1,u_2,u_3) be an arbitrary quadruple of integers.
We will describe how to efficiently find all A
whose entries of odd index are in , whose entries of even index
are in , and that satisfy

    (A)=u_0,   (A)=u_1,   (i⋆ A)=u_2,   (i⋆ A)=u_3 . 1

For each A_1∈ we form the vector 

    V_1  ((A_1),(A_1),(i⋆ A_1),(i⋆ A_1))

and for each A_2∈ we form the vector 

    V_2  (u_0-(A_2),u_1-(A_2),u_2-(i⋆ A_2),u_3-(i⋆ A_2)) .

We now show that the A that satisfy relationship (<ref>) are exactly those formed
by the A_1 and A_2 with V_1=V_2.


Let A_1∈ and A_2∈ and let A be the sequence
that satisfies the relationship A(z)=A_1(z)+A_2(z), i.e., A is a sequence
of length n with {±1,± i} entries.
Then A satisfies (<ref>) if and only if the vectors V_1 and V_2
defined as above satisfy V_1=V_2.



Consider the first entry of V_1 and V_2.  If they are equal, we have

    (A_1) = u_0 - (A_2) .

This implies that (A)=u_0 since (A)=(A_1)+(A_2).
Similarly, we derive (A)=u_1, (i⋆ A)=u_2, and (i⋆ A)=u_3,
as required.  Conversely, (A)=u_0 implies that the first entries of V_1
and V_2 are equal and the other equations from (<ref>) imply that
their other entries are also equal.


Therefore we have translated the problem of finding the sequences A that satisfy (<ref>)
into the problem of finding A_1∈ and A_2∈
that have matching vectors V_1 and V_2.
We can efficiently solve this problem using a string sorting
algorithm as described by (for example) <cit.>.
First, we sort the list  so that the
vectors V_1 formed above appear in lexicographically increasing
order when iterating through L_1∈.  Similarly, we sort
 in the same way so that the V_2 appear in lexicographically
increasing order.

We can then find all V_1 that match with V_2 by a linear
scan through the lists L_1 and L_2
and this requires only one pass
through each list.
For example, if at some point we find that V_1 is lexicographically greater
than V_2 then we try the next L_2 in the list  and if
we instead find that V_2 is lexicographically greater than V_1
then we try the next L_1 in the list .
If instead we find V_1=V_2 we iterate through the L_1 and L_2
that share the same value of V_1 and V_2 and save the A
formed by joining all such A_1 and A_2 in a new list .

It remains to determine the appropriate values of (u_0,u_1,u_2,u_3)
to use in (<ref>).
To do this we use a quadratic Diophantine equation solver and find all solutions
of the Diophantine equation

    u^2 + v^2 + x^2 + y^2 = 2n     in integers u, v, x, y .

Let U be the set of all pairs (u,v) for which
this equation is solvable.
By Lemma <ref> and operation E5
we know that if A is a complex Golay pair
then ((A),(A)) and ((i⋆ A),(i⋆ A))
must be members of U.
Therefore, we apply the above joining procedure
for all ((u_0,u_1),(u_2,u_3))∈ U^2.


At the conclusion of this stage we will have a list of sequences 
that could potentially be a member of a complex Golay pair.
By construction, the first member of all complex Golay pairs
(up to the equivalence described in Lemma <ref>)
of length n will be in .
To decrease the size of  even further we perform
additional filtering before adding sequences into ,
for example, we ensure that

    [](-1⋆ A),(-1⋆ A)  and  [](-i⋆ A),(-i⋆ A)

are both in U for each A added to .  We also
filter those A with A(z)^2 > 2n
for some z on the unit circle (see Section <ref>
for optimization details).



 §.§ Stage 2: Construct the second sequence B from A


In the second stage we take as input the list  generated
in the first stage, i.e., a list of the sequences A that were not filtered by
any of the filtering theorems we applied.  For each A∈ we attempt to
construct a second sequence B such that (A,B) is a complex Golay 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 are unknown and encoded using Boolean variables.
Because there are four possible values for each entry of B
we use two Boolean variables to encode each entry. 
Although the exact encoding used is arbitrary, we fixed the following
encoding in our implementation, where the variables v_2k and v_2k+1
represent b_k, the kth entry of B:

    [   v_2k v_2k+1    b_k;                    1;                   -1;                    i;                   -i ]


To encode the property that (A,B) is a complex Golay pair in our SAT instance
we add the conditions that define (A,B) to be a complex Golay pair, i.e.,

    N_A(s) + N_B(s) = 0     for     s = 1, …, n-1 .

These equations could be encoded using clauses in conjunctive normal form
(for example by constructing
logical circuits to perform complex multiplication and addition and then
converting those circuits into CNF clauses).  However, we found that a much
more efficient and convenient method was to use a programmatic SAT solver.

The concept of a programmatic SAT solver was first introduced
by <cit.> where a programmatic SAT solver was
shown to be more efficient than a standard SAT solver when solving
instances derived from RNA folding problems.
More recently, a programmatic SAT solver was also shown to be
useful when searching for Williamson matrices and good matrices by <cit.>.
Generally, programmatic SAT solvers perform well when there is additional
domain-specific knowledge known about the problem being solved.
Often this knowledge cannot easily be encoded into a SAT instance directly
but can be given to a programmatic SAT solver
to help guide the solver in its search.

Concretely, a programmatic SAT solver is compiled with a piece
of code that encodes a property any solution
must satisfy.  Periodically the SAT solver will run this code
while performing its search, and if the current partial
assignment violates a property that is 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 it is used to increase the efficiency
of the remainder of the search.
The reason that these clauses can be so useful is because they can
encode facts that the SAT solver would have no way of learning
otherwise, since the SAT solver has no knowledge of the domain of
the problem.

Not only does this paradigm allow the SAT solver to perform its
search more efficiently, it also allows instances
to be much more expressive.  Under this framework SAT instances do not
have to consist solely of Boolean formulas in conjunctive normal
form (the typical format of SAT instances) but can consist of
clauses in conjunctive normal form combined with a piece of
code that programmatically expresses clauses.
Increased expressiveness is also a feature of SMT (SAT modulo theories) solvers,
though SMT solvers typically require additional overhead
and only support a fixed number of theories such as those
specified in the SMT library <cit.>.  Additionally,
one can compile instance-specific programmatic SAT solvers
that are tailored to perform searches for a specific class of
problems.

For our purposes we use a programmatic
SAT solver tailored to search for
sequences B that when paired with a given sequence A
form a complex Golay pair.
Each instance will contain the 2n variables v_0, …, v_2n-1
that encode the entries of B as previously specified.
In detail, the code given to the SAT solver does the following:


  * Compute and store the values N_A(k) for k=1, …, n-1.

  * Initialize s to n-1.  This will be a variable that
controls which autocorrelation condition we are currently examining.

  * Examine the current partial assignment to v_0, v_1,
v_2n-2, and v_2n-1.  If all these values have been assigned
then we can determine the values of b_0 and b_n-1.
From these values we compute
N_B(s) = b_0b_n-1.
If N_A(s)+N_B(s)≠0 then (A,B) cannot be a complex Golay pair
(regardless of the values of b_1, …, b_n-2) and therefore
we learn a conflict clause saying that b_0 and b_n-1 cannot
both be assigned to their current values.  More explicitly,
if v_k^cur represents the literal v_k when v_k
is currently assigned to true and the literal ¬ v_k when
v_k is currently assigned to false
we learn the clause

    ¬ v_0^cur¬ v_1^cur¬ v_2n-2^cur¬ v_2n-1^cur

that says that at least one of {v_0, v_1, v_2n-2, v_2n-1} must have their value changed.

  * 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)≠0 then learn a clause preventing the values
of b_k that appear in the definition of N_B(s) from being assigned
the way that they currently are.  Continue to repeat this step until s=0.

  * If all values of B are assigned but no clauses
have been learned in steps 3–4 then output the complex Golay pair (A,B).
If an exhaustive search is desired, learn a clause preventing
the values of B from being assigned the way they currently are;
otherwise learn nothing and return control to the SAT solver.

For each A in the list  from stage 1 we run a SAT solver with the above
programmatic code. The list of all outputs (A,B) in step 5 shown
above now form a complete list of complex Golay pairs of length n
up to the equivalence given in Lemma <ref>.  In fact,
since Lemma <ref> says that we can set b_0=1 we can
assume that both v_0 and v_1 are always set to false.
In other words, we can add the two unit clauses ¬ v_0
and ¬ v_1 into our SAT instance without omitting any
complex Golay pairs up to equivalence.



 §.§ Postprocessing: Enumerating all complex Golay pairs


At the conclusion of the second stage we have obtained a list of complex
Golay pairs of length n such that every complex Golay pair of length n
is equivalent to some pair in our list.  However, because we have not
accounted for all the equivalences in Section <ref> some pairs
in our list may be equivalent to each other.  In some sense such pairs
should not actually be considered distinct, so to count how many distinct
complex Golay pairs exist in length n we would like
to find and remove pairs that are equivalent from the list.
Additionally, to verify the counts given by <cit.>
it is necessary to produce a list that contains
all complex Golay pairs.  We now describe an algorithm that
does both, i.e., it produces a list of all complex Golay pairs as 
well as a list of all inequivalent complex Golay pairs.

In detail, our algorithm performs the following steps:


  * Initialize Ω_all to be the empty set.  This variable will be a set that will be
populated with and eventually contain all complex Golay pairs of length n.

  * Initialize Ω_inequiv to be the empty set.  This variable
will be a set that will be populated with and eventually contain all inequivalent
complex Golay pairs of length n.

  * For each (A,B) in the set of complex Golay pairs generated in stage 2:


  * If (A,B) is already in Ω_all then skip this (A,B)
and proceed to the next pair.

  * Initialize Γ to be the set containing (A,B).
This variable will be a set that will be populated with and eventually
contain all complex Golay pairs equivalent to (A,B).

  * For every γ in Γ add
E1(γ), …, E5(γ) to Γ.
Continue to do this until every pair in Γ has been examined and no new
pairs are added to Γ.

  * Add (A,B) to Ω_inequiv and add all pairs in Γ to
Ω_all.


Additionally, pseudocode for this algorithm is given in Algorithm <ref>.


After running this algorithm listing the members of Ω_all
gives a list of all complex Golay pairs of length n and
listing the members of Ω_inequiv
gives a list of all inequivalent complex Golay pairs of length n.
At this point we can also construct the complete list of sequences
that appear in any complex Golay pair of length n.
To do this it suffices to add A and B to a new
set Ω_seqs for each (A,B)∈Ω_all.



 §.§ Optimizations


Although the method described will correctly enumerate all complex Golay
pairs of a given length n, for the benefit of potential implementors
we mention a few optimizations that we found helpful.

In the preprocessing step and stage 1 it is necessary to evaluate
a polynomial at points on the unit circle and determine its squared
absolute value.
The fastest way we found to do this used the discrete Fourier transform (DFT).
For example, let A' be the sequence , , or A
under consideration but padded with trailing zeros
so that A' is of length N.
By definition of the discrete Fourier transform we have that
the jth entry of (A')
is exactly A'(z_j) where z_jexp(2π ij/N).
Thus, we determine the values of A'(z_j)^2 by taking
the squared absolute values of the entries of (A').
If A'(z)^2>2n for some z then by Lemma <ref>
or Corollary <ref>
we can discard A' from consideration.  To guard against potential
inaccuracies introduced by the algorithms used to compute the DFT we
actually ensure that A'(z)^2>2n+ϵ for some tolerance ϵ
that is small but larger than the accuracy of the DFT
(e.g., ϵ=10^-3).


In stage 1 we solve the quadratic Diophantine equation
u^2 + v^2 + x^2 + y^2 = 2n in integers u, v, x, y.
In fact, we can also add the constraints

    u+v≡ n 2  and   x+y≡ n2

(the second is unnecessary, being implied by the first) because of the following lemma.

Suppose R and I are the  and  of a sequence
A∈±1,± i^n.  Then R+I≡ n2.


Let #_c denote the number of entries in A with value c.  Then

    R + I = (#_1 - #_-1) + (#_i - #_-i) ≡#_1 + #_-1 + #_i + #_-i2

since -1≡12.
The quantity on the right is 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^2 + I^2 + x^2 + y^2 = 2n
are solvable in integers (x,y) where R and I are given.
This can be 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 fast lookup.
At this point we also check if we can find a z on the unit circle with A(z)^2>2n.
Because the joining process is the bottleneck of our algorithm it is
important to make this filtering step as efficient as possible;
we did this by
computing the values of A(z)^2 via the expression
A_1(z)+A_2(z)^2 and precomputing the values of A_1(z) and A_2(z)
for A_1∈ and A_2∈ 
on points of the form z=exp(2π i j/32) with j=0,…,31.

We found that it was more efficient
to not check the condition for each j in ascending order (i.e., for each z
in ascending complex argument)
but to first perform the check on points z with larger spacing between them.
In our implementation we checked the condition for z of the form exp(2π ij/N)
with j odd and N=8, 16, and 32 in that order.
(This ignores checking the
condition when z∈{±1,± i} but that is desirable since
A(i^k)^2=(i^k⋆ A)^2+(i^k⋆ A)^2
and the sums-of-squares condition is a strictly stronger filtering method.)

The downside of precomputing the values of A_1(z) and A_2(z)
for A_1∈ and A_2∈ is that when the lists 
and  are large this requires a significant amount of
memory.  In our implementation storing the values of A_2(z)
for all A_2∈ in the
lengths 27 and 28 each required about 8GB of RAM, the
searches in lengths 25 and 26 each required about 2GB of RAM,
and the searches in lengths 23 and 24 each required about 0.5GB of RAM.
However, the amount of precomputation could be increased or decreased
depending on the amount of memory available.

In stage 1 we make a pass through the lists  and 
for each valid possibility of the values (u_0,u_1,u_2,u_3).
However, it is possible to reuse some work between passes.  For example,
it is only necessary to sort the lists  and  once.
Even though the values of V_2 depend on the values (u_0,u_1,u_2,u_3)
the relative ordering of  is unaffected since all V_2 are
shifted by the same amount in each pass.  It is even possible to
do some parts of the passes simultaneously: for example, if the
first coordinate of V_1 is greater than the first coordinate of V_2
for some u_0 then V_1 is lexicographically greater than V_2
in all passes with that value of u_0 and so we can iterate to the
next L_2 in all such passes.

At the conclusion of stage 1 we added an additional round of filtering
to the sequences in  that was not in our previous work.
Since this step was not
a bottleneck we used a more involved filtering process here;
in our implementation we ensured that A(z)^2≤2n
for the 2^10 points of the form z_j=exp(2π ij/2^10).  Additionally,
we keep track of the local maxima found by this process and use the
quadratic interpolation method described in Section <ref>
to find more accurate approximations to the local maxima of A(z)^2.

In stage 2 one can also include properties that complex Golay sequences
must satisfy in the SAT instances.
As an example of this, we state the following proposition that
was published by <cit.> at ISSAC 2018.

Let (A,B) be a complex Golay pair.  Then

    a_k a_n-k-1 b_k b_n-k-1 = ± 1     for    k=0, …, n-1.



To prove this, we use the following simple lemma.

Let c_k∈_4 for k=0, …, n-1.  Then

    ∑_k=0^n-1i^c_k=0     implies    ∑_k=0^n-1 c_k≡02.



Let #_c denote the number of c_k with value c.
Note that the sum on the left implies that #_0=#_2 and #_1=#_3
because the 1s must cancel with the -1s and the is must cancel with the -is.
Then
∑_k=0^n-1c_k=#_1+2#_2+3#_3=4#_1+2#_2≡02.

We now prove Proposition <ref>.

Let c_k, d_k∈_4 be such
that a_k=i^c_k and b_k=i^d_k.  Using this notation
the multiplicative equation from Proposition <ref> becomes
the additive congruence

    c_k + c_n-k-1 + d_k + d_n-k-1≡ 0 2 . 2

Since (A,B) is a complex Golay pair, the autocorrelation equations give us

    ∑_k=0^n-s-1[]i^c_k-c_k+s+i^d_k-d_k+s = 0

for s=1, …, n-1.  Using Lemma <ref> and the fact that
-1≡12 gives

    ∑_k=0^n-s-1[]c_k + c_k+s + d_k + d_k+s≡ 0 2

for s=1, …, n-1.
With s=n-1 (or s=1) one immediately derives (<ref>) for k=0.
Adding together these congruences for s=n-1 and s=n-2 
derives (<ref>) for k=1.
The congruences for s=n-2 and s=n-3 give (<ref>) for k=2 and
proceeding in this manner one derives (<ref>) for all k.


In short, Proposition <ref> tells us that an even number of
a_k, a_n-k-1, b_k, and b_n-k-1 are real for each k=0, …, n-1.
For example, if exactly one of a_k and a_n-k-1 is real then exactly
one of b_k and b_n-k-1 must also be real.
In this case, using our encoding from
Section <ref> we add the two binary clauses

    v_2k v_2(n-k-1)    and    ¬ v_2k¬ v_2(n-k-1)

to our SAT instance.
These clauses say that exactly one of v_2k and v_2(n-k-1) is true.
Conversely, if an even number of a_k and a_n-k-1 are real then
an even number of b_k and b_n-k-1 must also be real.  In this case we
add the two binary clauses

    v_2k¬ v_2(n-k-1)    and    ¬ v_2k v_2(n-k-1)

to our SAT instance.



§ RESULTS


In order to provide a verification of the counts given by <cit.>,
<cit.>, and <cit.>
we implemented the enumeration method described in Section <ref>.
The preprocessing step was performed by a C program
and used the mathematical library FFTW by <cit.> for computing
the values of f_0, …, f_127 as described in Section <ref>. 
Stage 1 was performed by a C++ program,
used FFTW for computing the values of A_1(z) and A_2(z),
and a Maple script by <cit.>
for determining the solvability of the Diophantine equations 
given in Section <ref>.  Stage 2 was performed
by the programmatic SAT solver MapleSAT by <cit.>.
The postprocessing step was performed by a Python script.



We ran our implementation on a cluster of machines running CentOS 7 and
using Intel Xeon E5-2683V4 processors running at 2.1 GHz and
using at most 8.5GB of RAM.
The lengths up to 22 were run on a single core,
the lengths 23 and 24 were run on 10 cores,
the lengths 25 and 26 were run on 100 cores, and the lengths
27 and 28 were run on 1000 cores.
The work was parallelized across N cores by partitioning
the list  into N sublists of approximately equal size.
Everything in the stages proceeded exactly as before
except that in stage 1 the list 
was about N times shorter than it would 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 CPU time
but completed in about 3 hours when distributed across 1000 cores.
The timings for the preprocessing step and the two stages of our algorithm are
given in Table <ref>; the timings for the postprocessing step
were negligible.  The times are given as the total amount of CPU time used
across all cores.
Our code is available online as a part of the MathCheck project
(see https://uwaterloo.ca/mathcheckuwaterloo.ca/mathcheck)
and the resulting enumeration of complex Golay pairs 
is available from <https://doi.org/10.5281/zenodo.1246337>.

The sizes of the lists  and  computed in the preprocessing step
and the size of the list  computed in stage 1
are given in Table <ref> for all lengths n≤28.
Without applying any filtering  would have size 4^n so
Table <ref> demonstrates the power of the criteria we used to
perform filtering; for example in length 28 less than 10^-8%
of the possible sequences A were added to .
A result of <cit.> implies that for z on the unit circle
the maximum value of A(z)^2
is nlog n + O(nloglog n)
for almost all ±1,± i-polynomials A of degree n.
In particular, the number of polynomials that do not satisfy
this is at most 4^n/(log n)^4, implying that
the size of  will be o(4^n)
though in practice we even have ≤2^n.

The generated SAT instances had 2n variables
(encoding the entries b_0, …, b_n-1),
two unit clauses (encoding b_0=1),
2n/2 binary clauses (encoding Proposition <ref>), and
n-1 programmatic constraints (encoding Definition <ref>).
Once the values of
b_k are known for k=0, …, n/2,
the programmatic constraints uniquely determine the remaining values of b_k.
Thus 4^n/2 is a crude upper bound on the number
of possible assignments to the values of b_0, …, b_n-1.

We now give an upper bound on the complexity of the parts of our algorithm that do
not run in polynomial time.  The preprocessing and postprocessing costs
are dominated by the costs of stage 1 and 2. The cost of stage 1
is at most 4^n, an upper bound on the number of generated SAT instances.
The cost of stage 2 is difficult to specify rigorously since it may
depend on the specific SAT solver that is used.  However, there are
at most 4^n/2 possible assignments for each instance, so using a solver
that checks each assignment once gives an overall cost of 4^3n/2.
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 <ref>.
The sizes of Ω_seqs match those given
by <cit.> in all cases,
the sizes of Ω_all
match those given by <cit.> for all n≤26
(the largest length they included)
and the sizes of Ω_inequiv
match those given by <cit.>
for all n≤19 (the largest length they exhaustively solved).







Because <cit.>, <cit.>, and <cit.> do not provide
implementations or timings for the enumerations they completed it is not possible
for us to compare the efficiency of our algorithm to previous algorithms.
However, in Table <ref> we compare our implementation's timings
to the timings we previously presented <cit.>.
Table <ref> shows that the improved version of our algorithm performs about
an order of magnitude faster in general.



§ FUTURE WORK


Besides increasing the length to which complex Golay pairs have been enumerated
there are a number of avenues for improvements that could be made in
future work.  As one example, we remark that
we have not exploited the algebraic structure
of complex Golay pairs revealed by <cit.>.
In particular, their work contains a theorem implying that if p≡34
is a prime dividing n and A is a member of
a complex Golay pair of length n then A(z) is not irreducible
over _p(i).  Ensuring that this property holds
could be added to the filtering conditions used in stage 1.
In fact, the authors relate the factorization of A(z) over _p(i)
to the factorization of B(z) over _p(i) for any complex Golay pair
(A,B).  This factorization could potentially be used to perform stage 2
more efficiently, possibly supplementing or replacing the SAT solver entirely,
though it is unclear if such a method would perform better than our method in practice.
In any case, it would not be possible to apply their theorem in all lengths
since the length n might only be divisible by primes p=2 or p≡14.

A second possible improvement could be to use the symbolic
form of f(θ) defined in Section <ref>
to help find the maximum of f(θ).  For example,
a consequence of Lemma <ref> and Euler's identity is that

    f(θ) = N_A'(0) + 2∑_s=1^n-1(N_A'(s))cos(sθ) + 2∑_s=1^n-1(N_A'(s))sin(sθ)

showing that f has a simple form as a sum of sines and cosines.
Finding the real roots of 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 approximation method in practice.
Note that our method is not guaranteed to find
the global maximum of f since we only apply the quadratic
interpolation optimization method to ranges that we
know contain a local maximum by examining the
initial sampling of f.  With a more careful branch-and-bound
algorithm it would be possible to guarantee that the global maximum
was found but the approximations that we found were effective enough
that we did not pursue this.

Another possible improvement could be obtained by deriving further properties
like Proposition <ref> that complex Golay pairs must satisfy.
For example, consider the following property
that could be viewed as a strengthening of Proposition <ref>:

    a_ka_n-k-1 = (-1)^n+1b_kb_n-k-1    for    k=1, …, n-2.

An examination of all complex Golay pairs up to length 28 reveals that
they all satisfy this property except for a single complex Golay pair
up to equivalence. 
The only pairs that don't satisfy this property are equivalent to 

    [][1, 1, 1, -1, 1, 1, -1, 1], [1, i, i, -1, 1, -i, -i, -1]

and were already singled out by <cit.> for being
the only known examples of what they call
“cross-over” Golay sequence pairs.
Since a counterexample exists to this property there is no hope of proving
it in general, but perhaps a suitable generalization could be proven.

Lastly, the running time analysis of our algorithm could be improved.
We have given a crude upper bound of 4^3n/2=8^n for the complexity
of the parts of our algorithm that do not run in polynomial time.  However,
this running time would be improved by tighter bounds on either
the size of  or the number of possible assignments that the SAT solver will try.



§ ACKNOWLEDGEMENTS


This work was made possible by the facilities of the Shared Hierarchical 
Academic Research Computing Network (SHARCNET) and Compute/Calcul Canada.
The authors would also like to thank the anonymous reviewers
for their comments on this paper and on our ISSAC 2018 submission.
Their careful reviews improved our work's clarity and presentation.



