The Tax Historian

Question

Introduction

There is a tax collector that has some trouble managing the taxes of his kingdom: the historical records have burnt down in a great fire.

He wants to find out how many possible pasts there could be in terms of where the current money was inherited from. Luckily, his kingdom is very simple.

The kingdom can be modelled by a 2D boolean matrix, where l represents someone who has inherited money, and O represents someone who has not. For example:

l O l l

O O O l

l O l O

O O O l

(It will always be a rectangle)

In the next generation, the kingdom is smaller (The wolves are strong!).

The next generation would look like this, superimposed on the previous generation (x is a placeholder for a descendent in the next generation)

l O l l
 x x x
O O O l
 x x x
l O l O
 x x x
O O O l

A descendent will look at the ancestors that are directly around them (So the top-left x will see {l, O, O, O}, called an Unaligned rectangular neighbourhood)

If only one ancestor has inherited money, the descendant will inherit money from them. If more than one ancestor has inherited money, they will squabble and the descendant will end up not inheriting money. If no one has inherited money, the descendant will not inherit any money.

(More than one descendant can inherit from one ancestor)

So, the next generation would look like:

​
 l l O

 l l O

 l l O
​

Challenge

Input

The current state of the generation, as an array of arrays of any two distinct values, where the inner arrays are all of the same length.

E.g., for the example above, it could be:

[
  [True, True, False],
  [True, True, False],
  [True, True, False]
]

Output

An integer representing the number of unique previous generations where the next generation is the input.

You can assume that the answer will always be less than 2^30 - 1. (or 1073741823).

The previous generation would be called a "preimage" and this challenge would be to count the preimages.

Scoring

This is a fastest-code challenge, so each submission will be tested on my computer, and the submission that takes the least time will be the winner.

Example Input and Output

(Where 1 is a descendent who inherited money, and 0 is a descendent who did not inherit money)

Input:

[[1, 0, 1],
 [0, 1, 0],
 [1, 0, 1]]

Output:

4

Input:

[[1, 0, 1, 0, 0, 1, 1, 1],
 [1, 0, 1, 0, 0, 0, 1, 0],
 [1, 1, 1, 0, 0, 0, 1, 0],
 [1, 0, 1, 0, 0, 0, 1, 0],
 [1, 0, 1, 0, 0, 1, 1, 1]]

Output:

254

Input:

[[1, 1, 0, 1, 0, 1, 0, 1, 1, 0],
 [1, 1, 0, 0, 0, 0, 1, 1, 1, 0],
 [1, 1, 0, 0, 0, 0, 0, 0, 0, 1],
 [0, 1, 0, 0, 0, 0, 1, 1, 0, 0]]

Output:

11567

Answer 1

C++ using the BuDDy library

This seemed like a nice excuse to play with binary decision diagrams. The kingdom is converted into a big boolean formula where we have to count the number of ways in which it can be satisfied. That can (sometimes) be done more efficiently than it sounds.

The kingdom must be given as a program constant as a flat array and explicitely given dimensions. (Nice input is left as an execise for the reader :-)

Here is the embarrassingly simple code:

#include <iostream>
#include <bdd.h>

// describe the kingdom here:

constexpr int ROWS = 4;
constexpr int COLS = 10;

constexpr int a[] = {
   1, 1, 0, 1, 0, 1, 0, 1, 1, 0,
   1, 1, 0, 0, 0, 0, 1, 1, 1, 0,
   1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
   0, 1, 0, 0, 0, 0, 1, 1, 0, 0,
};

// end of description

// check dimensions
static_assert(ROWS*COLS*sizeof(int)==sizeof(a),
          "ROWS*COLS must be the number of entries of a");

// dimensions of previous generation
constexpr int R1 = ROWS+1;
constexpr int C1 = COLS+1;

// condition that exactly one is true
bdd one(bdd a, bdd b, bdd c, bdd d){
  bdd q = a & !b & !c & !d;
  q |= !a & b & !c & !d;
  q |= !a & !b & c & !d;
  q |= !a & !b & !c & d;
  return q;
}

int main()
{
  bdd_init(1000000, 10000); // tuneable, but not too important
  bdd_setvarnum(R1*C1);
  bdd q { bddtrue };
  for(int j=COLS-1; j>=0; j--) // handle high vars first
    for (int i=ROWS-1; i>=0; i--){
      int x=i+R1*j;
      bdd p=one(bdd_ithvar(x), bdd_ithvar(x+1),
                bdd_ithvar(x+R1), bdd_ithvar(x+R1+1));
      if (!a[COLS*i+j])
        p = !p;
      q &= p;
    }
  std::cout << "There are " << bdd_satcount(q) << " preimages\n";
  bdd_done();
}

To compile with debian 8 (jessie), install libbdd-dev and do g++ -std=c++11 -o hist hist.cpp -lbdd. (Optimizing options make almost no difference because the real work is done in the library.)

Big examples can lead to messages about garbage collection. They could be suppressed, but I prefer to see them.

bdd_satcount returns a double, but that is good enough for the expected range of results. The same counting technique is possible with exact (big) integers.

The code is optimized for ROWS<COLS. If you have a lot more rows than columns, it might be a good idea to transpose the matrix.

Answer 2

Python 2.7

This is just a naïve first attempt. It's not particularly fast, but it is correct.

The first observation is that each cell is dependent on exactly four cells in the previous generation. We can represent those four cells as a 4-bit number (0-15). According to the rules, if exactly one neighboring cell in the previous generation is 1, then a given cell in the current generation will be 1, otherwise, it will be 0. Those correspond to the powers of two, namely, [1, 2, 4, 8]. When the four ancestors are represented as a 4-bit number, any other number will result in a 0 in the current generation. With this information, upon seeing a cell in the current generation, we can narrow down the possibilities of the neighborhood in the previous generation to one of four or one of twelve possibilities respectively.

I've chosen to represent the neighborhood as follows:

32
10

where 0 is the least-significant bit, and so on.

The second observation is that for two adjacent cells in the current generation, the two neighborhoods from the previous generation overlap:

32  32
10  10

or:

32
10

32
10

In the horizontal case, the 2 from the left neighborhood overlaps with the 3 from the right neighborhood, and similarly with the 0 on the left and the 1 on the right. In the vertical case, the 1 from the top neighborhood overlaps with the 3 from the bottom neighborhood, and similarly with the 0 on the top and the 2 on the bottom.

This overlap means that we can narrow down the possibilities for yet-unchosen neighborhoods based on what we've already chosen. The code works it's way from left to right, top to bottom, in a recursive depth-first search for possible preimages. The following diagram indicates which previous neighborhoods we have to consider when looking at the possible neighborhoods of a current cell:

f = free choice
h = only have to look at the neighborhood to the left
v = only have to look at the neighborhood to the top
b = have to look at both left and top neighborhoods

[f, h, h, h, h],
[v, b, b, b, b],
[v, b, b, b, b],
[v, b, b, b, b]

Here's the code:

def good_horizontal(left, right):
    if (left & 4) >> 2 != (right & 8) >> 3:
        return False
    if left & 1 != (right & 2) >> 1:
        return False
    return True


def good_vertical(bottom, top):
    if (bottom & 8) >> 3 != (top & 2) >> 1:
        return False
    if (bottom & 4) >> 2 != (top & 1):
        return False
    return True


ones = [1, 2, 4, 8]
zeros = [0, 3, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15]
h = {}
v = {}

for i in range(16):
    h[i] = [j for j in range(16) if good_horizontal(i, j)]
    v[i] = [j for j in range(16) if good_vertical(i, j)]


def solve(arr):
    height = len(arr)
    width = len(arr[0])

    if height == 1 and width == 1:
        if arr[0][0] == 1:
            return 4
        else:
            return 12
    return solve_helper(arr)


def solve_helper(arr, i=0, j=0, partial=None):
    height = len(arr)
    width = len(arr[0])

    if arr[i][j] == 1:
        poss = ones
    else:
        poss = zeros

    if i == height - 1 and j == width - 1:  # We made it to the end of this chain
        if height == 1:
            return sum([1 for p in poss if p in h[partial[-1][-1]]])
        else:
            return sum([1 for p in poss if partial[-2][-1] in v[p] and p in h[partial[-1][-1]]])

    if j == width - 1:
        new_i, new_j = i + 1, 0
    else:
        new_i, new_j = i, j + 1

    if i == 0:
        if j == 0:
            # first call
            return sum([solve_helper(arr, new_i, new_j, [[p]]) for p in poss])
        # still in the first row
        return sum([solve_helper(arr, new_i, new_j, [partial[0] + [p]]) for p in poss if p in h[partial[0][-1]]])
    if j == 0:  # starting a new row
        return sum([solve_helper(arr, new_i, new_j, [r for r in partial + [[p]]]) for p in poss if partial[i - 1][0] in v[p]])
    return sum([solve_helper(arr, new_i, new_j, [r for r in partial[:-1] + ([partial[-1] + [p]])]) for p in poss if p in h[partial[i][-1]] and partial[i - 1][j] in v[p]])

To run it:

test3 = [
    [1, 1, 0, 1, 0, 1, 0, 1, 1, 0],
    [1, 1, 0, 0, 0, 0, 1, 1, 1, 0],
    [1, 1, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 1, 0, 0, 0, 0, 1, 1, 0, 0]
]

expected3 = 11567

assert(solve(test3) == expected3)