Skip to main content
replaced http://codegolf.stackexchange.com/ with https://codegolf.stackexchange.com/
Source Link
  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus PrimeSuboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTMKennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnorxnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashmap implementation. To test 36/19 requires storing 3^18 sums, and 3^18 longs is almost 3GB without any overhead - I gave it 6GB of heap because 4GB wasn't enough; to go any further (i.e. test 38/20) within 8GB of RAM would require further optimisation to store ints rather than longs. With 20 bits required to say which subset produces the sum that would leave 12 bits plus the implicit bits from the bucket; I fear that there would be too many false collisions to get any hits.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashmap implementation. To test 36/19 requires storing 3^18 sums, and 3^18 longs is almost 3GB without any overhead - I gave it 6GB of heap because 4GB wasn't enough; to go any further (i.e. test 38/20) within 8GB of RAM would require further optimisation to store ints rather than longs. With 20 bits required to say which subset produces the sum that would leave 12 bits plus the implicit bits from the bucket; I fear that there would be too many false collisions to get any hits.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashmap implementation. To test 36/19 requires storing 3^18 sums, and 3^18 longs is almost 3GB without any overhead - I gave it 6GB of heap because 4GB wasn't enough; to go any further (i.e. test 38/20) within 8GB of RAM would require further optimisation to store ints rather than longs. With 20 bits required to say which subset produces the sum that would leave 12 bits plus the implicit bits from the bucket; I fear that there would be too many false collisions to get any hits.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
Bounty Ended with 100 reputation awarded by CommunityBot
added 314 characters in body
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168

16/9 20/11 22/12 28/15 30/16 32/17 34/18 3436/1819 (Java)

  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashsethashmap implementation specialised for. To test 36/19 requires storing 3^18 sums, and 3^18 longs is almost 3GB without any overhead longs- I gave it 6GB of heap because 4GB wasn't enough; to go any further (i.e. test 38/20) within 8GB of RAM would require further optimisation to store ints rather than longs. With 20 bits required to say which subset produces the sum that would leave 12 bits plus the implicit bits from the bucket; I fear that there would be too many false collisions to get any hits.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
import java.util.*;

// Aiming to find a solution for (2n, n+1).
public class PPCG41021_QRTernary_FixedDensity {
    private static final int N = 34;36;
    private static int density;
    private static long start;
    private static long nextProgressReport;

    public static void main(String[] args) {
        start = System.nanoTime();
        nextProgressReport = start + 5 * 60 * 1000000000L;

        // 0, -1, 1, -2, 2, ...
        for (int i = 0; i < N - 1; i++) {
            int off = i >> 1;
            if ((i & 1) == 1) off = ~off;
            density = (N >> 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j < N; j++) a[j] = j < N - density ? '0' : '1';
            c = 1;
            Bs[1] = N - density;
            Bt[1] = density;
            gen(N - density, density, 1);
            System.out.println("----");
        }

        System.out.println("Finished in " + (System.nanoTime() - start)/1000000 + " ms");
    }

    // bubble4.pdf
    private static int c;
    // Each block B_i is composed of two integers (s_i, t_i) representing the number of 0s and 1s respectively. B_1 is the rightmost block.
    // c is the number of blocks.
    // A binary string \alpha = a_1 a_2 \cdots a_m ...
    private static int[] Bs = new int[N + 1], Bt = new int[N + 1];
    private static char[] a = new char[N];
    private static void gen(int s, int t, int r) {
        if (s > 0 && t > 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i >= j; i--) {
                updateBlock(s, t, i);
                char tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                gen(s-1, t-i, testSuffix(r) ? c-1 : r);
                tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                restoreBlock(s, t, i);
            }
        }
        visit();
    }

    private static int oracle(int s, int t, int r) {
        int j = pseudoOracle(s, t, r);
        updateBlock(s, t, j);
        int p = testNecklace(testSuffix(r) ? c - 1 : r);
        restoreBlock(s, t, j);
        return p == N ? j : j + 1;
    }

    private static int pseudoOracle(int s, int t, int r) {
        if (s == 1) return t;
        if (c == 1) return s == 2 ? N / 2 : 1;
        if (s - 1 > Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) <= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) >> 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) <= 0) ? 0 : 1); 
        }
        if (s == Bs[r]) return t;
        throw new UnsupportedOperationException("Hit the case not covered by the paper or its accompanying code");
    }

    private static int testNecklace(int r) {
        if (density == 0 || density == N) return 1;
        int p = 0;
        for (int i = 0; i < c; i++) {
            if (r - i <= 0) r += c;
            if (cmpBlocks(c-i, r-i) < 0) return 0;
            if (cmpBlocks(c-i, r-1) > 0) return N;
            if (r < c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 < b1) return -1;
        if (a1 > b1) return 1;
        if (a2 < b2) return -1;
        if (a2 > b2) return 1;
        return 0;
    }

    private static int cmpBlocks(int i, int j) {
        return cmpPair(Bs[i], Bt[i], Bs[j], Bt[j]);
    }

    private static boolean testSuffix(int r) {
        for (int i = 0; i < r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) < 0) return false;
            if (cmpBlocks(c-1-i, r-1) > 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 && c > 1) {
            Bs[c-1]++;
            Bs[c] = s - 1;
        }
        else {
            Bs[c] = 1;
            Bt[c] = i;
            Bs[c+1] = s-1;
            Bt[c+1] = t-i;
            c++;
        }
    }

    private static void restoreBlock(int s, int t, int i) {
        if (i == 0 && (c > 0 || (Bs[1] != 1 || Bt[1] != 0))) {
            Bs[c-1]--;
            Bs[c] = s;
        }
        else {
            Bs[c-1] = s;
            Bt[c-1] = t;
            c--;
        }
    }

    private static long[] stats = new long[N/2+1];
    private static long visited = 0;
    private static void visit() {
        String word = new String(a);

        if (precedesReversal(word) && testTernary(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
        visited++;
        if (precedesReversal(word) && testTernary(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
        if (System.nanoTime() > nextProgressReport) {
            System.out.println("Progress: visited " + visited + "; stats " + Arrays.toString(stats) + " after " + (System.nanoTime() - start)/1000000 + " ms");
             nextProgressReport += 5 * 60 * 1000000000L;
        }
    }

    private static boolean precedesReversal(String w) {
        int n = w.length();
        StringBuilder rev = new StringBuilder(w);
        rev.reverse();
        rev.append(rev, 0, n);
        for (int i = 0; i < n; i++) {
            if (rev.substring(i, i + n).compareTo(w) < 0) return false;
        }
        return true;
    }

    private static boolean testTernary(String word) {
        int n = word.length();
        String rep = word + word;

        int base = 1;
        for (char ch : word.toCharArray()) base += ch & 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p<2^57 without overflowing a long.
        // We're computingstoring 2^n3^(n/2) ~= 2^(0.8*n) sums, so while n < 2835.56 we don't get *too* bad a probability of false reject.
        // Now that(In we'refact atthe nbirthday =paradox 28,assumes thoughindependence, it's time toand thinkour aboutvalues usingaren't twoindependent, primesso andwe're thebetter Chineseoff Remainderthan Theoremthat).
        long p = (1L << 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i < basis.length; i++) basis[i] = (basis[i-1] * base) % p;

        int rows = n / 2 + 1;
        long[] colVals = new long[n];
        for (int col = 0; col < n; col++) {
            for (int row = 0; row < rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) & 1)) % p;
            }
        }

        MapInt96Int32MapInt57Int27 setmap = new MapInt96Int32(2 * (long)Math.powMapInt57Int27(3, colVals.length / 2));
        // The first long is the sum mod p.
        // The second long will eventually contain a sum mod p2 in the lower half;
        // the top bit indicates whether it's from the left or the right, and the next n/2 bits indicate
        // which elements from that are non-zero.
        // The most significant element (or least significant? TBD) is always present only in the positive form.
        // Special-case the initial insertion.
        int[] oldLens = new int[setint[map.entries.length];
        int[] oldSupercounts = new int[1 << 10];
        {
            // count = 1
            for (int colk = 0; colk < n;n/2; col++k++) {
                longint l2val = 1 << (col25 <- n/2k);
 ? 0 : 0x8000000000000000L) | (1L << (62 - (col %     if (n/2))!map.put(colVals[k], val)); { stats[1]++; return false; }
                if (!setmap.addput(colVals[col]colVals[k + n/2], l2val + (1 << 26))) { stats[1]++; return false; }
            }
        }
        final long keyMask = (1L << 37) - 1;
        for (int count = 2; count <= n/2; count++) {
            int[] lens = setmap.counts.clone();
            int[] supercounts = map.supercounts.clone();
            for (int isup = 0; isup < lens.length;1 i++<< 10; sup++) {
                int stopunaccountedFor = 2supercounts[sup] *- lens[i];oldSupercounts[sup];
                for (int jsupi = 20; *supi < 1 << 10 && unaccountedFor > 0; supi++) {
                    int i = (sup << 10) + supi;
                    int stop = lens[i];
                    unaccountedFor -= stop - oldLens[i];
                    for (int j <= stop;oldLens[i]; j +=< 2stop; j++) {
                        long existingPexistingKV = setmap.entries[i][j];
                        long existingQexistingKey = set.entries[i][j((existingKV & keyMask) << 20) + 1];i;
                        int existingVal = (int)(existingKV >>> 37);
    
                        // For each possible prepend...
                        int half = existingQ(existingVal <>> 026) ?* n/22;
 : 0;                      // We have 27 bits of key, of which the top marks the half, so 26 bits. That means there are 6 bits at the top which we need to not count.
                        int k = LongInteger.numberOfLeadingZeros(existingQexistingVal << 16) - 1;
                        while (k >= 0) {
                        long l2   int newVal = existingQexistingVal | (1L1 << (6225 - k));
                            long pos = (existingPexistingKey + colVals[k + half]) % p;
                            if (pos << 1 > p) pos = p - pos;
                            if (pos == 0 || !setmap.addput(pos, l2newVal)) { stats[count]++; return false; }
                            long neg = (p - existingPexistingKey + colVals[k + half]) % p;
                            if (neg << 1 > p) neg = p - neg;
                            if (neg == 0 || !setmap.addput(neg, l2newVal)) { stats[count]++; return false; }
                            k--;
                        }
                    }
                }
            }
            oldLens = lens;
            oldSupercounts = supercounts;
        }

        stats[n/2]++;
        return true;
    }

    static class MapInt96Int32MapInt57Int27 {
        private long[][] entries;
        private int[] counts;

        public MapInt96Int32(long initSize) {
            int numChains = appropriateNumBuckets(initSize);
            entries = new long[numChains][];
            counts = new int[numChains];
      private int[] }supercounts;

        private static intpublic appropriateNumBucketsMapInt57Int27(long size) {
            int approxentries = 1new long[1 << (int)(0.6720][];
 * (Math.log(size) / Math.log(2)));
        counts = new int[1 return<< approx20];
 < 8 ? 8 : approx;      supercounts = new int[1 << 10];
        }

        public boolean addput(long l1key, longint l2val) {
            int bucket = bucket(l1int)(key ^& l2(entries.length - 1));
            long insert = (key >>> 20) | (((long)val) << 37);
            final long mask = (1L << 37) - 1;

            long[] chain = entries[bucket];
            if (chain == null) {
                chain = new long[16];
                entries[bucket] = chain;
                chain[0] = l1;
                chain[1] = l2;insert;
                counts[bucket]++;
                supercounts[bucket >> 10]++;
                return true;
            }

            int stop = 2 * counts[bucket];
            for (int i = 0; i < stop; i += 2i++) {
                if (chain[i] == l1 && (chain[i+1]chain[i] & 0xffffffffLmask) == (l2insert & 0xffffffffLmask)) {
                    return false;
                }
            }

            if (stop == chain.length) {
                long[] newChain = new long[chain.length < 512 ? chain.length << 1 : chain.length * 3 /+ 2];512];
                System.arraycopy(chain, 0, newChain, 0, chain.length);
                entries[bucket] = newChain;
                chain = newChain;
            }
            chain[stop] = l1;
            chain[stop+1] = l2;insert;
            counts[bucket]++;
            supercounts[bucket >> 10]++;
            return true;
        }

        protected int bucket(long obj) {
            return (int)(obj & (entries.length - 1));
        }
    }
}

The first onesone found areis

0000101000100101000110010111111111 after 12188760 ms
0000101010001100010010010111111111 after 13179066 ms
0000101010010011000010010111111111 after 13578748 ms000001001010110001000101001111111111

and that's the only hit in 15 hours.

4/3:    0111                       (plus 8 different 8/6)
9/6:    001001011                  (and 5 others)
11/7:   00010100111                (and 3 others)
13/8:   0001001101011              (and 5 others)
15/9:   000010110110111            (and 21 others)
16/9:   0000101110111011           (and 1 other)
20/11:  00000101111011110111       (and others)
22/12:  0000001100110011101011     (and others)
24/13:  000000101011101011101011   (and others)
26/14:  00000001101110010011010111 (and others)
28/15:  0000000010000111100111010111 (and others)
30/16:  000000001011001110011010101111 (and probably others)
32/17:  00001100010010100100101011111111 (and others)
34/18:  0000101000100101000110010111111111 (and others)

16/9 20/11 22/12 28/15 30/16 32/17 34/18 (Java)

  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashset implementation specialised for longs.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
import java.util.*;

// Aiming to find a solution for (2n, n+1).
public class PPCG41021_QRTernary_FixedDensity {
    private static final int N = 34;
    private static int density;
    private static long start;
    private static long nextProgressReport;

    public static void main(String[] args) {
        start = System.nanoTime();
        nextProgressReport = start + 5 * 60 * 1000000000L;

        // 0, -1, 1, -2, 2, ...
        for (int i = 0; i < N - 1; i++) {
            int off = i >> 1;
            if ((i & 1) == 1) off = ~off;
            density = (N >> 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j < N; j++) a[j] = j < N - density ? '0' : '1';
            c = 1;
            Bs[1] = N - density;
            Bt[1] = density;
            gen(N - density, density, 1);
            System.out.println("----");
        }

        System.out.println("Finished in " + (System.nanoTime() - start)/1000000 + " ms");
    }

    // bubble4.pdf
    private static int c;
    // Each block B_i is composed of two integers (s_i, t_i) representing the number of 0s and 1s respectively. B_1 is the rightmost block.
    // c is the number of blocks.
    // A binary string \alpha = a_1 a_2 \cdots a_m ...
    private static int[] Bs = new int[N + 1], Bt = new int[N + 1];
    private static char[] a = new char[N];
    private static void gen(int s, int t, int r) {
        if (s > 0 && t > 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i >= j; i--) {
                updateBlock(s, t, i);
                char tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                gen(s-1, t-i, testSuffix(r) ? c-1 : r);
                tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                restoreBlock(s, t, i);
            }
        }
        visit();
    }

    private static int oracle(int s, int t, int r) {
        int j = pseudoOracle(s, t, r);
        updateBlock(s, t, j);
        int p = testNecklace(testSuffix(r) ? c - 1 : r);
        restoreBlock(s, t, j);
        return p == N ? j : j + 1;
    }

    private static int pseudoOracle(int s, int t, int r) {
        if (s == 1) return t;
        if (c == 1) return s == 2 ? N / 2 : 1;
        if (s - 1 > Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) <= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) >> 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) <= 0) ? 0 : 1); 
        }
        if (s == Bs[r]) return t;
        throw new UnsupportedOperationException("Hit the case not covered by the paper or its accompanying code");
    }

    private static int testNecklace(int r) {
        if (density == 0 || density == N) return 1;
        int p = 0;
        for (int i = 0; i < c; i++) {
            if (r - i <= 0) r += c;
            if (cmpBlocks(c-i, r-i) < 0) return 0;
            if (cmpBlocks(c-i, r-1) > 0) return N;
            if (r < c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 < b1) return -1;
        if (a1 > b1) return 1;
        if (a2 < b2) return -1;
        if (a2 > b2) return 1;
        return 0;
    }

    private static int cmpBlocks(int i, int j) {
        return cmpPair(Bs[i], Bt[i], Bs[j], Bt[j]);
    }

    private static boolean testSuffix(int r) {
        for (int i = 0; i < r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) < 0) return false;
            if (cmpBlocks(c-1-i, r-1) > 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 && c > 1) {
            Bs[c-1]++;
            Bs[c] = s - 1;
        }
        else {
            Bs[c] = 1;
            Bt[c] = i;
            Bs[c+1] = s-1;
            Bt[c+1] = t-i;
            c++;
        }
    }

    private static void restoreBlock(int s, int t, int i) {
        if (i == 0 && (c > 0 || (Bs[1] != 1 || Bt[1] != 0))) {
            Bs[c-1]--;
            Bs[c] = s;
        }
        else {
            Bs[c-1] = s;
            Bt[c-1] = t;
            c--;
        }
    }

    private static long[] stats = new long[N/2+1];
    private static long visited = 0;
    private static void visit() {
        String word = new String(a);

        if (precedesReversal(word) && testTernary(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
        visited++;
        if (System.nanoTime() > nextProgressReport) {
            System.out.println("Progress: visited " + visited + "; stats " + Arrays.toString(stats) + " after " + (System.nanoTime() - start)/1000000 + " ms");
             nextProgressReport += 5 * 60 * 1000000000L;
        }
    }

    private static boolean precedesReversal(String w) {
        int n = w.length();
        StringBuilder rev = new StringBuilder(w);
        rev.reverse();
        rev.append(rev, 0, n);
        for (int i = 0; i < n; i++) {
            if (rev.substring(i, i + n).compareTo(w) < 0) return false;
        }
        return true;
    }

    private static boolean testTernary(String word) {
        int n = word.length();
        String rep = word + word;

        int base = 1;
        for (char ch : word.toCharArray()) base += ch & 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p<2^57 without overflowing a long.
        // We're computing 2^n sums, so while n < 28.5 we don't get *too* bad a probability of false reject.
        // Now that we're at n = 28, though, it's time to think about using two primes and the Chinese Remainder Theorem.
        long p = (1L << 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i < basis.length; i++) basis[i] = (basis[i-1] * base) % p;

        int rows = n / 2 + 1;
        long[] colVals = new long[n];
        for (int col = 0; col < n; col++) {
            for (int row = 0; row < rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) & 1)) % p;
            }
        }

        MapInt96Int32 set = new MapInt96Int32(2 * (long)Math.pow(3, colVals.length / 2));
        // The first long is the sum mod p.
        // The second long will eventually contain a sum mod p2 in the lower half;
        // the top bit indicates whether it's from the left or the right, and the next n/2 bits indicate
        // which elements from that are non-zero.
        // The most significant element (or least significant? TBD) is always present only in the positive form.
        // Special-case the initial insertion.
        int[] oldLens = new int[set.entries.length];
        {
            // count = 1
            for (int col = 0; col < n; col++) {
                long l2 = (col < n/2 ? 0 : 0x8000000000000000L) | (1L << (62 - (col % (n/2))));
                if (!set.add(colVals[col], l2)) { stats[1]++; return false; }
            }
        }
        for (int count = 2; count <= n/2; count++) {
            int[] lens = set.counts.clone();
            for (int i = 0; i < lens.length; i++) {
                int stop = 2 * lens[i];
                for (int j = 2 * oldLens[i]; j < stop; j += 2) {
                    long existingP = set.entries[i][j];
                    long existingQ = set.entries[i][j + 1];

                    // For each possible prepend...
                    int half = existingQ < 0 ? n/2 : 0;
                    int k = Long.numberOfLeadingZeros(existingQ << 1) - 1;
                    while (k >= 0) {
                        long l2 = existingQ | (1L << (62 - k));
                        long pos = (existingP + colVals[k + half]) % p;
                        if (pos << 1 > p) pos = p - pos;
                        if (!set.add(pos, l2)) { stats[count]++; return false; }
                        long neg = (p - existingP + colVals[k + half]) % p;
                        if (neg << 1 > p) neg = p - neg;
                        if (!set.add(neg, l2)) { stats[count]++; return false; }
                        k--;
                    }
                }
            }
            oldLens = lens;
        }

        stats[n/2]++;
        return true;
    }

    static class MapInt96Int32 {
        private long[][] entries;
        private int[] counts;

        public MapInt96Int32(long initSize) {
            int numChains = appropriateNumBuckets(initSize);
            entries = new long[numChains][];
            counts = new int[numChains];
        }

        private static int appropriateNumBuckets(long size) {
            int approx = 1 << (int)(0.67 * (Math.log(size) / Math.log(2)));
            return approx < 8 ? 8 : approx;
        }

        public boolean add(long l1, long l2) {
            int bucket = bucket(l1 ^ l2);

            long[] chain = entries[bucket];
            if (chain == null) {
                chain = new long[16];
                entries[bucket] = chain;
                chain[0] = l1;
                chain[1] = l2;
                counts[bucket]++;
                return true;
            }

            int stop = 2 * counts[bucket];
            for (int i = 0; i < stop; i += 2) {
                if (chain[i] == l1 && (chain[i+1] & 0xffffffffL) == (l2 & 0xffffffffL)) return false;
            }

            if (stop == chain.length) {
                long[] newChain = new long[chain.length < 512 ? chain.length << 1 : chain.length * 3 / 2];
                System.arraycopy(chain, 0, newChain, 0, chain.length);
                entries[bucket] = newChain;
                chain = newChain;
            }
            chain[stop] = l1;
            chain[stop+1] = l2;
            counts[bucket]++;
            return true;
        }

        protected int bucket(long obj) {
            return (int)(obj & (entries.length - 1));
        }
    }
}

The first ones found are

0000101000100101000110010111111111 after 12188760 ms
0000101010001100010010010111111111 after 13179066 ms
0000101010010011000010010111111111 after 13578748 ms
4/3:    0111                       (plus 8 different 8/6)
9/6:    001001011                  (and 5 others)
11/7:   00010100111                (and 3 others)
13/8:   0001001101011              (and 5 others)
15/9:   000010110110111            (and 21 others)
16/9:   0000101110111011           (and 1 other)
20/11:  00000101111011110111       (and others)
22/12:  0000001100110011101011     (and others)
24/13:  000000101011101011101011   (and others)
26/14:  00000001101110010011010111 (and others)
28/15:  0000000010000111100111010111 (and others)
30/16:  000000001011001110011010101111 (and probably others)
32/17:  00001100010010100100101011111111 (and others)

16/9 20/11 22/12 28/15 30/16 32/17 34/18 36/19 (Java)

  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashmap implementation. To test 36/19 requires storing 3^18 sums, and 3^18 longs is almost 3GB without any overhead - I gave it 6GB of heap because 4GB wasn't enough; to go any further (i.e. test 38/20) within 8GB of RAM would require further optimisation to store ints rather than longs. With 20 bits required to say which subset produces the sum that would leave 12 bits plus the implicit bits from the bucket; I fear that there would be too many false collisions to get any hits.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
import java.util.*;

// Aiming to find a solution for (2n, n+1).
public class PPCG41021_QRTernary_FixedDensity {
    private static final int N = 36;
    private static int density;
    private static long start;
    private static long nextProgressReport;

    public static void main(String[] args) {
        start = System.nanoTime();
        nextProgressReport = start + 5 * 60 * 1000000000L;

        // 0, -1, 1, -2, 2, ...
        for (int i = 0; i < N - 1; i++) {
            int off = i >> 1;
            if ((i & 1) == 1) off = ~off;
            density = (N >> 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j < N; j++) a[j] = j < N - density ? '0' : '1';
            c = 1;
            Bs[1] = N - density;
            Bt[1] = density;
            gen(N - density, density, 1);
            System.out.println("----");
        }

        System.out.println("Finished in " + (System.nanoTime() - start)/1000000 + " ms");
    }

    private static int c;
    private static int[] Bs = new int[N + 1], Bt = new int[N + 1];
    private static char[] a = new char[N];
    private static void gen(int s, int t, int r) {
        if (s > 0 && t > 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i >= j; i--) {
                updateBlock(s, t, i);
                char tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                gen(s-1, t-i, testSuffix(r) ? c-1 : r);
                tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                restoreBlock(s, t, i);
            }
        }
        visit();
    }

    private static int oracle(int s, int t, int r) {
        int j = pseudoOracle(s, t, r);
        updateBlock(s, t, j);
        int p = testNecklace(testSuffix(r) ? c - 1 : r);
        restoreBlock(s, t, j);
        return p == N ? j : j + 1;
    }

    private static int pseudoOracle(int s, int t, int r) {
        if (s == 1) return t;
        if (c == 1) return s == 2 ? N / 2 : 1;
        if (s - 1 > Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) <= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) >> 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) <= 0) ? 0 : 1); 
        }
        if (s == Bs[r]) return t;
        throw new UnsupportedOperationException("Hit the case not covered by the paper or its accompanying code");
    }

    private static int testNecklace(int r) {
        if (density == 0 || density == N) return 1;
        int p = 0;
        for (int i = 0; i < c; i++) {
            if (r - i <= 0) r += c;
            if (cmpBlocks(c-i, r-i) < 0) return 0;
            if (cmpBlocks(c-i, r-1) > 0) return N;
            if (r < c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 < b1) return -1;
        if (a1 > b1) return 1;
        if (a2 < b2) return -1;
        if (a2 > b2) return 1;
        return 0;
    }

    private static int cmpBlocks(int i, int j) {
        return cmpPair(Bs[i], Bt[i], Bs[j], Bt[j]);
    }

    private static boolean testSuffix(int r) {
        for (int i = 0; i < r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) < 0) return false;
            if (cmpBlocks(c-1-i, r-1) > 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 && c > 1) {
            Bs[c-1]++;
            Bs[c] = s - 1;
        }
        else {
            Bs[c] = 1;
            Bt[c] = i;
            Bs[c+1] = s-1;
            Bt[c+1] = t-i;
            c++;
        }
    }

    private static void restoreBlock(int s, int t, int i) {
        if (i == 0 && (c > 0 || (Bs[1] != 1 || Bt[1] != 0))) {
            Bs[c-1]--;
            Bs[c] = s;
        }
        else {
            Bs[c-1] = s;
            Bt[c-1] = t;
            c--;
        }
    }

    private static long[] stats = new long[N/2+1];
    private static long visited = 0;
    private static void visit() {
        String word = new String(a);

        visited++;
        if (precedesReversal(word) && testTernary(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
        if (System.nanoTime() > nextProgressReport) {
            System.out.println("Progress: visited " + visited + "; stats " + Arrays.toString(stats) + " after " + (System.nanoTime() - start)/1000000 + " ms");
             nextProgressReport += 5 * 60 * 1000000000L;
        }
    }

    private static boolean precedesReversal(String w) {
        int n = w.length();
        StringBuilder rev = new StringBuilder(w);
        rev.reverse();
        rev.append(rev, 0, n);
        for (int i = 0; i < n; i++) {
            if (rev.substring(i, i + n).compareTo(w) < 0) return false;
        }
        return true;
    }

    private static boolean testTernary(String word) {
        int n = word.length();
        String rep = word + word;

        int base = 1;
        for (char ch : word.toCharArray()) base += ch & 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p<2^57 without overflowing a long.
        // We're storing 3^(n/2) ~= 2^(0.8*n) sums, so while n < 35.6 we don't get *too* bad a probability of false reject.
        // (In fact the birthday paradox assumes independence, and our values aren't independent, so we're better off than that).
        long p = (1L << 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i < basis.length; i++) basis[i] = (basis[i-1] * base) % p;

        int rows = n / 2 + 1;
        long[] colVals = new long[n];
        for (int col = 0; col < n; col++) {
            for (int row = 0; row < rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) & 1)) % p;
            }
        }

        MapInt57Int27 map = new MapInt57Int27();
        // Special-case the initial insertion.
        int[] oldLens = new int[map.entries.length];
        int[] oldSupercounts = new int[1 << 10];
        {
            // count = 1
            for (int k = 0; k < n/2; k++) {
                int val = 1 << (25 - k);
                if (!map.put(colVals[k], val)) { stats[1]++; return false; }
                if (!map.put(colVals[k + n/2], val + (1 << 26))) { stats[1]++; return false; }
            }
        }
        final long keyMask = (1L << 37) - 1;
        for (int count = 2; count <= n/2; count++) {
            int[] lens = map.counts.clone();
            int[] supercounts = map.supercounts.clone();
            for (int sup = 0; sup < 1 << 10; sup++) {
                int unaccountedFor = supercounts[sup] - oldSupercounts[sup];
                for (int supi = 0; supi < 1 << 10 && unaccountedFor > 0; supi++) {
                    int i = (sup << 10) + supi;
                    int stop = lens[i];
                    unaccountedFor -= stop - oldLens[i];
                    for (int j = oldLens[i]; j < stop; j++) {
                        long existingKV = map.entries[i][j];
                        long existingKey = ((existingKV & keyMask) << 20) + i;
                        int existingVal = (int)(existingKV >>> 37);
    
                        // For each possible prepend...
                        int half = (existingVal >> 26) * n/2;
                        // We have 27 bits of key, of which the top marks the half, so 26 bits. That means there are 6 bits at the top which we need to not count.
                        int k = Integer.numberOfLeadingZeros(existingVal << 6) - 1;
                        while (k >= 0) {
                            int newVal = existingVal | (1 << (25 - k));
                            long pos = (existingKey + colVals[k + half]) % p;
                            if (pos << 1 > p) pos = p - pos;
                            if (pos == 0 || !map.put(pos, newVal)) { stats[count]++; return false; }
                            long neg = (p - existingKey + colVals[k + half]) % p;
                            if (neg << 1 > p) neg = p - neg;
                            if (neg == 0 || !map.put(neg, newVal)) { stats[count]++; return false; }
                            k--;
                        }
                    }
                }
            }
            oldLens = lens;
            oldSupercounts = supercounts;
        }

        stats[n/2]++;
        return true;
    }

    static class MapInt57Int27 {
        private long[][] entries;
        private int[] counts;
        private int[] supercounts;

        public MapInt57Int27() {
            entries = new long[1 << 20][];
            counts = new int[1 << 20];
            supercounts = new int[1 << 10];
        }

        public boolean put(long key, int val) {
            int bucket = (int)(key & (entries.length - 1));
            long insert = (key >>> 20) | (((long)val) << 37);
            final long mask = (1L << 37) - 1;

            long[] chain = entries[bucket];
            if (chain == null) {
                chain = new long[16];
                entries[bucket] = chain;
                chain[0] = insert;
                counts[bucket]++;
                supercounts[bucket >> 10]++;
                return true;
            }

            int stop = counts[bucket];
            for (int i = 0; i < stop; i++) {
                if ((chain[i] & mask) == (insert & mask)) {
                    return false;
                }
            }

            if (stop == chain.length) {
                long[] newChain = new long[chain.length < 512 ? chain.length << 1 : chain.length + 512];
                System.arraycopy(chain, 0, newChain, 0, chain.length);
                entries[bucket] = newChain;
                chain = newChain;
            }
            chain[stop] = insert;
            counts[bucket]++;
            supercounts[bucket >> 10]++;
            return true;
        }
    }
}

The first one found is

000001001010110001000101001111111111

and that's the only hit in 15 hours.

4/3:    0111                       (plus 8 different 8/6)
9/6:    001001011                  (and 5 others)
11/7:   00010100111                (and 3 others)
13/8:   0001001101011              (and 5 others)
15/9:   000010110110111            (and 21 others)
16/9:   0000101110111011           (and 1 other)
20/11:  00000101111011110111       (and others)
22/12:  0000001100110011101011     (and others)
24/13:  000000101011101011101011   (and others)
26/14:  00000001101110010011010111 (and others)
28/15:  0000000010000111100111010111 (and others)
30/16:  000000001011001110011010101111 (and probably others)
32/17:  00001100010010100100101011111111 (and others)
34/18:  0000101000100101000110010111111111 (and others)
34/18
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168

16/9 20/11 22/12 28/15 30/16 32/17 3234/1718 (Java)

This uses a number of ideas to reduce the search space and cost:. View the revision history for more details on earlier versions of the code.

  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. However, I use a k-of-n subset combinatorial Gray code (the "revolving door" technique) to iterate through them.
  • For small n I save a factor of two on the memory by not storing the sums of subsets of more than half of the columns, since they can be calculated from the sums of subsets of less than half of the columns.
  • However, this isn't enough of a memory saving for n=28 to run reliably in 3GB of memory, so whenI've now merged that idea with the k-of-n doesn't quick-reject I switch into a ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashset implementation specialised for longs.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
import java.util.*;

// Aiming to find a solution for (2n, n+1).
public class PPCG41021_FixedDensityPPCG41021_QRTernary_FixedDensity {
    private static final int N = 32;34;
    private static int density;
    private static long start;
    private static long nextProgressReport;

    public static void main(String[] args) {
        start = System.nanoTime();
        nextProgressReport = start + 5 * 60 * 1000000000L;

        // 0, -1, 1, -2, 2, ...
        for (int i = 0; i < N - 1; i++) {
            int off = i >> 1;
            if ((i & 1) == 1) off = ~off;
            density = (N >> 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j < N; j++) a[j] = j < N - density ? '0' : '1';
            c = 1;
            Bs[1] = N - density;
            Bt[1] = density;
            gen(N - density, density, 1);
            System.out.println("----");
        }

        System.out.println("Finished in " + (System.nanoTime() - start)/1000000 + " ms");
    }

    // bubble4.pdf
    private static int c;
    // Each block B_i is composed of two integers (s_i, t_i) representing the number of 0s and 1s respectively. B_1 is the rightmost block.
    // c is the number of blocks.
    // A binary string \alpha = a_1 a_2 \cdots a_m ...
    private static int[] Bs = new int[N + 1], Bt = new int[N + 1];
    private static char[] a = new char[N];
    private static void gen(int s, int t, int r) {
        if (s > 0 && t > 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i >= j; i--) {
                updateBlock(s, t, i);
                char tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                gen(s-1, t-i, testSuffix(r) ? c-1 : r);
                tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                restoreBlock(s, t, i);
            }
        }
        visit();
    }

    private static int oracle(int s, int t, int r) {
        int j = pseudoOracle(s, t, r);
        updateBlock(s, t, j);
        int p = testNecklace(testSuffix(r) ? c - 1 : r);
        restoreBlock(s, t, j);
        return p == N ? j : j + 1;
    }

    private static int pseudoOracle(int s, int t, int r) {
        if (s == 1) return t;
        if (c == 1) return s == 2 ? N / 2 : 1;
        if (s - 1 > Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) <= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) >> 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) <= 0) ? 0 : 1); 
        }
        if (s == Bs[r]) return t;
        throw new UnsupportedOperationException("Hit the case not covered by the paper or its accompanying code");
    }

    private static int testNecklace(int r) {
        if (density == 0 || density == N) return 1;
        int p = 0;
        for (int i = 0; i < c; i++) {
            if (r - i <= 0) r += c;
            if (cmpBlocks(c-i, r-i) < 0) return 0;
            if (cmpBlocks(c-i, r-1) > 0) return N;
            if (r < c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 < b1) return -1;
        if (a1 > b1) return 1;
        if (a2 < b2) return -1;
        if (a2 > b2) return 1;
        return 0;
    }

    private static int cmpBlocks(int i, int j) {
        return cmpPair(Bs[i], Bt[i], Bs[j], Bt[j]);
    }

    private static boolean testSuffix(int r) {
        for (int i = 0; i < r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) < 0) return false;
            if (cmpBlocks(c-1-i, r-1) > 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 && c > 1) {
            Bs[c-1]++;
            Bs[c] = s - 1;
        }
        else {
            Bs[c] = 1;
            Bt[c] = i;
            Bs[c+1] = s-1;
            Bt[c+1] = t-i;
            c++;
        }
    }

    private static void restoreBlock(int s, int t, int i) {
        if (i == 0 && (c > 0 || (Bs[1] != 1 || Bt[1] != 0))) {
            Bs[c-1]--;
            Bs[c] = s;
        }
        else {
            Bs[c-1] = s;
            Bt[c-1] = t;
            c--;
        }
    }

    private static long[] stats = new long[N/2+1];
    private static long visited = 0;
    private static void visit() {
        String word = new String(a);
        if (precedesReversal(word) && testHybrid(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
    }

    private static void revolvingDoorDifference(long s, int k, int n, int[] out) {
        int gap = Long.numberOfTrailingZeros(~s);
        int next = 1 + Long.numberOfTrailingZeros(s & ~((1L << gap) - 1));
        if (gap >= n) {
            out[0] = -1;
            return;
        }

        if ((precedesReversal(k - gapword) &&& 1testTernary(word) == 0) {
            if System.out.println(gap == 0) {
                out[0] = next - 1;
                out[1] = next - 2;
            }
        word + " after else" if+ (gap == 1System.nanoTime() {
                out[0] = 0;
                out[1] = 1;
            }
            else {
                out[0] = gap - 2;
                out[1] = gap;
            }
     start)/1000000 + " }ms");
        else {visited++;
            if ((s & System.nanoTime(1L << next)) ==> 0nextProgressReport) {
                if System.out.println(next == n) {
                    out[0] = -1;
                }
        "Progress: visited " + visited + "; stats else" if+ Arrays.toString(gap == 0stats) {
                    out[0] = next - 1;
                    out[1] = next;
                }
                else {
                 + " after out[0]" =+ gap(System.nanoTime() - 1;
                    out[1] = next;
                }
            }
            else {
               start)/1000000 out[0]+ =" next;ms");
                out[1] = gap;
      nextProgressReport += 5 * 60 * }1000000000L;
        }
    }

    private static boolean precedesReversal(String w) {
        int n = w.length();
        StringBuilder rev = new StringBuilder(w);
        rev.reverse();
        rev.append(rev, 0, n);
        for (int i = 0; i < n; i++) {
            if (rev.substring(i, i + n).compareTo(w) < 0) return false;
        }
        return true;
    }

    private static boolean testHybridtestTernary(String word) {
        int n = word.length();
        String rep = word + word;

        int base = 1;
        for (char ch : word.toCharArray()) base += ch & 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p<2^57 without overflowing a long.
        // We're computing 2^n sums, so while n < 28.5 we don't get *too* bad a probability of false reject.
        // (InNow factthat thewe're birthdayat paradoxn assumes= independence28, andthough, ourit's valuestime aren'tto independent,think soabout we'reusing nottwo tooprimes badlyand off)the Chinese Remainder Theorem.
        long p = (1L << 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i < basis.length; i++) basis[i] = (basis[i-1] * base) % p;

        // We're guaranteed a collision unless base^rows >= 2^n
        int rows = n / 2 + 1;
        long[] colVals = new long[n];
        for (int col = 0; col < n; col++) {
            for (int row = 0; row < rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) & 1)) % p;
            }
        }

        SetLongMapInt96Int32 sumsset = new SetLongMapInt96Int32(1L2 <<* (n - 1)long);
     Math.pow(3, colVals.length / sums.add(02));
        int[] delta = new int[2];
   // The first long is booleanthe fullyCheckedsum =mod true;p.
        int count;
  // The second long will eventually forcontain (counta =sum 1;mod 2p2 *in countthe <=lower n;half;
 count++) {
      // the top bit indicates whether ifit's (countfrom >the 10)left {
or the right, and the next n/2 bits indicate
        // Drop into the lesswhich shortcuttableelements butfrom alsothat lessare memorynon-intensive ternary mitmzero.
                fullyChecked = false;
// The most significant element (or least significant? TBD) is always present only in the positive break;form.
        // Special-case the initial }insertion.
        int[] oldLens = new int[set.entries.length];
 long mask = (1L << count) - 1;{
            long// sumcount = 0;1
            for (int icol = 0; icol < count; i++) sum = (sum + colVals[i]) % p;
            whilen; (truecol++) {
                if (!sums.add(sum)) return false;
             long l2 = revolvingDoorDifference(mask,col count,< n, delta);
                if (delta[0]/2 <? 0) break;
               : mask0x8000000000000000L) ^=| (1L << delta[0])(62 ^- (1Lcol <<% delta[1](n/2))));
                sum =if (sum + p!set.add(colVals[col], -l2)) colVals[delta[0]]{ +stats[1]++; colVals[delta[1]])return %false; p;}
            }
        }

         iffor (fullyChecked)int {
count = 2; count <= n/2; count++) {
     // The sets of more than half ofint[] thelens columns= areset.counts.clone();
 equal to the sum of all of the columns
   for (int i = 0; i < lens.length; i++) //{
 minus a set of less than half of the columns.
      int stop = 2 * lens[i];
 long all = 0;
            for (int ij = 0;2 i* oldLens[i]; j < n;stop; i++j += 2) all{
 = (all + colVals[i]) % p;
            if (sums long existingP = set.containsZeroentries[i][j];
 && sums                  long existingQ = set.contains(all))entries[i][j return+ false;1];

            for (long[] negSums      // For each possible prepend...
                    int half = existingQ < 0 ? n/2 : sums0;
                    int k = Long.entriesnumberOfLeadingZeros(existingQ << 1) {- 1;
                if    while (negSumsk ==>= null0) continue;{
                for        long l2 = existingQ | (1L << (62 - k));
                        long negSumpos := negSums(existingP + colVals[k + half]) {% p;
                        if (negSumpos ==<< 01 > p) break;pos = p - pos;
                        if (sums!set.containsadd(pos, l2)) { stats[count]++; return false; }
                        long neg = (p - existingP + allcolVals[k -+ negSumhalf]) % p;
                        if (neg << 1 > p) neg = p - neg;
                        if (!set.add(neg, l2)) { stats[count]++; return false; }
                        k--;
                    }
                }
            }
        }
        else {
            sums = null;
            SetLong setoldLens = new SetLong(2 * (long)Math.pow(3, colVals.length / 2));
            return addTernaryValues(set, colVals , 0, colVals.length / 2, p) && addTernaryValues(set, colVals, colVals.length / 2, colVals.length, p);lens;
        }

        return true;
    }

    private static boolean addTernaryValues(SetLong set, long[] colVals, int minIdx, int maxIdxExcl, long p) {
        // TODO Ternary Gray code
        int width = maxIdxExcl - minIdx;
        int count = (int)(0.5 + Math.pow(3, width));
        for (int ternary = 0; ternary < count; ternary++) {
            long sum = 0;
            int t = ternary;
            boolean nonTrivial = false;
            for (int i = 0; i < width; i++, t /= 3) {
                if (t % 3 != 1) nonTrivial = true;
                sum = (sum + p + colVals[minIdx + i] * ((t % 3) - 1)) % p;
            }
            if (!nonTrivial) continue;
            if (sum == 0) { return false; }
            // We produce each number in positive and negative forms, so a clash will occur twice.
            // Therefore we can save memory by only storing half of them.
            /stats[n/ (How can we avoid *generating* the other half?)
            if (2 * sum < p && !set.add(sum)) { return false; }
        }2]++;
        return true;
    }

    static class SetLong implements Iterable<Long>MapInt96Int32 {
        private long[][] entries;
        private boolean containsZero =int[] false;counts;

        public SetLongMapInt96Int32(long initSize) {
            int numChains = appropriateNumBuckets(initSize);
            entries = new long[numChains][];
            counts = new int[numChains];
        }

        private static int appropriateNumBuckets(long size) {
            int approx = 1 << (int)(0.67 * (Math.log(size) / Math.log(2)));
            return approx < 8 ? 8 : approx;
        }

        public boolean add(long obj) {
            if (objl1, ==long 0l2) {
                boolean rv = !containsZero;
               int containsZerobucket = true;
                return rv;
          bucket(l1 ^ }l2);

            int bucket = bucket(obj);
            long[] chain = entries[bucket];
            if (chain == null) {
                chain = new long[8];long[16];
                entries[bucket] = chain;
                chain[0] = obj;l1;
                chain[1] = l2;
                counts[bucket]++;
                return true;
            }

            booleanint spaceFoundstop = false;2 * counts[bucket];
            for (int i = 0; i < chain.length; i++) {
                long o = chain[i];
                ifstop; (oi ==+= 02) {
                   if (chain[i] = obj;
                    return true;
                }
           == l1 && (chain[i+1] & else0xffffffffL) if== (ol2 ==& obj0xffffffffL)) return false;
            }

            if (!spaceFoundstop == chain.length) {
                long[] newChain = new long[chain.length < 512 ? chain.length << 1 : chain.length * 3 / 2];
                System.arraycopy(chain, 0, newChain, 0, chain.length);
                newChain[chain.length] = obj;
                entries[bucket] = newChain;
                chain = newChain;
            }
            chain[stop] = l1;
            chain[stop+1] = l2;
            counts[bucket]++;
            return true;
        }

        protected int bucket(long obj) {
            return (int)(obj & (entries.length - 1));
        }

        public boolean contains(long obj) {
            if (obj == 0) return containsZero;

            long[] chain = entries[bucket(obj)];
            if (chain != null) {
                for (long o : chain)
                {
                    if (o == 0) break;
                    if (o == obj) return true;
                }
            }

            return false;
        }

        public Iterator<Long> iterator() {
            return new Iterator<Long>() {
                private long[] chain;
                private int nextChainIdx;
                private int idx;
                private boolean hasNext = true; // Assuming this is never called on an empty set!

                public boolean hasNext() {
                    return hasNext;
                }

                public Long next() {
                    if (chain == null) {
                        // First call.
                        while ((chain = entries[nextChainIdx++]) == null);
                        if (containsZero) return 0L;
                    }

                    long rv = chain[idx++];
                    if (idx == chain.length || chain[idx] == 0) {
                        idx = 0;
                        while (nextChainIdx < entries.length && (chain = entries[nextChainIdx++]) == null);
                        if (nextChainIdx == entries.length) {
                            hasNext = false;
                        }
                    }
                    return rv;
                }

                public void remove() {
                    throw new UnsupportedOperationException();
                }
            };
        }
    }
}
00001100010010100100101011111111 after 5641756 ms
00001010001001100011010011111111 after 9343591 ms
000001010100100110010100111111110000101000100101000110010111111111 after 990119612188760 ms
000011010010100010011000111111110000101010001100010010010111111111 after 1210972113179066 ms
000010100100111000101000111111110000101010010011000010010111111111 after 1294183313578748 ms
4/3:    0111                       (plus 8 different 8/6)
9/6:    001001011                  (and 5 others)
11/7:   00010100111                (and 3 others)
13/8:   0001001101011              (and 5 others)
15/9:   000010110110111            (and 21 others)
16/9:   0000101110111011           (and 1 other)
20/11:  00000101111011110111       (and others)
22/12:  0000001100110011101011     (and others)
24/13:  000000101011101011101011   (and others)
26/14:  00000001101110010011010111 (and others)
28/15:  0000000010000111100111010111 (and others)
30/16 :  000000001011001110011010101111 (and probably others)
32/17:  00001100010010100100101011111111 (and others)

16/9 20/11 22/12 28/15 30/16 32/17 (Java)

This uses a number of ideas to reduce the search space and cost:

  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. However, I use a k-of-n subset combinatorial Gray code (the "revolving door" technique) to iterate through them.
  • For small n I save a factor of two on the memory by not storing the sums of subsets of more than half of the columns, since they can be calculated from the sums of subsets of less than half of the columns.
  • However, this isn't enough of a memory saving for n=28 to run reliably in 3GB of memory, so when the k-of-n doesn't quick-reject I switch into a ternary subset meet-in-the-middle approach. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • I also improve memory requirements and performance by using a custom hashset implementation specialised for longs.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
import java.util.*;

// Aiming to find a solution for (2n, n+1).
public class PPCG41021_FixedDensity {
    private static final int N = 32;
    private static int density;
    private static long start;

    public static void main(String[] args) {
        start = System.nanoTime();

        // 0, -1, 1, -2, 2, ...
        for (int i = 0; i < N - 1; i++) {
            int off = i >> 1;
            if ((i & 1) == 1) off = ~off;
            density = (N >> 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j < N; j++) a[j] = j < N - density ? '0' : '1';
            c = 1;
            Bs[1] = N - density;
            Bt[1] = density;
            gen(N - density, density, 1);
        }

        System.out.println("Finished in " + (System.nanoTime() - start)/1000000 + " ms");
    }

    private static int c;
    private static int[] Bs = new int[N + 1], Bt = new int[N + 1];
    private static char[] a = new char[N];
    private static void gen(int s, int t, int r) {
        if (s > 0 && t > 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i >= j; i--) {
                updateBlock(s, t, i);
                char tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                gen(s-1, t-i, testSuffix(r) ? c-1 : r);
                tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                restoreBlock(s, t, i);
            }
        }
        visit();
    }

    private static int oracle(int s, int t, int r) {
        int j = pseudoOracle(s, t, r);
        updateBlock(s, t, j);
        int p = testNecklace(testSuffix(r) ? c - 1 : r);
        restoreBlock(s, t, j);
        return p == N ? j : j + 1;
    }

    private static int pseudoOracle(int s, int t, int r) {
        if (s == 1) return t;
        if (c == 1) return s == 2 ? N / 2 : 1;
        if (s - 1 > Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) <= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) >> 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) <= 0) ? 0 : 1); 
        }
        if (s == Bs[r]) return t;
        throw new UnsupportedOperationException("Hit the case not covered by the paper or its accompanying code");
    }

    private static int testNecklace(int r) {
        if (density == 0 || density == N) return 1;
        int p = 0;
        for (int i = 0; i < c; i++) {
            if (r - i <= 0) r += c;
            if (cmpBlocks(c-i, r-i) < 0) return 0;
            if (cmpBlocks(c-i, r-1) > 0) return N;
            if (r < c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 < b1) return -1;
        if (a1 > b1) return 1;
        if (a2 < b2) return -1;
        if (a2 > b2) return 1;
        return 0;
    }

    private static int cmpBlocks(int i, int j) {
        return cmpPair(Bs[i], Bt[i], Bs[j], Bt[j]);
    }

    private static boolean testSuffix(int r) {
        for (int i = 0; i < r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) < 0) return false;
            if (cmpBlocks(c-1-i, r-1) > 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 && c > 1) {
            Bs[c-1]++;
            Bs[c] = s - 1;
        }
        else {
            Bs[c] = 1;
            Bt[c] = i;
            Bs[c+1] = s-1;
            Bt[c+1] = t-i;
            c++;
        }
    }

    private static void restoreBlock(int s, int t, int i) {
        if (i == 0 && (c > 0 || (Bs[1] != 1 || Bt[1] != 0))) {
            Bs[c-1]--;
            Bs[c] = s;
        }
        else {
            Bs[c-1] = s;
            Bt[c-1] = t;
            c--;
        }
    }

    private static void visit() {
        String word = new String(a);
        if (precedesReversal(word) && testHybrid(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
    }

    private static void revolvingDoorDifference(long s, int k, int n, int[] out) {
        int gap = Long.numberOfTrailingZeros(~s);
        int next = 1 + Long.numberOfTrailingZeros(s & ~((1L << gap) - 1));
        if (gap >= n) {
            out[0] = -1;
            return;
        }

        if (((k - gap) & 1) == 0) {
            if (gap == 0) {
                out[0] = next - 1;
                out[1] = next - 2;
            }
            else if (gap == 1) {
                out[0] = 0;
                out[1] = 1;
            }
            else {
                out[0] = gap - 2;
                out[1] = gap;
            }
        }
        else {
            if ((s & (1L << next)) == 0) {
                if (next == n) {
                    out[0] = -1;
                }
                else if (gap == 0) {
                    out[0] = next - 1;
                    out[1] = next;
                }
                else {
                    out[0] = gap - 1;
                    out[1] = next;
                }
            }
            else {
                out[0] = next;
                out[1] = gap;
            }
        }
    }

    private static boolean precedesReversal(String w) {
        int n = w.length();
        StringBuilder rev = new StringBuilder(w);
        rev.reverse();
        rev.append(rev, 0, n);
        for (int i = 0; i < n; i++) {
            if (rev.substring(i, i + n).compareTo(w) < 0) return false;
        }
        return true;
    }

    private static boolean testHybrid(String word) {
        int n = word.length();
        String rep = word + word;

        int base = 1;
        for (char ch : word.toCharArray()) base += ch & 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p<2^57 without overflowing a long.
        // We're computing 2^n sums, so while n < 28.5 we don't get *too* bad a probability of false reject.
        // (In fact the birthday paradox assumes independence, and our values aren't independent, so we're not too badly off).
        long p = (1L << 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i < basis.length; i++) basis[i] = (basis[i-1] * base) % p;

        // We're guaranteed a collision unless base^rows >= 2^n
        int rows = n / 2 + 1;
        long[] colVals = new long[n];
        for (int col = 0; col < n; col++) {
            for (int row = 0; row < rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) & 1)) % p;
            }
        }

        SetLong sums = new SetLong(1L << (n - 1));
        sums.add(0);
        int[] delta = new int[2];
        boolean fullyChecked = true;
        int count;
        for (count = 1; 2 * count <= n; count++) {
            if (count > 10) {
                // Drop into the less shortcuttable but also less memory-intensive ternary mitm.
                fullyChecked = false;
                break;
            }
            long mask = (1L << count) - 1;
            long sum = 0;
            for (int i = 0; i < count; i++) sum = (sum + colVals[i]) % p;
            while (true) {
                if (!sums.add(sum)) return false;
                revolvingDoorDifference(mask, count, n, delta);
                if (delta[0] < 0) break;
                mask ^= (1L << delta[0]) ^ (1L << delta[1]);
                sum = (sum + p - colVals[delta[0]] + colVals[delta[1]]) % p;
            }
        }

         if (fullyChecked) {
            // The sets of more than half of the columns are equal to the sum of all of the columns
            // minus a set of less than half of the columns.
            long all = 0;
            for (int i = 0; i < n; i++) all = (all + colVals[i]) % p;
            if (sums.containsZero && sums.contains(all)) return false;

            for (long[] negSums : sums.entries) {
                if (negSums == null) continue;
                for (long negSum : negSums) {
                    if (negSum == 0) break;
                    if (sums.contains((p + all - negSum) % p)) return false;
                }
            }
        }
        else {
            sums = null;
            SetLong set = new SetLong(2 * (long)Math.pow(3, colVals.length / 2));
            return addTernaryValues(set, colVals , 0, colVals.length / 2, p) && addTernaryValues(set, colVals, colVals.length / 2, colVals.length, p);
        }

        return true;
    }

    private static boolean addTernaryValues(SetLong set, long[] colVals, int minIdx, int maxIdxExcl, long p) {
        // TODO Ternary Gray code
        int width = maxIdxExcl - minIdx;
        int count = (int)(0.5 + Math.pow(3, width));
        for (int ternary = 0; ternary < count; ternary++) {
            long sum = 0;
            int t = ternary;
            boolean nonTrivial = false;
            for (int i = 0; i < width; i++, t /= 3) {
                if (t % 3 != 1) nonTrivial = true;
                sum = (sum + p + colVals[minIdx + i] * ((t % 3) - 1)) % p;
            }
            if (!nonTrivial) continue;
            if (sum == 0) { return false; }
            // We produce each number in positive and negative forms, so a clash will occur twice.
            // Therefore we can save memory by only storing half of them.
            // (How can we avoid *generating* the other half?)
            if (2 * sum < p && !set.add(sum)) { return false; }
        }
        return true;
    }

    static class SetLong implements Iterable<Long> {
        private long[][] entries;
        private boolean containsZero = false;

        public SetLong(long initSize) {
            int numChains = appropriateNumBuckets(initSize);
            entries = new long[numChains][];
        }

        private static int appropriateNumBuckets(long size) {
            int approx = 1 << (int)(0.67 * (Math.log(size) / Math.log(2)));
            return approx < 8 ? 8 : approx;
        }

        public boolean add(long obj) {
            if (obj == 0) {
                boolean rv = !containsZero;
                containsZero = true;
                return rv;
            }

            int bucket = bucket(obj);
            long[] chain = entries[bucket];
            if (chain == null) {
                chain = new long[8];
                entries[bucket] = chain;
                chain[0] = obj;
                return true;
            }

            boolean spaceFound = false;
            for (int i = 0; i < chain.length; i++) {
                long o = chain[i];
                if (o == 0) {
                    chain[i] = obj;
                    return true;
                }
                else if (o == obj) return false;
            }

            if (!spaceFound) {
                long[] newChain = new long[chain.length < 512 ? chain.length << 1 : chain.length * 3 / 2];
                System.arraycopy(chain, 0, newChain, 0, chain.length);
                newChain[chain.length] = obj;
                entries[bucket] = newChain;
            }

            return true;
        }

        protected int bucket(long obj) {
            return (int)(obj & (entries.length - 1));
        }

        public boolean contains(long obj) {
            if (obj == 0) return containsZero;

            long[] chain = entries[bucket(obj)];
            if (chain != null) {
                for (long o : chain)
                {
                    if (o == 0) break;
                    if (o == obj) return true;
                }
            }

            return false;
        }

        public Iterator<Long> iterator() {
            return new Iterator<Long>() {
                private long[] chain;
                private int nextChainIdx;
                private int idx;
                private boolean hasNext = true; // Assuming this is never called on an empty set!

                public boolean hasNext() {
                    return hasNext;
                }

                public Long next() {
                    if (chain == null) {
                        // First call.
                        while ((chain = entries[nextChainIdx++]) == null);
                        if (containsZero) return 0L;
                    }

                    long rv = chain[idx++];
                    if (idx == chain.length || chain[idx] == 0) {
                        idx = 0;
                        while (nextChainIdx < entries.length && (chain = entries[nextChainIdx++]) == null);
                        if (nextChainIdx == entries.length) {
                            hasNext = false;
                        }
                    }
                    return rv;
                }

                public void remove() {
                    throw new UnsupportedOperationException();
                }
            };
        }
    }
}
00001100010010100100101011111111 after 5641756 ms
00001010001001100011010011111111 after 9343591 ms
00000101010010011001010011111111 after 9901196 ms
00001101001010001001100011111111 after 12109721 ms
00001010010011100010100011111111 after 12941833 ms
4/3:    0111                       (plus 8 different 8/6)
9/6:    001001011                  (and 5 others)
11/7:   00010100111                (and 3 others)
13/8:   0001001101011              (and 5 others)
15/9:   000010110110111            (and 21 others)
16/9:   0000101110111011           (and 1 other)
20/11:  00000101111011110111       (and others)
22/12:  0000001100110011101011     (and others)
24/13:  000000101011101011101011   (and others)
26/14:  00000001101110010011010111 (and others)
28/15:  0000000010000111100111010111 (and others)
30/16   000000001011001110011010101111 (and probably others)

16/9 20/11 22/12 28/15 30/16 32/17 34/18 (Java)

This uses a number of ideas to reduce the search space and cost. View the revision history for more details on earlier versions of the code.

  • It's clear that wlog we can consider only circulant matrices in which the first row is a Lyndon word: if the word is non-prime then it must have property X, and otherwise we can rotate without affecting the score or property X.
  • Based on heuristics from the observed short winners, I'm now iterating through the Lyndon words starting at the ones with 50% density (i.e. the same number of 0 and 1) and working out; I use the algorithm described in A Gray code for fixed-density necklaces and Lyndon words in constant amortized time, Sawada and Williams, Theoretical Computer Science 502 (2013): 46-54.
  • An empirical observation is that the values occur in pairs: each optimum Lyndon word that I've found scores the same as its reversal. So I get about a factor of two speedup by only considering one half of each such pair.
  • My original code worked with BigInteger to give an exact test. I get a significant speed improvement, at the risk of false negatives, by operating modulo a large prime and keeping everything in primitives. The prime I've chosen is the largest one smaller than 257, as that allows multiplying by the base of my notional vector representation without overflowing.
  • I've stolen Suboptimus Prime's heuristic that it's possible to get quick rejections by considering subsets in increasing order of size. I've now merged that idea with the ternary subset meet-in-the-middle approach to test for colliding subsets. (Credit to KennyTM for suggesting trying to adapt the approach from the integer subset problem; I think that xnor and I saw the way to do it pretty much simultaneously). Rather than looking for two subsets which can include each column 0 or 1 times and have the same sum, we look for one subset which can include each column -1, 0, or 1 times and sum to zero. This significantly reduces the memory requirements.
  • There's an extra factor of two saving in memory requirements by observing that since each element in {-1,0,1}^m has its negation also in {-1,0,1}^m it's only necessary to store one of the two.
  • I also improve memory requirements and performance by using a custom hashset implementation specialised for longs.
  • Since the weight of the evidence suggests that we should look at 2n/(n+1), I'm speeding things up by just testing that.
  • There's some unnecessary but reassuring statistical output.
import java.util.*;

// Aiming to find a solution for (2n, n+1).
public class PPCG41021_QRTernary_FixedDensity {
    private static final int N = 34;
    private static int density;
    private static long start;
    private static long nextProgressReport;

    public static void main(String[] args) {
        start = System.nanoTime();
        nextProgressReport = start + 5 * 60 * 1000000000L;

        // 0, -1, 1, -2, 2, ...
        for (int i = 0; i < N - 1; i++) {
            int off = i >> 1;
            if ((i & 1) == 1) off = ~off;
            density = (N >> 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j < N; j++) a[j] = j < N - density ? '0' : '1';
            c = 1;
            Bs[1] = N - density;
            Bt[1] = density;
            gen(N - density, density, 1);
            System.out.println("----");
        }

        System.out.println("Finished in " + (System.nanoTime() - start)/1000000 + " ms");
    }

    // bubble4.pdf
    private static int c;
    // Each block B_i is composed of two integers (s_i, t_i) representing the number of 0s and 1s respectively. B_1 is the rightmost block.
    // c is the number of blocks.
    // A binary string \alpha = a_1 a_2 \cdots a_m ...
    private static int[] Bs = new int[N + 1], Bt = new int[N + 1];
    private static char[] a = new char[N];
    private static void gen(int s, int t, int r) {
        if (s > 0 && t > 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i >= j; i--) {
                updateBlock(s, t, i);
                char tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                gen(s-1, t-i, testSuffix(r) ? c-1 : r);
                tmp = a[s - 1]; a[s - 1] = a[s+t-i - 1]; a[s+t-i - 1] = tmp;
                restoreBlock(s, t, i);
            }
        }
        visit();
    }

    private static int oracle(int s, int t, int r) {
        int j = pseudoOracle(s, t, r);
        updateBlock(s, t, j);
        int p = testNecklace(testSuffix(r) ? c - 1 : r);
        restoreBlock(s, t, j);
        return p == N ? j : j + 1;
    }

    private static int pseudoOracle(int s, int t, int r) {
        if (s == 1) return t;
        if (c == 1) return s == 2 ? N / 2 : 1;
        if (s - 1 > Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) <= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) >> 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) <= 0) ? 0 : 1); 
        }
        if (s == Bs[r]) return t;
        throw new UnsupportedOperationException("Hit the case not covered by the paper or its accompanying code");
    }

    private static int testNecklace(int r) {
        if (density == 0 || density == N) return 1;
        int p = 0;
        for (int i = 0; i < c; i++) {
            if (r - i <= 0) r += c;
            if (cmpBlocks(c-i, r-i) < 0) return 0;
            if (cmpBlocks(c-i, r-1) > 0) return N;
            if (r < c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 < b1) return -1;
        if (a1 > b1) return 1;
        if (a2 < b2) return -1;
        if (a2 > b2) return 1;
        return 0;
    }

    private static int cmpBlocks(int i, int j) {
        return cmpPair(Bs[i], Bt[i], Bs[j], Bt[j]);
    }

    private static boolean testSuffix(int r) {
        for (int i = 0; i < r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) < 0) return false;
            if (cmpBlocks(c-1-i, r-1) > 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 && c > 1) {
            Bs[c-1]++;
            Bs[c] = s - 1;
        }
        else {
            Bs[c] = 1;
            Bt[c] = i;
            Bs[c+1] = s-1;
            Bt[c+1] = t-i;
            c++;
        }
    }

    private static void restoreBlock(int s, int t, int i) {
        if (i == 0 && (c > 0 || (Bs[1] != 1 || Bt[1] != 0))) {
            Bs[c-1]--;
            Bs[c] = s;
        }
        else {
            Bs[c-1] = s;
            Bt[c-1] = t;
            c--;
        }
    }

    private static long[] stats = new long[N/2+1];
    private static long visited = 0;
    private static void visit() {
        String word = new String(a);

        if (precedesReversal(word) && testTernary(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
        visited++;
        if (System.nanoTime() > nextProgressReport) {
            System.out.println("Progress: visited " + visited + "; stats " + Arrays.toString(stats) + " after " + (System.nanoTime() - start)/1000000 + " ms");
             nextProgressReport += 5 * 60 * 1000000000L;
        }
    }

    private static boolean precedesReversal(String w) {
        int n = w.length();
        StringBuilder rev = new StringBuilder(w);
        rev.reverse();
        rev.append(rev, 0, n);
        for (int i = 0; i < n; i++) {
            if (rev.substring(i, i + n).compareTo(w) < 0) return false;
        }
        return true;
    }

    private static boolean testTernary(String word) {
        int n = word.length();
        String rep = word + word;

        int base = 1;
        for (char ch : word.toCharArray()) base += ch & 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p<2^57 without overflowing a long.
        // We're computing 2^n sums, so while n < 28.5 we don't get *too* bad a probability of false reject.
        // Now that we're at n = 28, though, it's time to think about using two primes and the Chinese Remainder Theorem.
        long p = (1L << 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i < basis.length; i++) basis[i] = (basis[i-1] * base) % p;

        int rows = n / 2 + 1;
        long[] colVals = new long[n];
        for (int col = 0; col < n; col++) {
            for (int row = 0; row < rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) & 1)) % p;
            }
        }

        MapInt96Int32 set = new MapInt96Int32(2 * (long)Math.pow(3, colVals.length / 2));
        // The first long is the sum mod p.
        // The second long will eventually contain a sum mod p2 in the lower half;
        // the top bit indicates whether it's from the left or the right, and the next n/2 bits indicate
        // which elements from that are non-zero.
        // The most significant element (or least significant? TBD) is always present only in the positive form.
        // Special-case the initial insertion.
        int[] oldLens = new int[set.entries.length];
        {
            // count = 1
            for (int col = 0; col < n; col++) {
                long l2 = (col < n/2 ? 0 : 0x8000000000000000L) | (1L << (62 - (col % (n/2))));
                if (!set.add(colVals[col], l2)) { stats[1]++; return false; }
            }
        }
        for (int count = 2; count <= n/2; count++) {
            int[] lens = set.counts.clone();
            for (int i = 0; i < lens.length; i++) {
                int stop = 2 * lens[i];
                for (int j = 2 * oldLens[i]; j < stop; j += 2) {
                    long existingP = set.entries[i][j];
                    long existingQ = set.entries[i][j + 1];

                    // For each possible prepend...
                    int half = existingQ < 0 ? n/2 : 0;
                    int k = Long.numberOfLeadingZeros(existingQ << 1) - 1;
                    while (k >= 0) {
                        long l2 = existingQ | (1L << (62 - k));
                        long pos = (existingP + colVals[k + half]) % p;
                        if (pos << 1 > p) pos = p - pos;
                        if (!set.add(pos, l2)) { stats[count]++; return false; }
                        long neg = (p - existingP + colVals[k + half]) % p;
                        if (neg << 1 > p) neg = p - neg;
                        if (!set.add(neg, l2)) { stats[count]++; return false; }
                        k--;
                    }
                }
            }
            oldLens = lens;
        }

        stats[n/2]++;
        return true;
    }

    static class MapInt96Int32 {
        private long[][] entries;
        private int[] counts;

        public MapInt96Int32(long initSize) {
            int numChains = appropriateNumBuckets(initSize);
            entries = new long[numChains][];
            counts = new int[numChains];
        }

        private static int appropriateNumBuckets(long size) {
            int approx = 1 << (int)(0.67 * (Math.log(size) / Math.log(2)));
            return approx < 8 ? 8 : approx;
        }

        public boolean add(long l1, long l2) {
            int bucket = bucket(l1 ^ l2);

            long[] chain = entries[bucket];
            if (chain == null) {
                chain = new long[16];
                entries[bucket] = chain;
                chain[0] = l1;
                chain[1] = l2;
                counts[bucket]++;
                return true;
            }

            int stop = 2 * counts[bucket];
            for (int i = 0; i < stop; i += 2) {
                if (chain[i] == l1 && (chain[i+1] & 0xffffffffL) == (l2 & 0xffffffffL)) return false;
            }

            if (stop == chain.length) {
                long[] newChain = new long[chain.length < 512 ? chain.length << 1 : chain.length * 3 / 2];
                System.arraycopy(chain, 0, newChain, 0, chain.length);
                entries[bucket] = newChain;
                chain = newChain;
            }
            chain[stop] = l1;
            chain[stop+1] = l2;
            counts[bucket]++;
            return true;
        }

        protected int bucket(long obj) {
            return (int)(obj & (entries.length - 1));
        }
    }
}
0000101000100101000110010111111111 after 12188760 ms
0000101010001100010010010111111111 after 13179066 ms
0000101010010011000010010111111111 after 13578748 ms
4/3:    0111                       (plus 8 different 8/6)
9/6:    001001011                  (and 5 others)
11/7:   00010100111                (and 3 others)
13/8:   0001001101011              (and 5 others)
15/9:   000010110110111            (and 21 others)
16/9:   0000101110111011           (and 1 other)
20/11:  00000101111011110111       (and others)
22/12:  0000001100110011101011     (and others)
24/13:  000000101011101011101011   (and others)
26/14:  00000001101110010011010111 (and others)
28/15:  0000000010000111100111010111 (and others)
30/16:  000000001011001110011010101111 (and probably others)
32/17:  00001100010010100100101011111111 (and others)
Ahem
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
added 1598 characters in body
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
added 785 characters in body
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
added 9913 characters in body
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
added 287 characters in body
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
Bugfix, performance improvement by quick-reject
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
added 319 characters in body
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading
Source Link
Peter Taylor
  • 43.1k
  • 4
  • 70
  • 168
Loading