# Fastest count of certain hypercubes with labeled vertices

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.

• TIO is not great for timing. You might have to commit to running the code on your PC for a fixed amount of time.
– Simd
Commented Feb 6 at 17:37
• You cannot use both fastest-code and fastest-algorithm. According to your description, the winning criterion is fastest-code. Commented Feb 6 at 17:42
• @Arnauld I have removed fastest-algorithm. Commented Feb 6 at 17:47
• @Simd I have modified the SCORE section and I hope it is acceptable now. Commented Feb 6 at 17:48
• Surely a rust coder wants to show off their language?
– Simd
Commented Feb 14 at 19:48

# 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()
return count

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

• Thank you. Do you think this can be executed on tio.run or other online platform? If not, what do you advise to install on my Windows PC? Commented Feb 7 at 16:42
• @FabiusWiesner install python 3, and run 'pip install python-sat' in the terminal. Or I guess you could run it in a notebook in Google Collab if that's too much hassle. Commented Feb 7 at 16:43
• This is very fast for n = 1..4. n = 5 isn't too slow either.
– Simd
Commented Feb 7 at 17:01
• What is the point of the 39 iterations where nothing is outputted?
– Simd
Commented Feb 7 at 17:04
• @Simd the 40 iterations are due to my request, it is just to have more resolution on the measured time. This code is very good and takes 23.5 seconds to compute (one time) from 1 to 5, however I prefer to leave the challenge as it is, there is enough resolution to discriminate other solutions. Commented Feb 7 at 17:35

# 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);
}

• Wonderful, thank you! This did the expected value 15028134 for n=6 in only 2 hours and 45 minutes. However, I fear that a more sophisticated algorithm is needed to try to compute the case n=7. Commented Feb 15 at 16:22
• I'm glad to see this in Rust. My answer was written with the vague intention to port to C. @FabiusWiesner I have an idea for n=7 but I don't know when / if I'll get round to it. Commented Feb 15 at 23:17
• @Level River St when I have time, I want to try to use #SAT solvers for the case n=7. Commented Feb 16 at 7:23

# 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

return ((1<<(1<<n))-1)//((1<<(1<<k))+1)

def swap_parts(n, k1, 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):
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:
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
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:
continue
good.append(x)
compatible = {}
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):
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
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:
continue
good.append(x)
compatible = {}
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):
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)

• This is very impressive!
– Simd
Commented Feb 16 at 11:32
• Thank you! I will try to do n=7 with this one, or at least evaluate how long it could take. Commented Feb 16 at 12:46
• @FabiusWiesner The updated version got the expected answer for n=7 (8378070864) in about an hour. (Numpy is much faster, and it also doesn't try to store the output list for the final hypercube, which for n=7 would normally require like 128 gigabytes of memory). Commented Feb 16 at 14:55
• +1 this is exactly what I was thinking of doing next. But it is ... long! Looks like a lot of work. Commented Feb 17 at 20:18

# 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!

• Thank you. Maybe it could be furtherly optimized by checking that the two (n-1)-faces (1882 x 1882) are valid on all the n axes before checking the rectangles. Commented Feb 12 at 8:50
• @FabiusWiesner thanks. I'm not 100% sure what you mean, because you need to do 3 checks per rectangle in order to be sure it's bad. However I did try changing the order in which the 3 checks for each rectangle are made to see if it makes a difference. I first carry out two checks that span between the 2 n-1 cubes (which should each have a 50% chance of showing a rectangle is OK), and then finally do a check within the same n-1 cube (which should have a greater than 50% chance of being OK.) Moving this final check to the beginning of the three checks didn't make the code noticeably faster. Commented Feb 15 at 23:31
• The other thing I ought to do is get rid of that product because it uses too much memory to run n=6 Commented Feb 15 at 23:32
• Basically I was thinking about the fact that what is done in the new dimension for 1882 x 1882 could be checked in all n dimensions, although I didn't evaluate its complexity, just an idea. Commented Feb 16 at 7:22