Matrix property X revisited (or the Joy of X)

11

4

This challenge is partly an algorithms challenge, partly an optimization challenge and partly simply a fastest code challenge.

A T matrix is fully specified by its first row r and first column c. Each remaining element of the matrix is just a copy of the element that is diagonally up and left. That is M[i,j] = M[i-1,j-1]. We will allow T matrices which are not square. 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 T 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 one or more columns is simply an element-wise summation of their columns. That is the sum of two or more columns containing x elements each is another column containing x elements. The sum of one column is trivially the column itself.

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 T matrix with binary entries and which does not have property X. For clarity, a matrix with binary entries has the property that each one of its entries is either 0 or 1.

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

All the answers at Find highest scoring matrix without property X are valid here but they are not optimal. There are T matrices without property X which are not cyclic.

For example, there is a 7 by 12 T matrix without property X but no such cyclic matrix.

21/11 would beat all the current answers from this and the previous challenge.

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.

Bonus The first answer with a score greater than 2 gets an immediate 200 point bounty award. Ton Hospel has now achieved this!


Current leader board

  • C++. Score 31/15 by Ton Hospel
  • Java. Score 36/19 by Peter Taylor
  • Haskell. Score 14/8 by alexander-brett

user9206

Posted 2015-04-26T10:01:42.697

Reputation:

By the " two non-empty sets of columns with non-identical indices" you mean two sets of columns that are disjoint? Or, to rephrase this, is {1 ,3}, {1, 5} a valid two subsets of columns? – pawel.boczarski – 2015-04-26T10:52:26.117

@pawel.boczarski No not disjoint. Just non-identical. So {1 ,3}, {1, 5} is valid. – None – 2015-04-26T10:53:40.757

Ok. What about M[i,1] - is it a "borrow" from the last column of M[i-1] (zero is not a valid matrix column index)? And actually this is "up and left" rather than "up and right". – pawel.boczarski – 2015-04-26T10:58:23.620

@pawel.boczarski "right" was a typo. Thanks. The first row and column can be set to anything you like. They define the whole rest of the matrix. Does that answer your question? – None – 2015-04-26T13:13:55.537

Ok, I got it. It was my fault to not have read carefully that the first column is also defined. – pawel.boczarski – 2015-04-26T14:28:52.827

For the bonus, do you mean the *first* answer with a score greater than 2 gets 100 bounty? Otherwise that might cost you... :) – trichoplax – 2015-04-26T18:30:52.867

@trichoplax I do! – None – 2015-04-26T18:31:41.520

@Lembik: (continuing the comments from the now delted answer) As one particular sum of columns you consider the sum of all columns in a particular set S. Do you also allow the empty set S={} as one of those sets in your definition? This would just eliminate all matrices that contian a column that consists only of zeros. If i is the index of the said column, then the matrix would have the property X with the sets {},{i}. – flawr – 2015-04-27T18:32:32.020

@flawr No I don't allow the empty set (see the definition of property X). However we know that a matrix (with 2 or more columns) with a column consisting only of zeros always has property X in any case as that column plus any other column equals that other column. – None – 2015-04-27T18:34:39.173

Ah ok, well then it does not matter whether you allow the empty sum or not, except for the matrix 1x1 matrix 0=) – flawr – 2015-04-27T18:35:57.037

Is there any reason for calling the space of matrices T rather than Toeplitz? – Peter Taylor – 2015-04-28T14:08:21.677

@PeterTaylor No, not really. I was hoping for a pun involving T and X but in the end couldn't come up with a good one. – None – 2015-04-28T14:10:03.220

@Lembik: I'm having trouble fully understanding property X. Previously you state that any matrix containing a column of zeros because it is functionally an identity property. This indicates that the same column can appear in both sets. If I have column A which is some non-zero column, and column B which consists only of zeros, then the Matrix has property X because S1(A) = S2(A+B). So the "non-identical indices" property is met by any difference between the sets? Is the number of columns in each comparison set limited to 2, or can it be larger? – CChak – 2015-04-28T14:17:11.407

@CChak Yes the same column can appear in both sets. The two sets don't have to be disjoint, just not-identical (because that would trivially be uninteresting). The number of columns in each comparison set can be as large as you like. – None – 2015-04-28T14:24:15.540

@Lembik: Thanks for the quick response, but I hit return and it posted before I was done. The real question I was getting to was is there a limit on the number of columns in the summation set? or is it simply 1 or 2? I am assuming that each column index can be used once and only once per set. Again because that would make the problem uninteresting. – CChak – 2015-04-28T14:27:19.560

@CChak I also replied to that I think (after an edit). Is it clear now? – None – 2015-04-28T14:27:59.467

Yep. I've got it. thanks! – CChak – 2015-04-28T14:32:39.560

For the vector sum, since the matrix is binary, is the addition also performed in the binary system? For example, is 1 + 1 == 2 or 1 + 1 == 0? – Reto Koradi – 2015-05-10T21:37:40.390

@RetroKoradi 1+1=2 – None – 2015-05-11T05:51:22.440

Answers

6

C++, Score 23/12 25/13 27/14 28/14 31/15

Finally a result with ratio > 2:

rows=15,cols=31
1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 0 1 1 0 
1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 0 1 1 
1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 0 1 
1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 0 
1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 1 
0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 1 
0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 
1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 1 
0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 0 
0 0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 0 
0 0 0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 0 
1 0 0 0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 1 
0 1 0 0 0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 0 
0 0 1 0 0 0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 0 
1 0 0 1 0 0 0 1 0 0 1 1 1 1 1 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 

I completely explored 1 to 14 rows. 15 would take too long to completely explore. The results are:

1/1   = 1
2/2   = 1
4/3   = 1.333
5/4   = 1.25
7/5   = 1.4
9/6   = 1.5
12/7  = 1.714
14/8  = 1.75
16/9  = 1.778
18/10 = 1.8
20/11 = 1.818
23/12 = 1.917
25/13 = 1.923
28/14 = 2

The code given below is an older version of the program. The newest version is at https://github.com/thospel/personal-propertyX.

/*
  Compile using something like:
    g++ -Wall -O3 -march=native -fstrict-aliasing -std=c++11 -g propertyX.cpp -lpthread -o propertyX
*/
#include <cstdint>
#include <climits>
#include <ctgmath>
#include <iostream>
#include <vector>
#include <array>
#include <chrono>
#include <mutex>
#include <atomic>
#include <thread>

using namespace std;

const int ELEMENTS = 2;

using uint    = unsigned int;
using Element = uint64_t;
using Column  = array<Element, ELEMENTS>;
using Set     = vector<Column>;
using Sum     = uint8_t;
using Index   = uint32_t;
using sec = chrono::seconds;

int const PERIOD = 5*60;
int const MAX_ROWS = 54;
int const COL_FACTOR = (MAX_ROWS+1) | 1;                // 55
int const ROW_ZERO   = COL_FACTOR/2;                    // 27
int const ROWS_PER_ELEMENT = CHAR_BIT * sizeof(Element) / log2(COL_FACTOR); //11
Element constexpr ELEMENT_FILL(Element v = ROW_ZERO, int n = ROWS_PER_ELEMENT) {
    return n ? ELEMENT_FILL(v, n-1) * COL_FACTOR + v : 0;
}
Element constexpr POWN(Element v, int n) {
    return n ? POWN(v, n-1)*v : 1;
}
Element const ELEMENT_TOP = POWN(COL_FACTOR, ROWS_PER_ELEMENT -1);
int const MAX_COLS = ROWS_PER_ELEMENT * ELEMENTS;       // 22

atomic<Index> col_next;
atomic<uint>  period;
chrono::steady_clock::time_point start;
mutex period_mutex;

uint ratio_row;
uint ratio_col;
mutex ratio_mutex;

auto const nr_threads = thread::hardware_concurrency();
// auto const nr_threads = 1;

struct State {
    State(uint cols);
    void process(Index i);
    void extend(uint row);
    void print(uint rows);
    Index nr_columns() const { return static_cast<Index>(1) << cols_; }

    Column last_;
    Element top_;
    int top_offset_;
    uint ratio_row_ = 0;
    uint ratio_col_ = 1;
    uint cols_;
    array<Sum, MAX_ROWS + MAX_COLS -1> side;
    vector<Set> sets;
};

ostream& operator<<(ostream& os, Column const& column) {
    for (int i=0; i<ELEMENTS; ++i) {
        auto v = column[i];
        for (int j=0; j<ROWS_PER_ELEMENT; ++j) {
            auto bit = v / ELEMENT_TOP;
            cout << " " << bit;
            v -= bit * ELEMENT_TOP;
            v *= COL_FACTOR;
        }
    }
    return os;
}

State::State(uint cols) : cols_{cols} {
    sets.resize(MAX_ROWS+2);
    for (int i=0; i<2; ++i) {
        sets[i].resize(2);
        for (int j=0; j < ELEMENTS; ++j) {
            sets[i][0][j] =  ELEMENT_FILL();
            sets[i][1][j] =  static_cast<Element>(-1) - ELEMENT_FILL(1);
        }
    }
    top_ = POWN(COL_FACTOR, (cols_-1) % ROWS_PER_ELEMENT);
    top_offset_ = ELEMENTS - 1 - (cols_-1) / ROWS_PER_ELEMENT;
}

void State::print(uint rows) {
    for (auto c=0U; c<cols_;c++) {
        for (auto r=0U; r<rows;r++) {
            cout << static_cast<int>(side[cols_-c+r-1]) << " ";
        }
        cout << "\n";
    }
    cout << "----------" << endl;
}

void check(uint cols, uint t) {
    State state(cols);

    Index nr_columns = state.nr_columns();
    while (1) {
        Index col = col_next++;
        if (col >= nr_columns) break;
        state.process(col);

        auto now = chrono::steady_clock::now();
        auto elapsed = chrono::duration_cast<sec>(now-start).count();
        if (elapsed >= period) {
            lock_guard<mutex> lock{period_mutex};
            if (elapsed >= period) {
                cout << "col=" << col << "/" << nr_columns << " (" << 100.*col/nr_columns << "% " << elapsed << " s)" << endl;
                period = (elapsed/PERIOD+1)*PERIOD;
            }
        }
    }
}

void State::process(Index col) {
    last_.fill(0);
    for (uint i=0; i<cols_; ++i) {
        Element bit = col >> i & 1;
        side[i] = bit;
        Element carry = 0;
        for (int j=0; j<ELEMENTS; ++j) {
            auto c = last_[j] % COL_FACTOR;
            last_[j] = last_[j] / COL_FACTOR + carry * ELEMENT_TOP;
            if (j == top_offset_ && bit) last_[j] += top_;
            carry = c;
        }
    }
    // cout << "col=" << col << ", value=" << last_ << "\n";
    extend(0);
}

void State::extend(uint row) {
    // cout << "Extend row " << row << " " << static_cast<int>(side[cols_+row-1]) << "\n";
    if (row >= MAX_ROWS) throw(range_error("row out of range"));

    // Execute subset sum. The new column is added to set {from} giving {to}
    // {sum} is the other set.
    auto const& set_from = sets[row];
    auto const& set_sum  = sets[row + 1];
    auto      & set_to   = sets[row + 2];
    if (set_to.size() == 0) {
        auto size = 3 * set_from.size() - 2;
        set_to.resize(size);
        for (int j=0; j<ELEMENTS; ++j)
            set_to[size-1][j] = static_cast<Element>(-1) - ELEMENT_FILL(1);
    }

    // Merge sort {set_from - last_} , {set_from} and {set_from + last_}
    auto ptr_sum    = &set_sum[1][0];
    auto ptr_low    = &set_from[0][0];
    auto ptr_middle = &set_from[0][0];
    auto ptr_high   = &set_from[0][0];
    Column col_low, col_high;
    for (int j=0; j<ELEMENTS; ++j) {
        col_low   [j] = *ptr_low++  - last_[j];
        col_high  [j] = *ptr_high++ + last_[j];
    }

    auto ptr_end = &set_to[set_to.size()-1][0];
    auto ptr_to  = &set_to[0][0];
    while (ptr_to < ptr_end) {
        for (int j=0; j<ELEMENTS; ++j) {
            if (col_low[j] < ptr_middle[j]) goto LOW;
            if (col_low[j] > ptr_middle[j]) goto MIDDLE;
        }
        // low == middle
        // cout << "low == middle\n";
        return;

      LOW:
        // cout << "LOW\n";
        for (int j=0; j<ELEMENTS; ++j) {
            if (col_low[j] < col_high[j]) goto LOW0;
            if (col_low[j] > col_high[j]) goto HIGH0;
        }
        // low == high
        // cout << "low == high\n";
        return;

      MIDDLE:
        // cout << "MIDDLE\n";
        for (int j=0; j<ELEMENTS; ++j) {
            if (ptr_middle[j] < col_high[j]) goto MIDDLE0;
            if (ptr_middle[j] > col_high[j]) goto HIGH0;
        }
        // middle == high
        // cout << "middle == high\n";
        return;

      LOW0:
        // cout << "LOW0\n";
        for (int j=0; j<ELEMENTS; ++j) {
            *ptr_to++  = col_low[j];
            col_low[j] = *ptr_low++ - last_[j];
        }
        goto SUM;

      MIDDLE0:
        // cout << "MIDDLE0\n";
        for (int j=0; j<ELEMENTS; ++j)
            *ptr_to++ = *ptr_middle++;
        goto SUM;

      HIGH0:
        // cout << "HIGH0\n";
        for (int j=0; j<ELEMENTS; ++j) {
            *ptr_to++ = col_high[j];
            col_high[j] = *ptr_high++ + last_[j];
        }
        goto SUM;
      SUM:
        for (int j=-ELEMENTS; j<0; ++j) {
            if (ptr_to[j] > ptr_sum[j]) {
                ptr_sum += ELEMENTS;
                goto SUM;
            }
            if (ptr_to[j] < ptr_sum[j]) goto DONE;
        }
        // sum == to
        for (int j=-ELEMENTS; j<0; ++j)
            if (ptr_to[j] != ELEMENT_FILL()) {
                // sum == to and to != 0
                // cout << "sum == to\n";
                // cout << set_sum[(ptr_sum - &set_sum[0][0])/ELEMENTS-1] << "\n";
                return;
            }
      DONE:;
    }
    // cout << "Wee\n";
    auto row1 = row+1;
    if (0)
        for (uint i=0; i<row1+2; ++i) {
            cout << "Set " << i << "\n";
            auto& set = sets[i];
            for (auto& column: set)
                cout << column << "\n";
        }

    if (row1 * ratio_col_ > ratio_row_ * cols_) {
        ratio_row_ = row1;
        ratio_col_ = cols_;
        lock_guard<mutex> lock{ratio_mutex};

        if (ratio_row_ * ratio_col > ratio_row * ratio_col_) {

            auto now = chrono::steady_clock::now();
            auto elapsed = chrono::duration_cast<sec>(now-start).count();
            cout << "cols=" << cols_ << ",rows=" << row1 << " (" << elapsed << " s)\n";
            print(row1);
            ratio_row = ratio_row_;
            ratio_col = ratio_col_;
        }
    }

    auto last = last_;

    Element carry = 0;
    for (int j=0; j<ELEMENTS; ++j) {
        auto c = last_[j] % COL_FACTOR;
        last_[j] = last_[j] / COL_FACTOR + carry * ELEMENT_TOP;
        carry = c;
    }

    side[cols_+row] = 0;
    extend(row1);

    last_[top_offset_] += top_;
    side[cols_+row] = 1;
    extend(row1);

    last_ = last;
}

void my_main(int argc, char** argv) {
    if (!col_next.is_lock_free()) cout << "col_next is not lock free\n";
    if (!period.  is_lock_free()) cout << "period is not lock free\n";

    int min_col = 2;
    int max_col = MAX_COLS;
    if (argc > 1) {
        min_col = atoi(argv[1]);
        if (min_col < 2)
            throw(range_error("Column must be >= 2"));
        if (min_col > MAX_COLS)
            throw(range_error("Column must be <= " + to_string(MAX_COLS)));
    }
    if (argc > 2) {
        max_col = atoi(argv[2]);
        if (max_col < min_col)
            throw(range_error("Column must be >= " + to_string(min_col)));
        if (max_col > MAX_COLS)
            throw(range_error("Column must be <= " + to_string(MAX_COLS)));
    }

    for (int cols = min_col; cols <= max_col; ++cols) {
        cout << "Trying " << cols << " columns" << endl;
        ratio_row = 0;
        ratio_col = 1;
        col_next = 0;
        period = PERIOD;
        start = chrono::steady_clock::now();
        vector<thread> threads;
        for (auto t = 1U; t < nr_threads; t++)
            threads.emplace_back(check, cols, t);
        check(cols, 0);
        for (auto& thread: threads)
            thread.join();
    }
}

int main(int argc, char** argv) {
    try {
        my_main(argc, argv);
    } catch(exception& e) {
        cerr << "Error: " << e.what() << endl;
        exit(EXIT_FAILURE);
    }
    exit(EXIT_SUCCESS);
}

Ton Hospel

Posted 2015-04-26T10:01:42.697

Reputation: 14 114

This is great. The great mystery is if you can ever get a score of 2. – None – 2016-11-05T14:19:57.307

By extrapolation, 28/14 ought to exist I think which would be really exciting. But is it just out of reach? – None – 2016-11-05T19:50:19.660

n=14 would take about 200 days with my current code on my 8-core CPU. The code can probably be sped up by 30% or so. After that I run out of ideas. And your extrapolation seems rather optimistic anyways... – Ton Hospel – 2016-11-05T19:56:21.780

I think the circulant 50 by 25 matrix with first row 01011011100010111101000001100111110011010100011010 might work. This was found by an optimization heuristic which might be useful. – None – 2016-11-05T20:10:17.990

140 hours for exhaustive cover of n=14 is very impressively fast I have to say. – None – 2016-11-06T07:04:46.667

This is pretty amazing. How did you do it? – None – 2016-11-19T12:50:41.440

2

Haskell 14/8 = 1.75

1 1 0 0 0 1 0 1 1 0 1 1 0 0
1 1 1 0 0 0 1 0 1 1 0 1 1 0
0 1 1 1 0 0 0 1 0 1 1 0 1 1
1 0 1 1 1 0 0 0 1 0 1 1 0 1
0 1 0 1 1 1 0 0 0 1 0 1 1 0
0 0 1 0 1 1 1 0 0 0 1 0 1 1
0 0 0 1 0 1 1 1 0 0 0 1 0 1
0 0 0 0 1 0 1 1 1 0 0 0 1 0

Previously 9/6 = 1.5

1 0 1 0 1 1 0 0 1
1 1 0 1 0 1 1 0 0
1 1 1 0 1 0 1 1 0
1 1 1 1 0 1 0 1 1
1 1 1 1 1 0 1 0 1
1 1 1 1 1 1 0 1 0

I wrote this, then looked at the answers to the other question and was... discouraged.

import Data.List
import Data.Hashable
import Control.Monad
import Control.Parallel.Strategies
import Control.Parallel
import qualified Data.HashSet as S

matrix§indices = [ matrix!!i | i<-indices ]

powerset :: [a] -> [[a]]
powerset = filterM (const [True, False])

hashNub :: (Hashable a, Eq a) => [a] -> [a]
hashNub l = go S.empty l
    where
      go _ []     = []
      go s (x:xs) = if x `S.member` s
        then go s xs
        else x : go (S.insert x s) xs

getMatrix :: Int -> Int -> [Int] -> [[Int]]
getMatrix width height vector = [ vector § [x..x+width-1] | x<-[0..height-1] ]

hasDuplicate :: (Hashable a, Eq a) => [a] -> Bool
hasDuplicate m = go S.empty m
    where
        go _ [] = False
        go s (x:xs) = if x `S.member` s
            then True
            else go (S.insert x s) xs

hasProperty :: [[Int]] -> Bool
hasProperty matrix =
    let
        base = replicate (length (matrix !! 0)) 0::[Int]
    in
        if elem base matrix then
            False
        else
            if hasDuplicate matrix then
                False
            else
                if hasDuplicate (map (foldl (zipWith (+)) base) (powerset matrix)) then
                    False
                else
                    True


pmap = parMap rseq

matricesWithProperty :: Int -> Int -> [[[Int]]]
matricesWithProperty n m =
    let
        base = replicate n 0::[Int]
    in
    filter (hasProperty) $
    map (getMatrix n m) $
    sequence [ [0,1] | x<-[0..n+m-1] ]

firstMatrixWithProperty :: Int -> Int -> [[Int]]
firstMatrixWithProperty n m = head $ matricesWithProperty n m

main = mapM (putStrLn. show) $ map (firstMatrixWithProperty 8) [1..]

alexander-brett

Posted 2015-04-26T10:01:42.697

Reputation: 1 485

Thank you! A Haskell answer is always very welcome. I think the first interesting case is 12/7. Can you get that? – None – 2015-04-28T13:44:03.930

I was running it on a dual-core 2009 laptop, so no :) I'll try again on a faster machine though – alexander-brett – 2015-04-28T14:03:54.513

Very nice. I just added a comment that 21/11 would beat all previous answers. – None – 2015-04-28T19:15:07.147

Could you explain what your code outputs exactly please? – None – 2015-04-29T11:27:50.527

In main, the number (in this case 8) is the height of the matrix, and the program iterates through widths [1..]. For each height/width combination, it prints out an array of columns of a valid matrix. – alexander-brett – 2015-04-30T06:56:55.280

I awarded you the bounty as the best new answer. Congratulations :) – None – 2015-05-05T18:45:57.633

1

C

Here's an answer which works and is significantly more memory-efficient than the Haskell. Unfortunately, my computer is still too slow to achieve better than 14/8 in any reasonable length of time.

Try compiling with gcc -std=c99 -O2 -fopenmp -o matrices.exe matrices.c and running with matrices.exe width height or similar. The output is an integer, the bits of which form the basis for the matrix in question, for instance:

$ matrices.exe 8 14
...
valid i: 1650223

Then, since 1650223 = 0b110010010111000101111, the matrix in question is:

0 0 1 1 1 0 1 0 0 1 0 0 1 1
0 ...
1 ...
0
1
1
1
1

If someone with 8 cores and time on their hands wants to run this for a while, I think some good things could result :)


#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

/*
 * BEGIN WIKIPEDIA CODE
 */
const long long m1  = 0x5555555555555555; //binary: 0101...
const long long m2  = 0x3333333333333333; //binary: 00110011..
const long long m4  = 0x0f0f0f0f0f0f0f0f; //binary:  4 zeros,  4 ones ...
const long long m8  = 0x00ff00ff00ff00ff; //binary:  8 zeros,  8 ones ...
const long long m16 = 0x0000ffff0000ffff; //binary: 16 zeros, 16 ones ...
const long long m32 = 0x00000000ffffffff; //binary: 32 zeros, 32 ones
const long long hff = 0xffffffffffffffff; //binary: all ones
const long long h01 = 0x0101010101010101; //the sum of 256 to the power of 0,1,2,3...
//This uses fewer arithmetic operations than any other known
//implementation on machines with fast multiplication.
//It uses 12 arithmetic operations, one of which is a multiply.
long long hamming(long long x) {
    x -= (x >> 1) & m1;             //put count of each 2 bits into those 2 bits
    x = (x & m2) + ((x >> 2) & m2); //put count of each 4 bits into those 4 bits
    x = (x + (x >> 4)) & m4;        //put count of each 8 bits into those 8 bits
    return (x * h01)>>56;  //returns left 8 bits of x + (x<<8) + (x<<16) + (x<<24) + ...
}
/*
 * END WIKIPEDIA CODE
 */

int main ( int argc, char *argv[] ) {
    int height;
    int width;

    sscanf(argv[1], "%d", &height);
    sscanf(argv[2], "%d", &width);

    #pragma omp parallel for
    for (
        /*
         * We know that there are 2^(h+w-1) T-matrices, defined by the entries
         * in the first row and first column. We'll let the long long i
         * represent these entries, with 1s represented by set bits.
         *
         * The first (0) and last (1) matrix we will ignore.
         */
        long long i = 1;
        i < (1 << (height+width-1))-1;
        i++
    ) {
        // Flag for keeping track as we go along.
        int isvalid = 1;

        /*
         * Start by representing the matrix as an array of columns, with each
         * non-zero matrix entry as a bit. This allows us to construct them and
         * check equality very quickly.
         */
        long *cols = malloc(sizeof(long)*width);
        long colmask = (1 << height)-1;
        for (int j = 0; j < width; j++) {
            cols[j] = (i >> j) & colmask;
            if (cols[j] == 0) {
                //check no zero rows
                isvalid = 0;
            } else {
                //check no duplicate rows
                for (int k = 0; k < j; k++) {
                    if (cols[j] == cols[k]) {
                        isvalid = 0;
                    }
                }
            }
        }

        if (isvalid == 1) {
            /*
             * We'll also represent the matrix as an array of rows, in a
             * similar manner.
             */
            long *rows = malloc(sizeof(long)*height);
            long rowmask = (1 << width)-1;
            for (int j = 0; j < height; j++) {
                rows[j] = (i >> j) & rowmask;
            }

            int *sums[(1 << width)];
            for (long j = 0; j < 1<<width; j++) {
                sums[j] = (int*)malloc(sizeof(int)*height);
            }

            for (
                /*
                 * The powerset of columns has size 2^width. Again with the
                 * long; this time each bit represents whether the
                 * corresponding row is a member of the subset. The nice thing
                 * about this is we can xor the permutation with each row,
                 * then take the hamming number of the resulting number to get
                 * the sum.
                 */
                long permutation = 1;
                (isvalid == 1) && (permutation < (1 << width)-1);
                permutation ++
            ) {
                for (int j = 0; j < height; j++) {
                    sums[permutation][j] = hamming( rows[j] & permutation);
                }
                for (int j = permutation-1; (isvalid == 1) && (j > -1); j--) {
                    if (memcmp(sums[j], sums[permutation], sizeof(int)*height) == 0) {
                        isvalid = 0;
                    }
                }
            }

            for (long j = 0; j < 1<<width; j++) {
                free(sums[j]);
            }

            free(rows);

        }

        if (isvalid == 1) {
            printf ("valid i: %ld\n", i);
        }

        free(cols);
    }

    return 0;
}

alexander-brett

Posted 2015-04-26T10:01:42.697

Reputation: 1 485

I get alexander-brett.c: In function ‘main’: alexander-brett.c:107:21: warning: implicit declaration of function ‘memcmp’ [-Wimplicit-function-declaration] if (memcmp(sums[j], sums[permutation], sizeof(int)*height) == 0) { ^ alexander-brett.c:122:13: warning: format ‘%ld’ expects argument of type ‘long int’, but argument 2 has type ‘long long int’ [-Wformat=] printf ("valid i: %ld\n", i); – None – 2015-05-01T16:50:06.063

How long does ./alexander-brett 8 14 take for you? – None – 2015-05-01T18:24:21.483

Hi Lembik, 8 14 got 5 answers in an hour for me on a quad core machine. I managed to compile with those headers on windows, it'd be wierd if memcmp went missing... – alexander-brett – 2015-05-01T18:26:44.737

I tried your code with 7 12 and one of the answers it outputs is valid i: 7481. In python bin(7481) is 0b1110100111001 which isn't long enough. Any idea what is going on? – None – 2015-05-19T19:36:34.650

Huh, that's pretty upsetting... possible it's supposed to start with a 0, ie 0b01110100111001? – alexander-brett – 2015-05-19T20:24:03.033

It's worse than an out by one problem. Could you maybe try it yourself to see? – None – 2015-05-19T21:23:07.497