target = 26
CMAX = (target-1)//2
# for odd target, one could reduce CMAX by one and
# compute the missing ones by a formula.
# target 33, CMAX = 16, runs in 6 min on my home computer with 8 GiB.
# target 31, CMAX = 15, runs in 2 min. (factor 3)
# n=50 might take a month, but the space is more crucial.
# More aggressive pruning of states might be required.
"""https://oeis.org/A346800
Number of fixed polyominoes with n cells that have a diagonal axis of
symmetry going from lower left to upper right.
"""
A346800 = [1, 0, 2, 1, 5, 4, 16, 13, 54, 46, 186, 167, 660, 612, 2384,
2267, 8726, 8464, 32278, 31822, 120419, 120338, 452420, 457320,
1709845, 1745438, 6494848, 6686929, 24779026, 25703792, 94899470,
99096382, 364680344] # more values are known
"""https://oeis.org/A351127
Number of polyominoes with n cells with maximal symmetry and center of
rotation about the center of a square.
"""
A351127 = [1, 0, 0, 0, 1, 0, 0, 1, 2, 0, 0, 1, 2, 0, 0, 4, 4, 0, 0, 7,
7, 0, 0, 16, 11, 0, 0, 29, 20, 0, 0, 67, 36, 0, 0, 119, 65, 0, 0, 264,
117, 0, 0, 478, 216, 0, 0, 1043, 396, 0, 0, 1910, 736, 0, 0, 4116,
1369, 0, 0, 7592, 2558, 0, 0, 16201, 4787, 0, 0, 30114, 9006]
"""
num_poly[c] counts partial polyominoes in the first c
columns, using the cells 0<=j<=i1:
target = int(sys.argv[1])
CMAX = (target-1)//2
def succ0(L,i):
"""L = L0,L1,...,L[i-1],L[i],...,L[n-1]
The kink is between L[i-1] and L[i].
L_i is replaced by 0 if this does not isolate L[i]'s component.
"""
if L[i]==0:
return L
left_part,prev,right_part = L[:i], L[i], L[i+1:]
if any(x==prev for x in left_part+right_part):
return rename_by_appearance(left_part+(0,)+right_part)
else: # component prev is isolated
return None
def succ1(L,i):
"""L = (L0,L1,...,L[i-1],L[i],...,L[n-1])
The new L_i is occupied.
"""
prev = L[i]
if i>0:
left = L[i-1]
else:
left = 0
if prev:
if left in (0,prev):
return L
else:
# two components are merged
old = max(left,prev)
new = min(left,prev)
# rename all appearances of old to new
return tuple(x if x=1 cells before L, of which only the first is occupied"
if r==1 and L[0]>0: return (L[0],) + L # L[0]==1
return (1,) + (r-1)*(0,) + tuple(x+1 if x>0 else 0 for x in L)
def add_after(L,r):
"add r>=1 cells at the end of L, of which only the last is occupied"
if r==1:
return succ1_end(L)
return L + (r-1)*(0,) + (max(L)+1,)
num_poly = [0,0]
num_poly[0] = { (1,): { (1,0) : 1 } # the singleton polyomino
} # column 0
num_poly[1] = { (1,1): { (1,1) : 1, # the vertical domino
(2,1) : 1, # the triomino
},
(1,0): { (1,1) : 1, # the horizontal domino}
} } # column 1
count_poly = defaultdict(int)
count_poly_w = defaultdict(int); count_poly_w[2]=0
def accumulate_results(counts):
# collect the results
length = 0
for state,nums in counts.items():
length += len(nums)
if max(state)==1: # connected
for (B,I),num in nums.items():
count_poly[B,I] += num
count_poly_w[B+2*I] += num
# B+2*I cells when reflected along the diagonal
#if B+2*I==13: print("UUU",state,B,I,num)
return length
for counts in num_poly:
accumulate_results(counts) # take into account the initial values
work = space = 0
def update(newstate, map_parameters, subtract=False):
# auxiliary routine
"map_parameters maps the original parameters (B,I) to new ones."
# nums and newcounts are global variables!
global work
nc = newcounts[newstate]
for (B,I),num in nums.items():
BB,II = map_parameters(B,I)
if BB+2*II<=target:
if subtract:
nc[BB,II] -= num
else:
nc[BB,II] += num
work += len(nums)
for c in range(2, CMAX+1):
"""
add column c cell by cell from the bottom.
Column c has c+1 cells.
Altogether, cells are added in this order.
.
10 .
6 9 .
3 5 8 .
1 2 4 7 .
"""
oldcounts = num_poly[c-1]
for r in range(c): # add cell (c,r)
newcounts = defaultdict(lambda : defaultdict(int))
for state,nums in oldcounts.items():
newstate = succ1(state,r)
update(newstate, lambda B,I : (B,I+1))
newstate = succ0(state,r)
if newstate:
update(newstate, lambda B,I : (B,I)) # identity
oldcounts = newcounts
# now add the diagonal cell (c,c)
newcounts = defaultdict(lambda : defaultdict(int))
for state,nums in oldcounts.items():
newstate = succ1_end(state)
update(newstate, lambda B,I : (B+1,I))
newstate = succ0_end(state)
update(newstate, lambda B,I : (B,I))
# now consider those polyominoes that touch the bottom and/or
# top border only in the last column.
# ADD those where the top cell (c,c) is the only occupied
# cell on the diagonal.
# ADD those where the bottom cell (c,0) is the only occupied
# cell on the lowest row.
# SUBTRACT those where both conditions hold.
#
# We count these classes by removing the respective top and bottom
# boundary cells.
# The resulting polyomino lives in a smaller translated wedge of size c0.
for c0 in range(c):
#print(c0,"IN")
#for x,v in sorted(newcounts.items()): print(x,c0,dict(v))
np = num_poly[c0] # these have c0+1 columns
# c-c0 cells need to be added.
# add all additional cells above c0 with a singleton occupied cell
for state,nums in np.items():
newstate = add_after(state,c-c0)
update(newstate, lambda B,I : (1, B+I))
# add all additional cells below with a singleton occupied cell
for state,nums in np.items():
newstate = add_before(state,c-c0)
update(newstate, lambda B,I : (B, I+1))
# both above and below:
for top in range(1,c-c0):
bottom = c-c0-top
# add "top" cells at the top and "bottom" cells at the bottom
#print(f"{c=} {c0=} {top=} {bottom=}")
for state,nums in np.items():
newstate = add_after(add_before(state, bottom), top)
update(newstate, lambda B,I : (1, B+I+1), subtract=True)
newcounts[(1,)+(c-1)*(0,)+(2,)][1,1] -= 1 # extra correction:
# two isolated cells at the extreme end and nothing else.
space += accumulate_results(newcounts) # collect and accumulate the results
num_poly.append(dict(newcounts))
# could convert dicts to lists to save space?
#print(c,"::::")
#for x,v in sorted(newcounts.items()): print(x,dict(v))
print("****",f"{c=} {CMAX=}",list(sorted(count_poly_w.items())))
print(end=f"work {work}, space {space}. ")
memstat(extensive=True)