"""
Redelmeier meets dynamic programming.


Redelmeier's algorithm with the untried set organized as a queue.
Prototype implementation for enumerating polyominoes up to size nmax.
The queue is implemented as an array.

Combination with DP-idea.
At a certain "level", the state is stored, and paths leading to the
same state are coalesced.

"level" has several possibities:
* the number of binary decisions so far
* the number of POSITIVE decisions (=n=size of S) (*chosen here)

"state" means a triplet of (S union X, U, |S|)
(if |S| is not constant)

Potential improvement:
Points that are not reachable from U can be added to X.

Here it the selection strategy of the next splitting vertex from U CAN make
a difference.

* queue (*here*)
* stack
* some other (heuristic) order, e.g. the point closest to the origin
  (in order to keep the thing "round" and gather more similar shapes)

depthmax=16 13475 13589 3;          1% savings (for dpmax=16, python is faster than pypy)
depthmax=20 211168 216122 9;      2.3% savings
depthmax=24 3251381 3446914 32;   5.5% savings
depthmax=26 12743522 13773039 75; 7.5% savings, real 116m 7,755s

"""

from collections import defaultdict

boundaryconfigurations =  defaultdict(int)

depthmax = 16

nmax = 10
nmax = 6
nmax = 8

nmax = max(nmax, depthmax)

PRINT_SOLUTIONS = False

queuelength = nmax * 5 + 10 # polyomino plus neighbors plus safety buffer
Q = [0] * queuelength
occupied_or_adjacent = defaultdict(bool)
count = defaultdict(int)

for x in range(nmax):
    occupied_or_adjacent[x,0,-1] = occupied_or_adjacent[-x,0,0] = True
for x in range(-nmax+1,nmax):
    for y in range(1,nmax):
        occupied_or_adjacent[x,y,-1] = occupied_or_adjacent[x,-y,0] = True
# lower border and starting cell

def construct(stackbegin, stackend, n, depth=0):
    """Current polyomino has n cells.
    UNTRIED points are stored in Q[stackbegin] ... Q[stackend-1]."""
    #if depth > depthmax:
    #    return
    count[n] += 1
    if depth == depthmax:
        Untried = frozenset(Q[stackbegin:stackend])
        Used = frozenset(pt for pt,v in occupied_or_adjacent.items() if v).difference(Untried)
        boundaryconfigurations[Untried,Used,n] += 1
        return
    if 0 and n>=nmax: # different version: level = number of included cells
        Untried = frozenset(Q[stackbegin:stackend])
        Used = frozenset(pt for pt,v in occupied_or_adjacent.items() if v).difference(Untried)
        boundaryconfigurations[Untried,Used] += 1
        return
    for i in range(stackbegin, stackend):
        if depth+1+i-stackbegin > depthmax:
            break
        #print(f"{n=} {i=} {stackbegin}:{stackend} {Q[stackbegin:stackend]}",)
        x,y,z = Q[i]
        # include the cell (x,y,z):
        new_neighbors = [nbr for nbr in ((x-1,y,z),(x,y-1,z),
                                         (x+1,y,z),(x,y+1,z),
                                         (x,y,z-1),(x,y,z+1))
                         if not occupied_or_adjacent[nbr]]
        for k,nbr in enumerate(new_neighbors):
            occupied_or_adjacent[nbr]=True
            Q[stackend+k] = nbr

        # recursive call:
        construct(i+1, stackend+len(new_neighbors), n+1, depth+1+i-stackbegin)
        
        for nbr in new_neighbors:
            occupied_or_adjacent[nbr]=False # undo the mark; nbr becomes "unseen"
        
Q[0] = (0,0,0) # The starting square (0,0) is put on the queue
construct(0, 1, 0) # start the enumeration
for i,val in sorted(count.items()):
    print(f"{i:2} {val:6}") # results
compressed, total = len(boundaryconfigurations), sum(boundaryconfigurations.values())
multiplicity = max(boundaryconfigurations.values())
print(f"{depthmax=}, {compressed=}, {total=}, saving={100*(1-compressed/total):.2}%, max.{multiplicity=}")
#for k,v in boundaryconfigurations.items():
#    if v>1: print (v,k)