Find highest scoring matrix without property X

14

2

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)

user9206

Posted 2014-11-05T07:19:32.653

Reputation:

Ah, property X is on columns, not rows. – Optimizer – 2014-11-05T07:26:46.750

As written, the 1 by 2 matrix 01 has property X because the set of the first column has the same vector sum as the empty set. Perhaps you meant nonempty sets of columns? I think it's cleaner not to change it though. – xnor – 2014-11-05T07:50:11.363

@xnor Good point. I have edited the question as you suggested. – None – 2014-11-05T09:09:22.650

Wow, do you actually know whether there is a solution with score higher than 2? – justhalf – 2014-11-05T09:16:09.640

2The easiest reading of the rules is still that 01 has property X: (1) = (0) + (1). If you want to exclude that then you should say that the two sets of columns must be disjoint. – Peter Taylor – 2014-11-05T09:19:59.067

@PeterTaylor That was just a mistake by me. I have deleted the offending line in the question. Thank you for spotting this. – None – 2014-11-05T09:21:19.943

@Lembik: Btw, why do you write 12/8 instead of 4/3? – justhalf – 2014-11-05T09:21:45.823

@justhalf 12/8 > 4/3 :) I am saying that you can find an 8 by 12 matrix without property X. – None – 2014-11-05T09:29:02.037

Oops, haha, miscalculated. Sorry for the confusion! – justhalf – 2014-11-05T09:35:19.537

pretty sure its immposible to get a higher score then 2 – Vajura – 2014-11-05T09:45:47.097

I believe we can get infinite score. – justhalf – 2014-11-05T09:48:09.360

Hm you might be right i was trying to get a general solution with linear algebra with linear independance but the rules dont apply actualy which i just noticed. I am very sure tho you can get a solution which will show if its infinite or a finite with a specific number – Vajura – 2014-11-05T09:56:19.553

Oops, mistake in assumption. My method doesn't generate infinite result, but at least it finds something better than 12/8. – justhalf – 2014-11-05T10:04:57.347

Wait... can the matrix contain numbers different than 0 and 1? – Gerwin – 2014-11-05T15:16:24.363

No it can not. . – None – 2014-11-05T15:17:28.997

I tried to prove a finite upper bound by pigeonhole, but am a log factor off. The idea is that a wide enough matrix has more column-subsets of a given size then possible sums of that many column vectors, so there must be a repeated sum. – xnor – 2014-11-05T18:51:41.977

@xnor I don't believe there is a finite upper bound, but I don't know how to prove that either. Which is why I hope the contest will prove fun :) – None – 2014-11-05T19:13:33.207

I can at least show that obtaining a score of x requires a matrix whose size is asymptotically exponential in x, so matrices achieving large scores would have to be enormous. – xnor – 2014-11-05T19:23:17.447

@xnor That is very interesting and I would love to see it. Of course at this point a score of 3 would be very impressive. Does your method say how large a matrix would be needed for that? – None – 2014-11-05T19:25:13.130

Unfortunately, it only gives trivial bounds for scores like 3. For a score of 5, it gives a minimum of 4x20, and for a score of 10, it gives 62x620. – xnor – 2014-11-05T20:22:02.307

1

This question will give much insight on this problem (on how hard it is to check property X, which is NP-hard, unfortunately) http://mathoverflow.net/questions/157634/number-of-vectors-so-that-no-two-subset-sums-are-equal

– justhalf – 2014-11-06T08:51:58.113

@justhalf This is interesting but not directly relevant. The problem they consider does not include the restriction that the matrix is cyclic which might make it a lot easier, who knows. In particular, we don't know if the problem of telling if a cyclic matrix has property X is NP-hard or not. My guess is that it is not. – None – 2014-11-06T11:43:52.520

Yes, it doesn't deal with the cyclic nature. Let's see whether anyone can find a polynomial time algorithm to check property X on cyclic matrices. – justhalf – 2014-11-06T11:52:41.890

3Currently we are just brute-forcing all the 2^m column combinations to check property X. If we could somehow devise a "meet in the middle" scheme (see the "subset sum" problem) this could probably reduce that to m * 2^(m/2)... – kennytm – 2014-11-07T20:11:39.913

Answers

11

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 < 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.

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)

Peter Taylor

Posted 2014-11-05T07:19:32.653

Reputation: 41 901

This is a good improvement. It seems using Lyndon words means you only need to check roughly 2^n/n binary strings for the first row, instead of 2^n. – None – 2014-11-05T20:15:09.847

As you are using each digit of BigInteger as a matrix cell, won't there be wrong answer when n > 10? – kennytm – 2014-11-05T20:15:22.483

@KennyTM, note that the second parameter is the radix. There is a small bug: I should be using n rather than rows, although it's fail-safe in the sense that it would discard valid solutions rather than accept invalid ones. It also doesn't affect the results. – Peter Taylor – 2014-11-05T21:50:00.087

Can you explain more about the "non-prime iff no property X" that you seem to imply in your sentences? – justhalf – 2014-11-06T00:18:53.647

@justhalf, I don't think that's implied and I certainly don't intend to imply it. A prime word can have property X or not: the point is that rotating it doesn't change whether or not it has property X, so we can consider the Lyndon word in the rotational set as a representative element. – Peter Taylor – 2014-11-06T07:18:34.437

I see, so the purpose of that sentence is only to claim that the only first rows that we need to brute force are all those Lyndon words. – justhalf – 2014-11-06T07:23:05.863

1I think we're practically limited to this score, since the property X checking is very time consuming, unless we found another equivalent condition which can be evaluated faster. That's why I was so eager to see that "non-prime" implies property X =D – justhalf – 2014-11-06T07:55:35.560

@justhalf, I reckon that I can check up to MAX_N=24 in 12 hours on a fast computer using one thread, but that seems to be about as far as I can go because of memory constraints. – Peter Taylor – 2014-11-06T08:17:57.543

It isn't at all obvious to me why "we can consider only circulant matrices in which the first row is a Lyndon word". I certainly believe you, but would you indulge me in explaining why property X implies Lyndon wordity? I'd appreciate it. :) – COTO – 2014-11-09T16:59:59.597

@COTO, that's not the implication. There are three types of words: non-prime words, Lyndon words, and rotations of Lyndon words. The only one for which there's an implication is that non-prime implies that every column has an identical twin (or maybe even more identical siblings) and hence has property X. – Peter Taylor – 2014-11-09T17:16:21.180

OK, thanks. I meant to say "the absence of property X implies Lyndon wordity" since you've stated the contrapositive: that the absence of Lyndon wordity implies property X. – COTO – 2014-11-09T17:36:48.667

@COTO, I haven't stated that, and it's not true. It would imply that every matrix has property X. – Peter Taylor – 2014-11-09T17:37:57.633

"[W]e can consider only circulant matrices in which the first row is a Lyndon word". Ergo if the first row is not a Lyndon word, we cannot consider the resulting circulant matrix. Ergo the resulting circulant matrix must have property X, otherwise we would have to consider it. What am I getting confused? – COTO – 2014-11-09T17:51:03.103

@COTO, "otherwise we can rotate". I'm partitioning the matrices into equivalence classes and only testing one representative from each equivalence class. – Peter Taylor – 2014-11-09T17:54:32.993

Ah, OK. I read something different out of your sentence structure, but I get it now. – COTO – 2014-11-09T17:59:32.920

Wah, 28/15, a good result! I'm curious to see the detailed results for each n. Is there a pattern that can be observed in what kind of Lyndon words satisfy the requirement? If we can find any, then we might be able to find a constructive proof. – justhalf – 2014-11-10T01:26:41.867

@PeterTaylor I've been trying to understand your revolvingDoorDifference function and I can say that that's some high level sorcery :) Is that your own implementation or did you find it somewhere? – Suboptimus Prime – 2014-11-11T16:26:56.713

1

@SuboptimusPrime, I found it at http://people.math.sfu.ca/~kya17/teaching/math343/16-343.pdf and fixed a bug. Interestingly the algorithm I'm now using to iterate through Lyndon words is one of a class of related algorithms which also does k-of-n subsets, so I might be able to refactor and share some code.

– Peter Taylor – 2014-11-11T16:36:08.310

@PeterTaylor Can I borrow your revolving door implementation? It's seems better than the algorithm I'm currently using. – Suboptimus Prime – 2014-11-11T16:47:29.793

@SuboptimusPrime, go for it. We're already both stealing from each other, no reason to stop now ;) – Peter Taylor – 2014-11-11T17:12:09.597

@PeterTaylor Cool, thanks. Calculating the sum of columns on each iteration was the second biggest bottleneck for my program. That still leaves the biggest bottleneck which is the hashset performance, even though my hashset is almost twice as fast as the default one. – Suboptimus Prime – 2014-11-11T17:30:52.727

9

Python 2 - 21/12

In the process of proving that a 2-(3/n) always exists for any n

Inspired by this question, I used De Bruijn Sequence to brute force the possible matrices. And after bruteforcing for n=6,7,8,9,10, I found a pattern that the highest solution is always in the shape of (n, 2n-3).

So I created another method to bruteforce just that shape of matrix, and use multiprocessing to speed things up, since this task is highly distributable. In 16-core ubuntu, it can find solution for n=12 in around 4 minutes:

Trying (0, 254)
Trying (254, 509)
Trying (509, 764)
Trying (764, 1018)
Trying (1018, 1273)
Trying (1273, 1528)
Trying (1528, 1782)
Trying (1782, 2037)
Trying (2037, 2292)
Trying (2292, 2546)
Trying (2546, 2801)
Trying (2801, 3056)
Trying (3056, 3310)
Trying (3820, 4075)
Trying (3565, 3820)
Trying (3310, 3565)
(1625, 1646)
[[0 0 0 1 0 0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0]
 [0 0 1 0 0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0]
 [0 1 0 0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0]
 [1 0 0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0]
 [0 0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 1]
 [0 1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0]
 [1 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0]
 [0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1]
 [1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0]
 [1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 1]
 [1 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 1 1]
 [1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 1 0 1 1 1]]
(12, 21)
Score: 1.7500

real    4m9.121s
user    42m47.472s
sys 0m5.780s

The bulk of computation goes to checking property X, which requires checking all subsets (there are 2^(2n-3) subsets)

Note that I rotate the first row to the left, not to the right as in the question. But these are equivalent since you can just reverse the whole matrix. =)

The code:

import math
import numpy as np
from itertools import combinations
from multiprocessing import Process, Queue, cpu_count

def de_bruijn(k, n):
    """
    De Bruijn sequence for alphabet k
    and subsequences of length n.
    """
    alphabet = list(range(k))
    a = [0] * k * n
    sequence = []
    def db(t, p):
        if t > n:
            if n % p == 0:
                for j in range(1, p + 1):
                    sequence.append(a[j])
        else:
            a[t] = a[t - p]
            db(t + 1, p)
            for j in range(a[t - p] + 1, k):
                a[t] = j
                db(t + 1, t)
    db(1, 1)
    return sequence

def generate_cyclic_matrix(seq, n):
    result = []
    for i in range(n):
        result.append(seq[i:]+seq[:i])
    return np.array(result)

def generate_cyclic_matrix_without_property_x(n=3, n_jobs=-1):
    seq = de_bruijn(2,n)
    seq = seq + seq[:n/2]
    max_idx = len(seq)
    max_score = 1
    max_matrix = np.array([[]])
    max_ij = (0,0)
    workers = []
    queue = Queue()
    if n_jobs < 0:
        n_jobs += cpu_count()+1
    for i in range(n_jobs):
        worker = Process(target=worker_function, args=(seq,i*(2**n-2*n+3)/n_jobs, (i+1)*(2**n-2*n+3)/n_jobs, n, queue))
        workers.append(worker)
        worker.start()
    (result, max_ij) = queue.get()
    for worker in workers:
        worker.terminate()
    return (result, max_ij)

def worker_function(seq,min_idx,max_idx,n,queue):
    print 'Trying (%d, %d)' % (min_idx, max_idx)
    for i in range(min_idx, max_idx):
        j = i+2*n-3
        result = generate_cyclic_matrix(seq[i:j], n)
        if has_property_x(result):
            continue
        else:
            queue.put( (result, (i,j)) )
            return

def has_property_x(mat):
    vecs = zip(*mat)
    vector_sums = set()
    for i in range(1, len(vecs)+1):
        for combination in combinations(vecs, i):
            vector_sum = tuple(sum(combination, np.array([0]*len(mat))))
            if vector_sum in vector_sums:
                return True
            else:
                vector_sums.add(vector_sum)
    return False

def main():
    import sys
    n = int(sys.argv[1])
    if len(sys.argv) > 2:
        n_jobs = int(sys.argv[2])
    else:
        n_jobs = -1
    (matrix, ij) = generate_cyclic_matrix_without_property_x(n, n_jobs)
    print ij
    print matrix
    print matrix.shape
    print 'Score: %.4f' % (float(matrix.shape[1])/matrix.shape[0])

if __name__ == '__main__':
    main()

Old answer, for reference

The optimal solution so far (n=10):

(855, 872)
[[1 1 0 1 0 1 0 0 1 1 1 1 0 1 1 1 0]
 [1 0 1 0 1 0 0 1 1 1 1 0 1 1 1 0 1]
 [0 1 0 1 0 0 1 1 1 1 0 1 1 1 0 1 1]
 [1 0 1 0 0 1 1 1 1 0 1 1 1 0 1 1 0]
 [0 1 0 0 1 1 1 1 0 1 1 1 0 1 1 0 1]
 [1 0 0 1 1 1 1 0 1 1 1 0 1 1 0 1 0]
 [0 0 1 1 1 1 0 1 1 1 0 1 1 0 1 0 1]
 [0 1 1 1 1 0 1 1 1 0 1 1 0 1 0 1 0]
 [1 1 1 1 0 1 1 1 0 1 1 0 1 0 1 0 0]
 [1 1 1 0 1 1 1 0 1 1 0 1 0 1 0 0 1]]
(10, 17)
Score: 1.7000

For n=7:

(86, 97)
[[0 1 1 1 0 1 0 0 1 1 1]
 [1 1 1 0 1 0 0 1 1 1 0]
 [1 1 0 1 0 0 1 1 1 0 1]
 [1 0 1 0 0 1 1 1 0 1 1]
 [0 1 0 0 1 1 1 0 1 1 1]
 [1 0 0 1 1 1 0 1 1 1 0]
 [0 0 1 1 1 0 1 1 1 0 1]]
(7, 11)
Score: 1.5714

A solution with the shape as described by OP (n=8):

(227, 239)
[[0 1 0 1 1 1 1 1 0 1 1 0]
 [1 0 1 1 1 1 1 0 1 1 0 0]
 [0 1 1 1 1 1 0 1 1 0 0 1]
 [1 1 1 1 1 0 1 1 0 0 1 0]
 [1 1 1 1 0 1 1 0 0 1 0 1]
 [1 1 1 0 1 1 0 0 1 0 1 1]
 [1 1 0 1 1 0 0 1 0 1 1 1]
 [1 0 1 1 0 0 1 0 1 1 1 1]]
(8, 12)
Score: 1.5000

But a better one (n=8):

(95, 108)
[[0 1 1 0 0 1 0 0 0 1 1 0 1]
 [1 1 0 0 1 0 0 0 1 1 0 1 0]
 [1 0 0 1 0 0 0 1 1 0 1 0 1]
 [0 0 1 0 0 0 1 1 0 1 0 1 1]
 [0 1 0 0 0 1 1 0 1 0 1 1 0]
 [1 0 0 0 1 1 0 1 0 1 1 0 0]
 [0 0 0 1 1 0 1 0 1 1 0 0 1]
 [0 0 1 1 0 1 0 1 1 0 0 1 0]]
(8, 13)
Score: 1.6250

It also found another optimal solution at n=9:

(103, 118)
[[0 1 0 1 1 1 0 0 0 0 1 1 0 0 1]
 [1 0 1 1 1 0 0 0 0 1 1 0 0 1 0]
 [0 1 1 1 0 0 0 0 1 1 0 0 1 0 1]
 [1 1 1 0 0 0 0 1 1 0 0 1 0 1 0]
 [1 1 0 0 0 0 1 1 0 0 1 0 1 0 1]
 [1 0 0 0 0 1 1 0 0 1 0 1 0 1 1]
 [0 0 0 0 1 1 0 0 1 0 1 0 1 1 1]
 [0 0 0 1 1 0 0 1 0 1 0 1 1 1 0]
 [0 0 1 1 0 0 1 0 1 0 1 1 1 0 0]]
(9, 15)
Score: 1.6667

The code is as follows. It's just brute force, but at least it can find something better than OP's claim =)

import numpy as np
from itertools import combinations

def de_bruijn(k, n):
    """
    De Bruijn sequence for alphabet k
    and subsequences of length n.
    """
    alphabet = list(range(k))
    a = [0] * k * n
    sequence = []
    def db(t, p):
        if t > n:
            if n % p == 0:
                for j in range(1, p + 1):
                    sequence.append(a[j])
        else:
            a[t] = a[t - p]
            db(t + 1, p)
            for j in range(a[t - p] + 1, k):
                a[t] = j
                db(t + 1, t)
    db(1, 1)
    return sequence

def generate_cyclic_matrix(seq, n):
    result = []
    for i in range(n):
        result.append(seq[i:]+seq[:i])
    return np.array(result)

def generate_cyclic_matrix_without_property_x(n=3):
    seq = de_bruijn(2,n)
    max_score = 0
    max_matrix = []
    max_ij = (0,0)
    for i in range(2**n+1):
        for j in range(i+n, 2**n+1):
            score = float(j-i)/n
            if score <= max_score:
                continue
            result = generate_cyclic_matrix(seq[i:j], n)
            if has_property_x(result):
                continue
            else:
                if score > max_score:
                    max_score = score
                    max_matrix = result
                    max_ij = (i,j)
    return (max_matrix, max_ij)

def has_property_x(mat):
    vecs = zip(*mat)
    vector_sums = set()
    for i in range(1, len(vecs)):
        for combination in combinations(vecs, i):
            vector_sum = tuple(sum(combination, np.array([0]*len(mat))))
            if vector_sum in vector_sums:
                return True
            else:
                vector_sums.add(vector_sum)
    return False

def main():
    import sys
    n = int(sys.argv[1])
    (matrix, ij) = generate_cyclic_matrix_without_property_x(n)
    print ij
    print matrix
    print matrix.shape
    print 'Score: %.4f' % (float(matrix.shape[1])/matrix.shape[0])

if __name__ == '__main__':
    main()

justhalf

Posted 2014-11-05T07:19:32.653

Reputation: 2 070

A great start :) – None – 2014-11-05T10:20:11.530

2@Lembik Now I can beat almost (limited by computational time) anyone claiming any score below 2. =) – justhalf – 2014-11-05T12:53:58.203

In that case, can you beat 19/10 ? – None – 2014-11-05T14:59:08.753

@Lembik I don't think I can. It requires n >= 31, which implies I'd need to check up to 2^(2n-3) = 2^59 combinations of 31-dimensional vector. Won't finish in our lifetime =D – justhalf – 2014-11-05T15:03:42.147

I'm not sure where de Bruijn sequences fit into things - can't you just brute-force by counting? But your brute force doesn't actually seem to be complete, because I think I've got a better result than the one you found. – Peter Taylor – 2014-11-05T17:12:55.260

2Can you prove that you can always obtain a matrix of n*(2n-3) – xnor – 2014-11-05T19:00:10.247

At least the formula doesn't work for n=5: The best ratio is 6/5 or 7/6. – kennytm – 2014-11-05T20:08:26.610

@xnor: No, I haven't been able to prove it, unfortunately. This is just a conjecture after seeing a pattern, and the code is just to brute-force it. – justhalf – 2014-11-06T00:13:45.277

@PeterTaylor: I used de Bruijn sequence because I initially thought that the vectors of all possible combinations don't have property X (which is why I thought the answer is infinite). But after I found that that's not the case, I still continued using de Bruijn sequence. You seem to have found the truly optimal solution, though. Congrats! =). And btw using this sequence I only need to brute-force 2^n instead of 2^(2n-3), although like you said, it's not complete. – justhalf – 2014-11-06T00:15:20.223

@KennyTM: Yes, the formula only works for n>=6. For n<=5, we have: 1/1, 2/2, 4/3, 5/4, 6/5. – justhalf – 2014-11-06T00:16:28.523

Since the score is smaller for low n (6/5 instead of 7/5), and given the Java answer is correct it is higher for big n (20/11 instead of 19/11, 22/12 instead of 21/12), I believe the matrix width is superlinear in n, and the score may be able to exceed 2. – kennytm – 2014-11-06T07:46:49.190

@KennyTM: I can't really follow your reasoning. What do you mean by "higher for big n"? My answer also gets higher for bigger n, just that its limit is 2. – justhalf – 2014-11-06T07:53:22.653

@justhalf I mean the Java program shows for higher n, the pattern becomes 2n-2. Perhaps it becomes 2n-1, 2n, ... 2n+x for even higher n, which breaks the limit "2". – kennytm – 2014-11-06T08:06:22.907

I see. Let's see whether that's true. =). Btw @xnor, I found another answer that seems to suggest there will always be solution for 2n-3: http://math.stackexchange.com/a/672555/92861

– justhalf – 2014-11-06T08:08:56.420

Oops, yeah, that's true, sorry for that. I edited my answer. – justhalf – 2014-11-07T23:50:41.230

7

24/13 26/14 28/15 30/16 32/17 (C#)

Edit: Deleted outdated info from my answer. I'm using mostly the same algorithm as Peter Taylor (Edit: looks like he is using a better algorithm now), although I've added some of my own optimizations:

  • I've implemented "meet in the middle" strategy of searching for column sets with the same vector sum (suggested by this KennyTM's comment). This strategy improved memory usage a lot, but it's rather slow, so I've added the HasPropertyXFast function, that quickly checks if there are small set with equal sums before using "meet in the middle" approach.
  • While iterating through column sets in the HasPropertyXFast function, I start from checking column sets with 1 column, then with 2, 3 and so on. The function returns as soon as the first collision of column sums is found. In practice it means that I usually have to check just a few hundreds or thousands of column sets rather than millions.
  • I'm using long variables to store and compare entire columns and their vector sums. This approach is at least an order of magnitude faster than comparing columns as arrays.
  • I've added my own implementation of hashset, optimized for long data type and for my usage patterns.
  • I'm reusing the same 3 hashsets for the entire lifetime of the application to reduce the number of memory allocations and improve performance.
  • Multithreading support.

Program output:

00000000000111011101010010011111
10000000000011101110101001001111
11000000000001110111010100100111
11100000000000111011101010010011
11110000000000011101110101001001
11111000000000001110111010100100
01111100000000000111011101010010
00111110000000000011101110101001
10011111000000000001110111010100
01001111100000000000111011101010
00100111110000000000011101110101
10010011111000000000001110111010
01001001111100000000000111011101
10100100111110000000000011101110
01010010011111000000000001110111
10101001001111100000000000111011
11010100100111110000000000011101
Score: 32/17 = 1,88235294117647
Time elapsed: 02:11:05.9791250

Code:

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

class Program
{
    const int MaxWidth = 32;
    const int MaxHeight = 17;

    static object _lyndonWordLock = new object();

    static void Main(string[] args)
    {
        Stopwatch sw = Stopwatch.StartNew();
        double maxScore = 0;
        const int minHeight = 17; // 1
        for (int height = minHeight; height <= MaxHeight; height++)
        {
            Console.WriteLine("Row count = " + height);
            Console.WriteLine("Time elapsed: " + sw.Elapsed + "\r\n");

            int minWidth = Math.Max(height, (int)(height * maxScore) + 1);
            for (int width = minWidth; width <= MaxWidth; width++)
            {
#if MULTITHREADING
                int[,] matrix = FindMatrixParallel(width, height);
#else
                int[,] matrix = FindMatrix(width, height);
#endif
                if (matrix != null)
                {
                    PrintMatrix(matrix);
                    Console.WriteLine("Time elapsed: " + sw.Elapsed + "\r\n");
                    maxScore = (double)width / height;
                }
                else
                    break;
            }
        }
    }

#if MULTITHREADING
    static int[,] FindMatrixParallel(int width, int height)
    {
        _lyndonWord = 0;
        _stopSearch = false;

        int threadCount = Environment.ProcessorCount;
        Task<int[,]>[] tasks = new Task<int[,]>[threadCount];
        for (int i = 0; i < threadCount; i++)
            tasks[i] = Task<int[,]>.Run(() => FindMatrix(width, height));

        int index = Task.WaitAny(tasks);
        if (tasks[index].Result != null)
            _stopSearch = true;

        Task.WaitAll(tasks);
        foreach (Task<int[,]> task in tasks)
            if (task.Result != null)
                return task.Result;

        return null;
    }

    static volatile bool _stopSearch;
#endif

    static int[,] FindMatrix(int width, int height)
    {
#if MULTITHREADING
        _columnSums = new LongSet();
        _left = new LongSet();
        _right = new LongSet();
#endif

        foreach (long rowTemplate in GetLyndonWords(width))
        {
            int[,] matrix = new int[width, height];
            for (int x = 0; x < width; x++)
            {
                int cellValue = (int)(rowTemplate >> (width - 1 - x)) % 2;
                for (int y = 0; y < height; y++)
                    matrix[(x + y) % width, y] = cellValue;
            }

            if (!HasPropertyX(matrix))
                return matrix;

#if MULTITHREADING
            if (_stopSearch)
                return null;
#endif
        }

        return null;
    }

#if MULTITHREADING
    static long _lyndonWord;
#endif

    static IEnumerable<long> GetLyndonWords(int length)
    {
        long lyndonWord = 0;
        long max = (1L << (length - 1)) - 1;
        while (lyndonWord <= max)
        {
            if ((lyndonWord % 2 != 0) && PrecedesReversal(lyndonWord, length))
                yield return lyndonWord;

#if MULTITHREADING
            lock (_lyndonWordLock)
            {
                if (_lyndonWord <= max)
                    _lyndonWord = NextLyndonWord(_lyndonWord, length);
                else
                    yield break;

                lyndonWord = _lyndonWord;
            }
#else
            lyndonWord = NextLyndonWord(lyndonWord, length);
#endif
        }
    }

    static readonly int[] _lookup =
    {
        32, 0, 1, 26, 2, 23, 27, 0, 3, 16, 24, 30, 28, 11, 0, 13, 4, 7, 17,
        0, 25, 22, 31, 15, 29, 10, 12, 6, 0, 21, 14, 9, 5, 20, 8, 19, 18
    };

    static int NumberOfTrailingZeros(uint i)
    {
        return _lookup[(i & -i) % 37];
    }

    static long NextLyndonWord(long w, int length)
    {
        if (w == 0)
            return 1;

        int currentLength = length - NumberOfTrailingZeros((uint)w);
        while (currentLength < length)
        {
            w += w >> currentLength;
            currentLength *= 2;
        }

        w++;

        return w;
    }

    private static bool PrecedesReversal(long lyndonWord, int length)
    {
        int shift = length - 1;

        long reverse = 0;
        for (int i = 0; i < length; i++)
        {
            long bit = (lyndonWord >> i) % 2;
            reverse |= bit << (shift - i);
        }

        for (int i = 0; i < length; i++)
        {
            if (reverse < lyndonWord)
                return false;

            long bit = reverse % 2;
            reverse /= 2;
            reverse += bit << shift;
        }

        return true;
    }

#if MULTITHREADING
    [ThreadStatic]
#endif
    static LongSet _left = new LongSet();
#if MULTITHREADING
    [ThreadStatic]
#endif
    static LongSet _right = new LongSet();

    static bool HasPropertyX(int[,] matrix)
    {
        long[] matrixColumns = GetMatrixColumns(matrix);
        if (matrixColumns.Length == 1)
            return false;

        return HasPropertyXFast(matrixColumns) || MeetInTheMiddle(matrixColumns);
    }

    static bool MeetInTheMiddle(long[] matrixColumns)
    {
        long[] leftColumns = matrixColumns.Take(matrixColumns.Length / 2).ToArray();
        long[] rightColumns = matrixColumns.Skip(matrixColumns.Length / 2).ToArray();

        if (PrepareHashSet(leftColumns, _left) || PrepareHashSet(rightColumns, _right))
            return true;

        foreach (long columnSum in _left.GetValues())
            if (_right.Contains(columnSum))
                return true;

        return false;
    }

    static bool PrepareHashSet(long[] columns, LongSet sums)
    {
        int setSize = (int)System.Numerics.BigInteger.Pow(3, columns.Length);
        sums.Reset(setSize, setSize);
        foreach (long column in columns)
        {
            foreach (long sum in sums.GetValues())
                if (!sums.Add(sum + column) || !sums.Add(sum - column))
                    return true;

            if (!sums.Add(column) || !sums.Add(-column))
                return true;
        }

        return false;
    }

#if MULTITHREADING
    [ThreadStatic]
#endif
    static LongSet _columnSums = new LongSet();

    static bool HasPropertyXFast(long[] matrixColumns)
    {
        int width = matrixColumns.Length;

        int maxColumnCount = width / 3;
        _columnSums.Reset(width, SumOfBinomialCoefficients(width, maxColumnCount));

        int resetBit, setBit;
        for (int k = 1; k <= maxColumnCount; k++)
        {
            uint columnMask = (1u << k) - 1;
            long sum = 0;
            for (int i = 0; i < k; i++)
                sum += matrixColumns[i];

            while (true)
            {
                if (!_columnSums.Add(sum))
                    return true;
                if (!NextColumnMask(columnMask, k, width, out resetBit, out setBit))
                    break;
                columnMask ^= (1u << resetBit) ^ (1u << setBit);
                sum = sum - matrixColumns[resetBit] + matrixColumns[setBit];
            }
        }

        return false;
    }

    // stolen from Peter Taylor
    static bool NextColumnMask(uint mask, int k, int n, out int resetBit, out int setBit)
    {
        int gap = NumberOfTrailingZeros(~mask);
        int next = 1 + NumberOfTrailingZeros(mask & (mask + 1));

        if (((k - gap) & 1) == 0)
        {
            if (gap == 0)
            {
                resetBit = next - 1;
                setBit = next - 2;
            }
            else if (gap == 1)
            {
                resetBit = 0;
                setBit = 1;
            }
            else
            {
                resetBit = gap - 2;
                setBit = gap;
            }
        }
        else
        {
            if (next == n)
            {
                resetBit = 0;
                setBit = 0;
                return false;
            }

            if ((mask & (1 << next)) == 0)
            {
                if (gap == 0)
                {
                    resetBit = next - 1;
                    setBit = next;
                }
                else
                {
                    resetBit = gap - 1;
                    setBit = next;
                }
            }
            else
            {
                resetBit = next;
                setBit = gap;
            }
        }

        return true;
    }

    static long[] GetMatrixColumns(int[,] matrix)
    {
        int width = matrix.GetLength(0);
        int height = matrix.GetLength(1);

        long[] result = new long[width];
        for (int x = 0; x < width; x++)
        {
            long column = 0;
            for (int y = 0; y < height; y++)
            {
                column *= 13;
                if (matrix[x, y] == 1)
                    column++;
            }

            result[x] = column;
        }

        return result;
    }

    static int SumOfBinomialCoefficients(int n, int k)
    {
        int result = 0;
        for (int i = 0; i <= k; i++)
            result += BinomialCoefficient(n, i);
        return result;
    }

    static int BinomialCoefficient(int n, int k)
    {
        long result = 1;
        for (int i = n - k + 1; i <= n; i++)
            result *= i;
        for (int i = 2; i <= k; i++)
            result /= i;
        return (int)result;
    }

    static void PrintMatrix(int[,] matrix)
    {
        int width = matrix.GetLength(0);
        int height = matrix.GetLength(1);

        for (int y = 0; y < height; y++)
        {
            for (int x = 0; x < width; x++)
                Console.Write(matrix[x, y]);
            Console.WriteLine();
        }

        Console.WriteLine("Score: {0}/{1} = {2}", width, height, (double)width / height);
    }
}


class LongSet
{
    private static readonly int[] primes =
    {
        17, 37, 67, 89, 113, 149, 191, 239, 307, 389, 487, 613, 769, 967, 1213, 1523, 1907,
        2389, 2999, 3761, 4703, 5879, 7349, 9187, 11489, 14369, 17971, 22469, 28087, 35111,
        43889, 54869, 68597, 85751, 107197, 133999, 167521, 209431, 261791, 327247, 409063,
        511333, 639167, 798961, 998717, 1248407, 1560511, 1950643, 2438309, 3047909,
        809891, 4762367, 5952959, 7441219, 9301529, 11626913, 14533661, 18167089, 22708867,
        28386089, 35482627, 44353297, 55441637, 69302071, 86627603, 108284507, 135355669,
        169194593, 211493263, 264366593, 330458263, 413072843, 516341057, 645426329,
        806782913, 1008478649, 1260598321
    };

    private int[] _buckets;
    private int[] _nextItemIndexes;
    private long[] _items;
    private int _count;
    private int _minCapacity;
    private int _maxCapacity;
    private int _currentCapacity;

    public LongSet()
    {
        Initialize(0, 0);
    }

    private int GetPrime(int capacity)
    {
        foreach (int prime in primes)
            if (prime >= capacity)
                return prime;

        return int.MaxValue;
    }

    public void Reset(int minCapacity, int maxCapacity)
    {
        if (maxCapacity > _maxCapacity)
            Initialize(minCapacity, maxCapacity);
        else
            ClearBuckets();
    }

    private void Initialize(int minCapacity, int maxCapacity)
    {
        _minCapacity = GetPrime(minCapacity);
        _maxCapacity = GetPrime(maxCapacity);
        _currentCapacity = _minCapacity;

        _buckets = new int[_maxCapacity];
        _nextItemIndexes = new int[_maxCapacity];
        _items = new long[_maxCapacity];
        _count = 0;
    }

    private void ClearBuckets()
    {
        Array.Clear(_buckets, 0, _currentCapacity);
        _count = 0;
        _currentCapacity = _minCapacity;
    }

    public bool Add(long value)
    {
        int bucket = (int)((ulong)value % (ulong)_currentCapacity);
        for (int i = _buckets[bucket] - 1; i >= 0; i = _nextItemIndexes[i])
            if (_items[i] == value)
                return false;

        if (_count == _currentCapacity)
        {
            Grow();
            bucket = (int)((ulong)value % (ulong)_currentCapacity);
        }

        int index = _count;
        _items[index] = value;
        _nextItemIndexes[index] = _buckets[bucket] - 1;
        _buckets[bucket] = index + 1;
        _count++;

        return true;
    }

    private void Grow()
    {
        Array.Clear(_buckets, 0, _currentCapacity);

        const int growthFactor = 8;
        int newCapacity = GetPrime(_currentCapacity * growthFactor);
        if (newCapacity > _maxCapacity)
            newCapacity = _maxCapacity;
        _currentCapacity = newCapacity;

        for (int i = 0; i < _count; i++)
        {
            int bucket = (int)((ulong)_items[i] % (ulong)newCapacity);
            _nextItemIndexes[i] = _buckets[bucket] - 1;
            _buckets[bucket] = i + 1;
        }
    }

    public bool Contains(long value)
    {
        int bucket = (int)((ulong)value % (ulong)_buckets.Length);
        for (int i = _buckets[bucket] - 1; i >= 0; i = _nextItemIndexes[i])
            if (_items[i] == value)
                return true;

        return false;
    }

    public IReadOnlyList<long> GetValues()
    {
        return new ArraySegment<long>(_items, 0, _count);
    }
}

Configuration file:

<?xml version="1.0" encoding="utf-8" ?>
<configuration>
  <runtime>
    <gcAllowVeryLargeObjects enabled="true" />
  </runtime>
</configuration>

Suboptimus Prime

Posted 2014-11-05T07:19:32.653

Reputation: 403

In some regards you seem to have pessimised rather than optimised. The only thing which really looks like an optimisation is allowing bits to clash by using ulong and letting the shift wrap instead of using BigInteger. – Peter Taylor – 2014-11-06T17:06:24.327

@PeterTaylor The most important optimisation is in HasPropertyX function. The function returns as soon as the first collision of column sums is found (unlike your scoreLyndonWord function). I've also sorted the column masks in such a way that we first check the column sets that are more likely to collide. These two optimisations improved the performance by an order of magnitude. – Suboptimus Prime – 2014-11-06T17:24:06.407

Although performance changes are often surprising, in principle the early abort shouldn't give more than a factor of 2, and GetSumOfColumns adds an extra loop which I would expect to cost more than that factor of 2. The mask sorting sounds interesting: maybe you could edit the answer to talk a bit about it? (And at some point I will experiment with an alternative way to do the early abort: the reason I can't do it is that HashSet doesn't support concurrent iteration and modification, but I have ideas to avoid the need for an iterator). – Peter Taylor – 2014-11-06T17:31:28.790

@PeterTaylor Added some explanation to the answer, hopefully it clears things up a bit. – Suboptimus Prime – 2014-11-06T18:33:21.023

@SuboptimusPrime: Hmm, yeah, the summing based on columns is basically what I do in my answer as well. But you used faster language (C# > Python in speed) and faster algorithm (Peter's > mine). I saw in another thread that you can use Gray Code as the mask to do incremental sum, so you only need to add/remove one vector each iteration. Perhaps you want to implement that? EDIT: ah, but your limit is the memory, so I guess 26 columns is the maximum for your algorithm so far? – justhalf – 2014-11-07T02:08:05.400

@justhalf, I've done Grey code, but it doesn't seem to make enough of a difference. Well, overall with some other optimisations I have a 7.5-fold improvement for MAX_N=22, but I think the biggest improvement comes from finding a better way to not use BigInteger. – Peter Taylor – 2014-11-07T07:40:48.867

@justhalf Using Gray code means that you have to iterate through column sets in a different order. Even though each iteration is faster, you end up doing way more iterations before you find two column sets with the same vector sum. I've tried it and it was about 20 times slower than my current algorithm. – Suboptimus Prime – 2014-11-07T16:16:40.367

@justhalf Btw, the maximum for my algorithm is actually 27 columns, but there are no 27x14 matrices without property X, so I can't improve my score without optimizing the memory usage of my program. – Suboptimus Prime – 2014-11-07T16:23:22.083

2@justhalf, using a Gray-esque approach for iterating over the subsets of a fixed size is actually worthwhile. It allowed me to find a 26/14 in under 9 minutes and 34 of them in two hours, at which point I aborted. Currently testing to see whether I can get 28/15 in a reasonable time. – Peter Taylor – 2014-11-07T18:00:55.510

30/16 is great! Seems we are creeping to 2.. but can we actually get to it? Do you think 16 is optimal for 30, 15 is optimal for 28 etc.? – None – 2014-11-09T19:36:18.707

@Lembik I'm pretty sure that 15 rows is optimal for 28 columns although there could be a 29/15 matrix without property X. I intentionally skipped 29/15 and went straight for 30/16 since 29/15 was taking too long. – Suboptimus Prime – 2014-11-09T19:45:26.950

@SuboptimusPrime It be great to find any ratio that contradicts the hypothesis that increasing the number or columns by 2 always increases the number of rows by 1. If the limit of the score really is 2 I feel that a) would be very surprising and b) ought to have a proof. How long does each matrix take to test in the 30/16 case? – None – 2014-11-09T19:48:49.427

@Lembik I did't really check how long it takes to test a single matrix but it's probably just a tiny fraction of a second. But it took over 3 hours to find a 30/16 matrix that doesn't have property X since the number of matrices grows exponentially as number of columns increases. – Suboptimus Prime – 2014-11-09T20:09:55.020

If each test is quick and the memory problems are ok, it might be worth trying a few thousand random 40/20 matrices where the first row has roughly equal numbers of 1s and 0s. – None – 2014-11-09T20:11:41.523

@Lembik I suspect that there is like a billion of 40/20 matrices that do have property X for each one that doesn't, so the probability of finding one is not too high. This is just a guess, but I can't check it anyway since the current version of my program can not work with 40/20 matrices. – Suboptimus Prime – 2014-11-09T20:28:15.567

@Lembik, my code takes about 20 seconds to check a 28/15. I'm searching for 29/15, currently at 35 hours with no success. I'm also searching for 30/16, with some additional heuristics to try to speed it up, but no hits yet in just under 6 hours. – Peter Taylor – 2014-11-09T23:18:38.580

Are you able to explore 31/16 and/or 32/16 exhaustively? – None – 2014-11-11T20:55:29.290

@Lembik Probably not. I've tried to search 31/16 matrices for 10 hours and didn't find anything. Exhaustive search would probably take much longer because even 29/15 didn't finish in 10 hours. – Suboptimus Prime – 2014-11-12T05:11:00.320

1@Lembik, I explored 29/15 exhaustively in 75.5 hours. 31/16 would take about 3 times as long - so more than a week. Both of us have made some optimisations since I started running that test of 29/15, so maybe it would be down to a week now. There's nothing stopping you from compiling my code or SuboptimusPrime's code and running it yourself if you have a computer which you can leave on for that long. – Peter Taylor – 2014-11-12T10:47:09.300