Fastest count of certain hypercubes with labeled vertices

Question

CHALLENGE

This is a fastest-code challenge.

Count how many n-dimensional hypercubes with n=1,2,3,4 exist, with vertices labeled with either 1 or 0, such that there does not exist any rectangle formed by 4 vertices such that one diagonal is 1-1 and the other diagonal is 0-0. To be more clear, the forbidden rectangles have the form:

10 or 01
01    10

Note that we are talking about all rectangles, not just the faces of the hypercube. For example in the cube:

   0*---1*
  /|   /|     
 0----0 |
 | 0--|-0
 |/   |/
 1*---0*

there is an invalid rectangle formed by the four vertices marked with a star (*), although all faces are valid.

The list of all known expected results is at OEIS A000609. This is actually a conjecture, but a probable one. For more information and some "hints", see the BACKGROUND section.

The count for n=1,2,3,4 must be repeated in a loop 40 times (just to have more resolution on the measured time), and the results for n=1,2,3,4 (4,14,104,1882 respectively) printed on the first round of the loop only.

You should provide a link for tio.run.

Here is my best in C (about 30 s).

BACKGROUND (optional)

I have asked a question related to this challenge at mathoverflow. See the links within that question and especially the one related to the generation of the rectangles.

Another possible idea is using the inclusion-exclusion principle: see this question.

SCORE

I will make a first rough placement list based on the time shown in the tio.run link provided, in the "Debug" section, "User time". Some time later a more precise one with times tested on my machine, since performance on tio.run may vary over time.

RECORDS

There will be another record list if someone is able to compute the case n = 5 or more in a reasonable time (let's say less than 10 minutes for n = 5 and 3 days for n = 6). In this case it may be that I won't double check, so I will trust your claim. I personally computed the case n = 5 in about 5 minutes with a slightly modified version of the above linked C software (which took 8 minutes).

STANDINGS

FASTEST CODE	USER	LANGUAGE	BEST TIME ON MY MACHINE
1	138 Aspen	Rust	0.1364747 s
2	gsitcia	Python	0.1462836 s
3	user1502040	Python	1.0249261 s
4	Level River St	Ruby	1.3259241 s
5	Me	C	25.888948 s

"n RECORD"	USER	LANGUAGE	RECORD and TIME MY PC
1	gsitcia	Python	n = 7 in 1h 19m 26s
2	gsitcia	Python	n = 6 in 1.946 s
3	138 Aspen	Rust	n = 6 in 2h 45m 24s
4	user1502040	Python	n = 6 in 22h 6m
5	gsitcia	Python	n = 5 in 0.156 s
6	138 Aspen	Rust	n = 5 in 1.947 s
7	user1502040	Python	n = 5 in 25.372 s
8	Level River St	Ruby	n = 5 in 27.900 s
9	Me	C	n = 5 in ~480 s

The time for "n RECORD" is for a single computation of the specified value for n.

Answer 1

Python + python-sat (about 1 second on my machine)

import itertools
from pysat.solvers import Glucose4

def count_cubes(n):
    clauses = []
    for (x, y, z) in itertools.combinations(range(1 << n), 3):
        e0_p = y & (~x)
        e0_n = x & (~y)
        e0 = e0_p | e0_n
        e1_p = z & (~x)
        e1_n = x & (~z)
        e1 = e1_p | e1_n
        if e0 & e1:
            continue
        w = x ^ (e0 | e1)
        x, y, z, w = x + 1, y + 1, z + 1, w + 1
        clauses.append([x, w, -y, -z])
        clauses.append([-x, -w, y, z])
    if not clauses:
        assert n == 1
        clauses.append([-1, 1, -2, 2])
    with Glucose4(bootstrap_with=clauses, incr=True) as m:
        count = 0
        while m.solve():
            count += 1
            model = m.get_model()
            m.add_clause([-v for v in model])
    return count

for t in range(40):
    for n in range(1, 5):
        r = count_cubes(n)
        if t == 0:
            print(n, r)

Answer 2

Rust, n=4 repeated 40 times in 0.6 seconds on TIO

A port of @Level River St's Ruby answer in Rust.

---

Try it online!

const N: usize = 4;

fn main() {
    let mut print_output = true;
    let mut a: Vec<u64> = vec![0];
    rec_function(1, &mut a, &mut print_output);
    
    for _ in 0..40 {
        print_output = false;
        rec_function(1, &mut a, &mut print_output);
    }
}

fn rec_function(n: usize, a: &mut Vec<u64>, print_output: &mut bool) {
    if n > N {
        return;
    }

    let mut u: Vec<(u64, u64, u64, u64)> = vec![];

    for i in 0..6_u64.pow((n - 1) as u32) {
        let mut j = i * 6;
        let (mut p, mut q, mut r, mut s) = (0, 0, 0, 0);
        let mut t = 0;
        for _ in 0..n {
            let k = j % 6;
            p = p << 1 | (k & 1);
            q = q << 1 | (41 >> k & 1);
            r = r << 1 | (25 >> k & 1);
            s = s << 1 | (26 >> k & 1);
            if k > 3 {
                t = k;
            }
            j /= 6;
        }
        if t > 4 {
            u.push((p, q, r, s));
        }
    }

    let mut b: Vec<u64> = vec![];
    let pairs = [0, (1 << (1 << (n - 1))) - 1];

    for &c in pairs.iter() {
        for &d in a.iter() {
            for &e in a.iter() {
                let e = (e ^ c) | d << (1 << (n - 1));
                let mut h = false;
                for &(p, q, r, s) in &u {
                    let g = e >> p & 1;
                    h = g != (e >> q & 1) && g == (e >> r & 1) && g != (e >> s & 1);
                    if h {
                        break;
                    }
                }
                if !h {
                    b.push(e);
                }
            }
        }
    }

    if *print_output {
        println!("{}", b.len() * 2);
    }
    rec_function(n + 1, &mut b, print_output);
}

Answer 3

Python 3 (+Numpy) (n=4 repeated 40 times in 0.15 seconds, n=6 in 2.5 seconds on TIO)

import itertools
import functools
from collections import Counter
import time
import numpy as np

@functools.lru_cache(None)
def gen_rectangles(n):
    x = []
    for t in itertools.product([
        (0,0,0),
        (1,1,1),
        (0,0,1),
        (1,1,0),
        (0,1,0),
        (1,0,1),
    ], repeat=n):
        a = b = c = 0
        for a1, b1, c1 in t:
            a, b, c = a*2+a1, b*2+b1, c*2+c1
        if a < b < c:
            x.append((a, b, c, a^b^c))
    return x

def mask(n, k):
    return ((1<<(1<<n))-1)//((1<<(1<<k))+1)

def swap_parts(n, k1, k2):
    m1 = mask(n, k1)
    m2 = mask(n, k2)
    m12 = m1 ^ m2
    ma =  m1 & ~m2
    mb = ~m1 &  m2
    d = 2**k2 - 2**k1
    return m12, ma, mb, d

def swap(x, m12, ma, mb, d):
    xs = x & m12
    a = x & ma
    b = x & mb
    return x ^ xs ^ (a>>d) ^ (b<<d)

def flip_parts(n, k):
    m = mask(n, k)
    k1 = 1 << k
    return m, k1

def flip(x, m, k1):
    a = x & m
    b = (x >> k1) & m
    return (a << k1) | b

@functools.lru_cache(None)
def parts(n):
    return [flip_parts(n, k) for k in range(n)], [swap_parts(n, k1, k2) for k2 in range(n) for k1 in range(k2)]

def reduce(x, n):
    fp, sp = parts(n)
    for p in fp:
        x1 = flip(x, *p)
        if x1 < x:
            return x1, 0, p
    for p in sp:
        x1 = swap(x, *p)
        if x1 < x:
            return x1, 1, p

def make_order(d):
    d = {a:b[0] for a,b in d.items()}
    C = Counter(d.values())
    empty = d.keys() - d.values()
    l = []
    while empty:
        i = empty.pop()
        j = d.get(i, None)
        if j is None: continue
        l.append(i)
        k = C[j] - 1
        if k == 0:
            empty.add(j)
        else:
            C[j] = k
    return l

@functools.lru_cache(None)
def gen_cubes(n):
    if n == 0:
        return [0, 1]
    if n == 1:
        return [0, 1, 2, 3]
    prev = gen_cubes(n-1)
    half = 2 ** (n-2)
    k = 2 ** (n-1)
    full = 2 ** k - 1
    n1 = n - 1
    bad = {}
    good = []
    for x in prev:
        bc = bin(x).count('1')
        if bc > half: 
            bad[x] = full ^ x, 2, full
            continue
        t = reduce(x, n1)
        if t:
            bad[x] = t
            continue
        good.append(x)
    order = make_order(bad)
    compatible = {}
    masks = []
    for a,b,c,d in gen_rectangles(n):
        if b >= k or c < k: continue
        a, b, c, d = 1 << a, 1 << b, 1 << (c-k), 1 << (d-k)
        masks.append((a|b, c|d, a, d, b, c))
    num_good = sum(1 for b in good for t in masks if b & t[1] in (t[3], t[5]))
    prev = np.array(prev, dtype=np.uint64)
    for b in good:
        l = prev
        for ma, mb, pa1, pb1, pa2, pb2 in masks:
            q = b & mb
            if q == pb1:
                pa = pa1
            elif q == pb2:
                pa = pa2
            else:
                continue
            l = l[l&ma != pa]#[a for a in l if a&ma != pa]
        compatible[b] = list(map(int, l))
    for x in reversed(order):
        pre, disc, p = bad[x]
        l = compatible[pre]
        if disc == 0:
            l = [flip(i, *p) for i in l]
        elif disc == 1:
            l = [swap(i, *p) for i in l]
        else:
            l = [p ^ i for i in l]
        compatible[x] = l
    out = []
    for x, l in compatible.items():
        x <<= k
        out.extend(x|i for i in l)
    return out

def count_cubes(n):
    if n == 0:
        return 2
    if n == 1:
        return 4
    prev = gen_cubes(n-1)
    half = 2 ** (n-2)
    k = 2 ** (n-1)
    full = 2 ** k - 1
    n1 = n - 1
    bad = {}
    good = []
    for x in prev:
        bc = bin(x).count('1')
        if bc > half: 
            bad[x] = full ^ x, 2, full
            continue
        t = reduce(x, n1)
        if t:
            bad[x] = t
            continue
        good.append(x)
    order = make_order(bad)
    compatible = {}
    masks = []
    for a,b,c,d in gen_rectangles(n):
        if b >= k or c < k: continue
        a, b, c, d = 1 << a, 1 << b, 1 << (c-k), 1 << (d-k)
        masks.append((a|b, c|d, a, d, b, c))
    num_good = sum(1 for b in good for t in masks if b & t[1] in (t[3], t[5]))
    prev = np.array(prev, dtype=np.uint64)
    for b in good:
        l = prev
        for ma, mb, pa1, pb1, pa2, pb2 in masks:
            q = b & mb
            if q == pb1:
                pa = pa1
            elif q == pb2:
                pa = pa2
            else:
                continue
            l = l[l&ma != pa]#[a for a in l if a&ma != pa]
        compatible[b] = len(l)
    for x in reversed(order):
        pre, disc, p = bad[x]
        compatible[x] = compatible[pre]
    return sum(compatible.values())

if __name__ == '__main__':
    t = time.time()
    for i in range(7):
        print(i, count_cubes(i))
    print(f'{time.time() - t:.1f} seconds 0-6')
    t = time.time()
    for _ in range(40):
        gen_cubes.cache_clear()
        parts.cache_clear()
        gen_rectangles.cache_clear()
        for i in range(5):
            gen_cubes(i)
    print(f'{time.time() - t:.3f} seconds 0-4 x40')

Try it online!

Based on @Level River St's Ruby answer. The main difference is that this answer exploits additional symmetry.

As a simple example, consider a single square of the hypercube. There are 3 cases:

• All four vertices are colored the same (2 cases like this (flip))
• A single vertex is colored differently (8 cases like this (flip/rotate))
• 2 adjacent vertices are colored 1 (4 cases like this (rotate only))

This answer considers the symmetries of the n-1 hypercube, drastically reducing the amount of work necessary (in the case of n=6, it reduces 94572 cases to only 63 cases)

Answer 4

Ruby, n=5 in 25 seconds on TIO

This code has been updated to enable n=6 to be reached without memory issues.

Specifically [0,(1<<2**n/2)-1].product(a,a).each has been modified to two nested loops. Facility for counting the iterations of the outer loop #p y has been added but commented out. The number of iterations is equal to the number of solutions for the previous n as explained below.

Timing facilities have also been added. On my machine this program uses about 33% of my 6 core processor and for n=6 it handles about 1000 iterations of the outer loop per hour, which means it would take 4 days to complete 94572 iterations for n=6. I'm glad to see others have used the ideas in this answer. There are some minor optimisations that could be made, but a different approach such as taking advantage of symmetry (see gsitcia's answer) is needed for major progress.

N=5
starttime=Time.now
f=->n=1,a=[0]{return if n>N
#generate rectangles
u=[]
(6**(n-1)).times{|i|
  j=i*6
  p=q=r=s=0
  t=0
  n.times{
    k=j%6
    p=p<<1|(k&1); q=q<<1|(41>>k&1); r=r<<1|(25>>k&1); s=s<<1|(26>>k&1)
    t=k if k>3 
    j/=6
  }
  u<<[p,q,r,s] if t>4
}
#check colourings
b=[]
[0,(1<<2**n/2)-1].product(a).each_with_index{|z,y|c,d=z;#p y
  a.each{|e|
    e=(e^c)|d<<2**n/2
    h=false
    u.each{|w|p,q,r,s=w 
      g=e>>p&1
      h = g!=e>>q&1 && g==e>>r&1 && g!=e>>s&1 
      break if h
    }
  b<<e unless h
  }
} 
p [n,b.size*2]
f[n+1,b]
}
f[]
p [starttime, Time.now, Time.now-starttime]

Try it online!

Ruby, n=4 repeated 40 times in 1.3 seconds on TIO

N=4
printoutput=true
f=->n=1,a=[0]{return if n>N
#generate rectangles
u=[]
(6**(n-1)).times{|i|
  j=i*6
  p=q=r=s=0
  t=0
  n.times{
    k=j%6
    p=p<<1|(k&1); q=q<<1|(41>>k&1); r=r<<1|(25>>k&1); s=s<<1|(26>>k&1)
    t=k if k>3 
    j/=6
  }
  u<<[p,q,r,s] if t>4
}
#check colourings
b=[]
[0,(1<<2**n/2)-1].product(a,a).each{|z|c,d,e=z
  e=(e^c)|d<<2**n/2
  h=false
  u.each{|w|p,q,r,s=w 
    g=e>>p&1
    h = g!=e>>q&1 && g==e>>r&1 && g!=e>>s&1 
    break if h
  }
  b<<e unless h
} 
p b.size*2 if printoutput
f[n+1,b]
}
40.times{f[];printoutput=false}

Try it online!

Explanation

This is a recursive function which takes the solution for the previous n, makes a cartesian product of two copies, calculates all rectangles that bridge between them, then checks the colourings of these bridging rectangles. Because the solution for the previous n is known to be free of illegal colourings, it is only necessary to check rectangles that span the new dimension. This approach means there are only 1882 x 1882 = 3541924 colourings to check for n=5 instead of 2**32 = 4 billion. Also we only calculate colourings where the most significant bit is 0 which halves the number of calculations; the remaining colourings are the same but with all the 0's and 1's flipped. Precalculating the rectangles before checking colourings also helps.

Generating rectangles

Paraphrasing this post referenced in the question, each rectangle can be represented as 4 points [p, q=p+x, r=p+x+y, s=p+y] where x and y are n-dimensional vectors with values of -1, 0 or +1 in each dimension. For the shape to be a rectangle x and y must be orthogonal, which means that for every dimension i, the product x[i]*y[i] must be zero. p[i] can be 0 or +1 but if either x[i] or y[i] is nonzero p[i] is constrained to only one possibility. This gives the following possibilities:

 k    =     0   1   2   3   4   5
 p[i] =     0   1   0   1   0   1
 x[i] =    +1  -1   0   0   0   0
 y[i] =     0   0   0   0  +1  -1

The value of j is treated as a base 6 number with n digits, the last of which is zero. Because the least significant digit of j generates the most significant bit of p,q,r and s (representing the new dimension n) this ensures that we only generate rectangles where x[n] is +1 ensuring that the generated rectangles span the new dimension.

To avoid double counting, we keep track of the last value of k in t where k>3 (hence y[i]<>0) and only accept rectangles where t==5 meaning the last nonzero value of y[i] was negative (not positive.) We also reject cases where y[i] is zero for all dimensions, since if y is zero, the shape is a line, not a rectangle.

Commented code

This is the same as the code for n=5 in 25 second above. The comments do not have significant impact on runtime.

# Set a constant N for recursion depth, and record  the start time
N=5
starttime=Time.now

# Define a recursive function f where n is the number of dimensons
# a is the set of solutions for the previous recursion. Each solution is represented by an integer
# Initially this is [0,1] for n=0 but we consider only colourings where the last bit is zero.

f=->n=1,a=[0]{return if n>N

# generate rectangles
u=[]                         # Make an array to store rectangles as arrays of 4 integers corresponding to indices of vertices
(6**(n-1)).times{|i|         # Iterate through all base 6 numbers with n-1 digits
  j=i*6                      # Multiply 6 by so that the last digit of j is zero. This ensures we only generate rectangles where x=+1 in the new dimension 
  p=q=r=s=0                  # Set up variabes for the 4 vertices
  t=0                        # Set up variable to record a code for value of y
  n.times{
    k=j%6                    # Extract the last digit of j and decode the values for the last bit of p,q,r,s shifting any existing bits left
    p=p<<1|(k&1)
    q=q<<1|(41>>k&1)
    r=r<<1|(25>>k&1)
    s=s<<1|(26>>k&1)
    t=k if k>3               # If y<>0, record the value of k in t: 4 for positive, 5 for negative
    j/=6                     # Shift j one base 6 digit right
  }
  u<<[p,q,r,s] if t>4        # If y is nonzero in at least one dimension, and the last nonzero dimension was negative, add the rectangle to the list
}

#check colourings
b=[]                         # Make an array to store solutions for dimension n
[0,(1<<2**n/2)-1].product(a).# Consider all possible solutions for n-1 from a with 0 in the most significant bit. Consider also an XOR mask for flipping bits
  each_with_index{|z,y|      # Iterate through all possibilities. y gives the index for use in monitoring progress
  c,d=z                      # c is the XOR mask for bit flipping. d is the n-1 cube for the high order bits 
  #p y                       # (Uncomment to display iteration progress)
  a.each{|e|                 # Iterate through all possible solutions for n-1 with 0 in the most significant bit. e is the n-1 cube for the low order bits.
    e=(e^c)|d<<2**n/2        # Proposed solution is e | d<<2**n/2 . c is either 0 or (1<<2**n/2)-1 which is used to flip the low order bits.
    h=false                  # h=false to indicate no illegal rectangles found
    u.each{|w|p,q,r,s=w      # Iterate through rectangles where w is an array of integers [p,q,r,s] corresponding to indices of vertices
      g=e>>p&1               # Find the value g representing the colour (0 or 1) of vertex p   
      h = g!=e>>q&1 &&       # If g is the same as the diagonally opposite vertex r and different from the adjacent vertices p and s...
       g==e>>r&1 && g!=e>>s&1 
      break if h             # ...h becomes true showing the rectangle is illegal. Break.
    }
  b<<e unless h              # If no illegal rectangles were found add e to the list of solutions
  }
} 
p [n,b.size*2]               # Print the value of n and the number of solutions found 
f[n+1,b]                     # (double the size of b since b only contains solutions where the last bit is zero)
}
f[]                          # Call the function, then print the start time, end time and elapsed time.
p [starttime, Time.now, Time.now-starttime]

Try it online!