Find highest scoring matrix without property X This challenge is partly an algorithms challenge, partly an optimization challenge and partly simply a fastest code challenge.

A cyclic matrix is fully specified by its first row r. The remaining rows are each cyclic permutations of the row r with offset equal to the row index. We will allow cyclic matrices which are not square so that they are simply missing some of their last rows. We do however always assume that the number of rows is no more than the number of columns. For example, consider the following 3 by 5 cyclic matrix.

10111
11011
11101


We say a matrix has property X if it contains two non-empty sets of columns with non-identical indices which have the same (vector) sum. The vector sum of two columns is simply an element-wise summation of the two columns.  That is the sum of two columns containing x elements each is another column containing x elements. 

The matrix above trivially has property X as the first and last columns are the same. The identity matrix never has property X. 

If we just remove the last column of the matrix above then we get an example which does not have property X and would give a score of 4/3.

1011
1101
1110


The task

The task is to write code to find the highest scoring cyclic matrix whose entries are all 0 or 1 and which does not have property X.

Score 

Your score will be the number columns divided by the number of rows in your best scoring matrix.

Tie Breaker

If two answers have the same score, the one submitted first wins.

In the (very) unlikely event that someone finds a method to get unlimited scores, the first valid proof of such a solution will be accepted.  In the even more unlikely event that you can find a proof of optimality of a finite matrix I will of course award the win too.

Hint

Getting a score of 12/8 is not too hard. 

Languages and libraries

You can use any language which has a freely available compiler/interpreter/etc. for Linux and any libraries which are also freely available for Linux.

Leading entries


36/19 by Peter Taylor (Java)
32/17 by Suboptimus Prime (C#)
21/12 by justhalf (Python 2)
16/9 20/11 22/12 28/15 30/16 32/17 34/18 36/19 (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 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 &lt; N - 1; i++) {
            int off = i &gt;&gt; 1;
            if ((i &amp; 1) == 1) off = ~off;
            density = (N &gt;&gt; 1) + off;

            // Iterate over Lyndon words of length N and given density.
            for (int j = 0; j &lt; N; j++) a[j] = j &lt; 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 &gt; 0 &amp;&amp; t &gt; 0) {
            int j = oracle(s, t, r);
            for (int i = t - 1; i &gt;= 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 &gt; Bs[r] + 1) return 0;
        if (s - 1 == Bs[r] + 1) return cmpPair(s-1, t, Bs[c-1]+1, Bt[c-1]) &lt;= 0 ? 0 : 1;
        if (s - 1 == Bs[r]) {
            if (s == 2) return Math.max(t - Bt[r], (t+1) &gt;&gt; 1);
            return Math.max(t - Bt[r], (cmpPair(s-1, t, Bs[c-1] + 1, Bt[c-1]) &lt;= 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 &lt; c; i++) {
            if (r - i &lt;= 0) r += c;
            if (cmpBlocks(c-i, r-i) &lt; 0) return 0;
            if (cmpBlocks(c-i, r-1) &gt; 0) return N;
            if (r &lt; c) p += Bs[r-i] + Bt[r-i];
        }
        return p;
    }

    private static int cmpPair(int a1, int a2, int b1, int b2) {
        if (a1 &lt; b1) return -1;
        if (a1 &gt; b1) return 1;
        if (a2 &lt; b2) return -1;
        if (a2 &gt; 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 &lt; r; i++) {
            if (c - 1 - i == r) return true;
            if (cmpBlocks(c-1-i, r-i) &lt; 0) return false;
            if (cmpBlocks(c-1-i, r-1) &gt; 0) return true;
        }
        return false;
    }

    private static void updateBlock(int s, int t, int i) {
        if (i == 0 &amp;&amp; c &gt; 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 &amp;&amp; (c &gt; 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) &amp;&amp; testTernary(word)) System.out.println(word + " after " + (System.nanoTime() - start)/1000000 + " ms");
        if (System.nanoTime() &gt; 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 &lt; n; i++) {
            if (rev.substring(i, i + n).compareTo(w) &lt; 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 &amp; 1;

        // Operating base b for b up to 32 implies that we can multiply by b modulo p&lt;2^57 without overflowing a long.
        // We're storing 3^(n/2) ~= 2^(0.8*n) sums, so while n &lt; 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 &lt;&lt; 57) - 13;
        long[] basis = new long[n];
        basis[0] = 1;
        for (int i = 1; i &lt; 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 &lt; n; col++) {
            for (int row = 0; row &lt; rows; row++) {
                colVals[col] = (colVals[col] + basis[row] * (rep.charAt(row + col) &amp; 1)) % p;
            }
        }

        MapInt57Int27 map = new MapInt57Int27();
        // Special-case the initial insertion.
        int[] oldLens = new int[map.entries.length];
        int[] oldSupercounts = new int[1 &lt;&lt; 10];
        {
            // count = 1
            for (int k = 0; k &lt; n/2; k++) {
                int val = 1 &lt;&lt; (25 - k);
                if (!map.put(colVals[k], val)) { stats[1]++; return false; }
                if (!map.put(colVals[k + n/2], val + (1 &lt;&lt; 26))) { stats[1]++; return false; }
            }
        }
        final long keyMask = (1L &lt;&lt; 37) - 1;
        for (int count = 2; count &lt;= n/2; count++) {
            int[] lens = map.counts.clone();
            int[] supercounts = map.supercounts.clone();
            for (int sup = 0; sup &lt; 1 &lt;&lt; 10; sup++) {
                int unaccountedFor = supercounts[sup] - oldSupercounts[sup];
                for (int supi = 0; supi &lt; 1 &lt;&lt; 10 &amp;&amp; unaccountedFor &gt; 0; supi++) {
                    int i = (sup &lt;&lt; 10) + supi;
                    int stop = lens[i];
                    unaccountedFor -= stop - oldLens[i];
                    for (int j = oldLens[i]; j &lt; stop; j++) {
                        long existingKV = map.entries[i][j];
                        long existingKey = ((existingKV &amp; keyMask) &lt;&lt; 20) + i;
                        int existingVal = (int)(existingKV &gt;&gt;&gt; 37);

                        // For each possible prepend...
                        int half = (existingVal &gt;&gt; 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 &lt;&lt; 6) - 1;
                        while (k &gt;= 0) {
                            int newVal = existingVal | (1 &lt;&lt; (25 - k));
                            long pos = (existingKey + colVals[k + half]) % p;
                            if (pos &lt;&lt; 1 &gt; 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 &lt;&lt; 1 &gt; 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 &lt;&lt; 20][];
            counts = new int[1 &lt;&lt; 20];
            supercounts = new int[1 &lt;&lt; 10];
        }

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

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

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

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


The first one found is

000001001010110001000101001111111111


and that's the only hit in 15 hours.

Smaller winners:

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)