target = 33 #26
CMAX = (target-1)//2
# for odd target, one could reduce CMAX by one and
# compute the missing ones by a formula:
# 2^n + sum_{i=0}^{n-2} (n-1-i)2^i
# target 31, CMAX = 15, runs in 3 min on my home computer with 8 GiB.
# target 33, CMAX = 16, runs in 10 min, takes 6GiB memory. (factor 3)
# n=50 might take a month, but the space is more crucial.
# More aggressive pruning of states might be required.
"""
**** c=16 CMAX=16
work 426724899, space 11561867. memory use: 5.983 GiB, cpu time: 684.21 sec
total gaps = 102198032, # intervals=196393813 avg. gap 0.52
"""
"""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]
CENTER_IN_SQUARE = False # D8 symmetries with center on a cell corner,A351127
#CENTER_IN_SQUARE = True # D8 symmetries with center on a cell center, A346800
CENTER_SQUARE_OCCUPIED = True
CENTER_SQUARE_OCCUPIED = False
# Two cases that can be treated separately by running
# with different starting conditions.
# The central square is however not counted in n.
# All other cells are counted with 1/4 of their actual number.
# The values that are stored for n are actually valid for 4n or 4n+1.
# If CENTER_SQUARE_OCCUPIED, the two flags are redundant.
if CENTER_IN_SQUARE:
if CENTER_SQUARE_OCCUPIED:
explanation = "Center (0,0) is occupied."
else:
explanation = "Center (0,0) is free."
else:
explanation = "Center is at a corner."
"""
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,)
count_poly_w = defaultdict(int)
def accumulate_results(counts):
# collect the results
length = 0
for (state,touch_bottom,touch_top),nums in counts.items():
length += len(nums)
if max(state)==1 and touch_bottom and touch_top: # connected
for n,num in nums.items():
count_poly_w[n] += num
#if B+2*I==13: print("UUU",state,B,I,num)
return length
sum_gap = num_intervals = work = 0
def update(newstate, parameter_shift=0):
# auxiliary routine
# nums and newcounts are global variables!
global sum_gap, num_intervals, work
# check nums:
if nums:
interval_gap = max(nums) - min(nums) - len(nums) + 1
sum_gap += interval_gap
#if interval_gap>1:
# print(state, nums.keys())
# very often, the keys contain only even numbers
work += len(nums)
num_intervals += 1
nc = None
for n,num in nums.items():
n += parameter_shift
if n<=target:
if nc is None:
nc = newcounts[newstate] # avoid empty dicts
nc[n] += num
if CENTER_IN_SQUARE:
if CENTER_SQUARE_OCCUPIED:
# The central square is however never counted.
newcounts = { ((1,),True,True): { 0 : 1 } } # the singleton polyomino
else:
newcounts = { ((0,),False,False): { 0 : 1 } } # the empty polyomino
else:
newcounts = { ((1,),True,True): { 1 : 1 }, # the singleton polyomino
((0,),False,False): { 0 : 1 } # the empty polyomino
} # column 0
count_poly_w[1]=1
count_poly_w[2]=0
print(explanation)
for c in range(1, 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 = newcounts # = num_poly[c-1]
for r in range(c): # add cell (c,r)
newcounts = defaultdict(lambda : defaultdict(int))
if CENTER_IN_SQUARE and r==0:
extra_count = 1 # B -> B+1
else:
extra_count = 2 # I -> I+1
for (state,touch_bottom,touch_top), nums in oldcounts.items():
newstate = succ1(state,r)
update((newstate,touch_bottom or r==0,touch_top), extra_count)
newstate = succ0(state,r)
if newstate:
update((newstate,touch_bottom,touch_top)) # n unchanged
oldcounts = newcounts
# now add the diagonal cell (c,c)
newcounts = defaultdict(lambda : defaultdict(int))
for (state,touch_bottom,touch_top), nums in oldcounts.items():
newstate = succ1_end(state)
update((newstate,touch_bottom,True), 1) # B -> B+1
newstate = succ0_end(state)
update((newstate,touch_bottom,touch_top))
space = accumulate_results(newcounts) # collect and accumulate the results
# 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)
print(f"total gaps = {sum_gap}, # intervals={num_intervals} "
f"avg. gap {sum_gap/num_intervals:1.2f}")
print(explanation)
OK = "Values OK."
for n in range(1,target+1):
print(end=f"{n}:{count_poly_w[n]}")
try:
if CENTER_IN_SQUARE:
if CENTER_SQUARE_OCCUPIED:
truth = A351127[4*n]
else:
truth = A351127[4*n-1]
else:
truth=A346800[n-1]
if count_poly_w[n]==truth:
print(end=". ")
else:
print(end=f" d({count_poly_w[n]-truth}) ")
OK = ""
except IndexError:
print(end=" ")
if OK:
OK = "Values OK as far as known"
print("\n"+OK)