# Impress Donald Knuth by counting polyominoes on the hyperbolic plane

This challenge is inspired by a talk about Schläfli symbols, etc that I gave in a Geometry seminar. While I was putting together this challenge, I saw that Donald Knuth himself was interested in (some subset of) this problem. In October 2016, he commented on a related OEIS sequence:

If [the OEIS author] is wrong about the hyperbolic {4,5} pentominoes, the next number is probably mistaken too. I don't have [time] right now to investigate further.

Successful completion of this challenge will have you investigating something that Donald Knuth might have investigated if only he had more time, and will result in new additions (and perhaps a rare correction) to the On-Line Encyclopedia of Integer Sequences.

# Challenge

This challenge will have you create a function that counts "free polyforms" on the hyperbolic plane. In particular, you will write a function that takes three positive integer parameters p, q, and n and counts the number of $$\n\$$-cell "free polyforms" on the tiling of the hyperbolic plane given by the Schläfli symbol $$\\{p,q\}\$$.

Shortest code wins.

# Definitions

The Schläfli symbol $$\\{p,q\}\$$ describes a tiling of the hyperbolic plane by $$\p\$$-gons, where each vertex touches exactly $$\q\$$ of the polygons. As an example, see the Wikipedia page for the $$\\{4,5\}\$$ tiling that Donald references above.

A free polyform is a collection of regular polygons that meet at their edges, counted up to rotation and reflection.

# Input

You can assume that the values of p and q which define the tiling indeed describe an actual tiling of the hyperbolic plane. This means that $$\p \geq 3\$$, and

• when $$\p = 3\$$, $$\q \geq 7\$$,
• when $$\p = 4\$$, $$\q \geq 5\$$,
• when $$\p = 5\$$, $$\q \geq 4\$$,
• when $$\p = 6\$$, $$\q \geq 4\$$, and
• when $$\p \geq 7\$$, $$\q \geq 3\$$.

# Data

OEIS sequence A119611 claims that f(4,5,n) = A119611(n), but Donald Knuth disputes the reasoning for the value of $$\A119611(5)\$$. (When I counted by hand, I got Knuth's answer, and I've included it in the table below.)

$$\p\$$ $$\q\$$ $$\n\$$ $$\f(p,q,n)\$$
3 7 1 1
3 7 2 1
3 7 3 1
3 7 4 3
3 7 5 4
3 7 6 12
3 7 7 27
3 9 8 82
4 5 3 2
4 5 4 5
4 5 5 16
6 4 3 3
7 3 1 1
7 3 2 1
7 3 3 3
8 3 3 4
9 3 3 4

Note: these values are computed by hand, so let me know if you suspect any mistakes.

# Final notes

The output of this program will result in quite a lot of new, interesting sequences for the OEIS. You are of course free to author any such sequences—but if you're not interested, I'll add the values you compute to the Encylopedia with a link to your answer.

• Is this code-golf or code-challenge?
– qwr
Feb 26 '20 at 4:20
• It's a code-golf challenge: the shortest code wins. But I expect that the shortest solution might be quite long, as in this answer to a previous code-golf challenge. Feb 26 '20 at 5:06
• Can you provide a reference implementation or some pseudocode? Feb 27 '20 at 21:44
• I don't know this stuff so I probably won't be able to do it. Can you maybe tell how you did it by hand? Feb 27 '20 at 22:06
• Monomino. Do do do-do-do. ♫♪ Mar 3 '20 at 22:09

# GAP and its kbmag package, 711 682 658 bytes

Note that the kbmag package consists not only of GAP code, it contains C programs that have to be compiled (see the package's README file).

LoadPackage("kbmag");I:=function(p,q,n)local F,H,R,r,s,x,c;F:=FreeGroup(2);s:=F.1;r:=F.2;R:=KBMAGRewritingSystem(F/[s^2,r^p,(s*r)^q]);AutomaticStructure(R);H:=SubgroupOfKBMAGRewritingSystem(R,[r]);AutomaticStructureOnCosets(R,H);x:=w->ReducedCosetRepresentative(R,H,w);c:=function(n,U,S,P)local N,Q,Z;if n=0 then Z:=Set(U,t->Set(U,p->(p/t)));return 1/Size(SetX(Union(Z,Set(Z,Q->Set(Q,q->(MappedWord(q,[s,r],[s,r^-1]))))),[1..p],{Q,i}->Set(Q,q->x(q*r^i))));fi;if P=[]then return 0;fi;N:=P;Q:=P{[2..Size(P)]};Z:=Filtered(Set([1..p],i->x(s*r^i*N)),w->not w in S);return c(n,U,S,Q)+c(n-1,Union(U,[N]),Union(S,Z),Union(Q,Z));end;return c(n,[],[r/r],[r/r]);end;


This is the result of removing indentation and newlines from this version, and some inlining:

LoadPackage("kbmag");
I:=function(p,q,n)
local F,G,H,R,r,s,x,c;
F:=FreeGroup(2);
s:=F.1;r:=F.2;
G:=F/[s^2,r^p,(s*r)^q];
R:=KBMAGRewritingSystem(G);
AutomaticStructure(R);
H:=SubgroupOfKBMAGRewritingSystem(R,[r]);
AutomaticStructureOnCosets(R,H);
x:=w->ReducedCosetRepresentative(R,H,w);
c:=function(n,U,S,P)
local N,Q,Z;
if n=0 then
Z:=Set(U,t->Set(U,p->(p/t)));
Z:=Union(Z,Set(Z,Q->Set(Q,q->(MappedWord(q,[s,r],[s,r^-1])))));
Z:=SetX(Z,[1..p],{Q,i}->Set(Q,q->x(q*r^i)));
return 1/Size(Z);
fi;
if P=[]then return 0;fi;
N:=P;Q:=P{[2..Size(P)]};
Z:=Filtered(Set([1..p],i->x(s*r^i*N)),w->not w in S);
return c(n,U,S,Q)+c(n-1,Union(U,[N]),Union(S,Z),Union(Q,Z));
end;
return c(n,[],[r/r],[r/r]);
end;


If the line containing {Q,i}-> doesn't work, your GAP is too old. You can then replace that line with:

Z:=SetX(Z,[1..p],function(Q,i)return Set(Q,q->x(q*r^i));end);


Several of the Set operations could be slightly faster List operations (the improved version at least uses that it is a set for even more golfing and a little speed compensation), but that would cost one byte each time.

And yes, Knuth's and your result is confirmed:

gap> Read("i.gap");
─────────────────────────────────────────────────────────────────────────────
by Derek Holt (https://homepages.warwick.ac.uk/staff/D.F.Holt/).
Homepage: https://gap-packages.github.io/kbmag
─────────────────────────────────────────────────────────────────────────────
gap> I(4,5,5);
16
gap> I(4,5,6);
55
gap> I(4,5,7);
224
gap> I(4,5,8);
978
gap> I(4,5,9);
4507
gap> I(4,5,10);
21430


The $$\n=7\$$ computation already takes several minutes. My computations also agree with the other results in the table.

• Pretty awesome!
– user9207
Feb 28 '20 at 18:43
• I'm blown away by this. Even if someone writes something faster for the current bounty, I'll make another one to give to you. Feb 28 '20 at 20:42
• This appears to work for polyhedra and tilings of the plane too? Feb 28 '20 at 20:53
• @PeterKagey Yes, since all the groups are automatic. Feb 28 '20 at 21:34

# bc, 673 670 bytes

-3 bytes thanks to @ceilingcat

define y(){t[z=u[x[++n]=e]/p]=1;for(o=p;o--;){if(!u[a=z*p+o]){v(u[u[a]=c=s]=a);s+=p;v(c)}}}define j(e,i){return(e%p+p+i)%p+e/p*p}define v(e){auto l;for(l=j(e,b=1);y=u[e];++b){e=j(y,-1)};for(;y=u[l];++b){l=j(y,1)};if(b==q)v(u[u[l]=e]=l)}p=read();q=read();g=read();u[e=s=2*p]=p;y();e=p;for(y();e==p;){while(n-1){if(t[e/p]*!t[u[e]/p]){c=y();for(o=6;o-=2;){for(d=s;d--;){i[a=w=2]=f=u[d];if(t[d/p]*t[f/p]){i[b=1]=d;m[d/p]=m[f/p]=++r;for(z=0;w-z++;){for(f=i[z];f=j(f,o-3)!=i[z];a*=2){if(t[y=u[f]/p]){b+=a;if(m[y]<r){i[++w]=u[f];m[y]=r}}}};if(c<b)c=b}}};for(k[h]=c;k[o++]-c;){};if(o>h){++h;if(n<g){e=p;break};++l};t[u[e]/p]=0;--n};while(++e>s){t[u[e=x[n--]]/p]=0}}};l;if(g<3)1



Try it online!

This is a golfed version of the C code I'm about to post. Here, you have to supply the parameters (p, q and n) on standard input when bc is running. The code spits out lots of numbers while it calculates (because statements such as i+=1 are silent, but ++i is a value and is printed) -- the answer is the final number.

n<3 has to be handled as a special case. Infinite planar grids are handled OK (which gives us lots of test cases). Finite grids (the Platonic solids, but also p=2 and/or q=2) are not handled reliably.

Here we are counting "free" polyominoes, unique up to translation, rotation and reflection. If you had asked the same question but disallowing reflection, I could have saved a few bytes and even beaten Christian's excellent answer as it stands. (I don't know if that answer would be shorter for the not-free question.) Here's the not-free code: Try it online! (655 652 bytes, answering a different question).

Here's the main code unwrapped:

 define y(){
t[z=u[x[++n]=e]/p]=1;
for(o=p; o--;){
if(!u[a=z*p+o]){
v(u[u[a]=c=s]=a);
s+=p;
v(c)
}
}
}
define j(e,i){
return(e%p+p+i)%p+e/p*p
}
define v(e){
auto l;
for(l=j(e,b=1); y=u[e]; ++b){
e=j(y,-1)
};
for(; y=u[l]; ++b){
l=j(y,1)
};
if(b==q)v(u[u[l]=e]=l)
}
u[e=s=2*p]=p;
y();
e=p;
for(y(); e==p;){
while(n-1){
if(t[e/p]*!t[u[e]/p]){
c=y();
for(o=6; o-=2;){
for(d=s; d--;){
i[a=w=2]=f=u[d];
if(t[d/p]*t[f/p]){
i[b=1]=d;
m[d/p]=m[f/p]=++r;
for(z=0; w-z++;){
for(f=i[z]; f=j(f,o-3)!=i[z]; a*=2){
if(t[y=u[f]/p]){
b+=a;
if(m[y]<r){
i[++w]=u[f];
m[y]=r
}
}
}
};
if(c<b)c=b
}
}
};
for(k[h]=c; k[o++]-c;){
};
if(o>h){
++h;
if(n<g){
e=p;
break
};
++l
};
t[u[e]/p]=0;
--n
};
while(++e>s){
t[u[e=x[n--]]/p]=0
}
}
};
l;
if(g<3)1


# C (gcc), 753745 705 bytes

-8 bytes thanks to @ceilingcat

-40 bytes, @ceilingcat strikes again!

#define X
#define S o=p;F[z=N[O[++n]=e]/p]=1 W o--)N[a=z*p+o]||v(c,s+=p,v(N[N[a]=c=s]=a))
#define W;while(
#define T atoi(*++V)
long C X,a,b,c,p,q,o,N X,F X,s,O X,A,n,P,e,m X,Q X,d,f,j,w,z,r;I(e,i){return(e%p+p+i)%p+e/p*p;}v(e,l){l=I(e,b=1)W N[e])++b,e=I(N[e],-1)W N[l])l=I(N[l],1),++b;b-q||v(N[N[l]=e]=l);}main(M,V)int**V;{p=T;q=T;M=T;N[e=s=2*p]=S;e=S;G:e=p W~-n){if(F[e/p]&!F[N[e]/p]){S;c=0;o=4 W o){d=s W d--)if(F[d/p]&F[(Q[a=w=2]=j=N[d])/p]){Q[b=1]=d;m[d/p]=m[j/p]=++r;z=0 W w-z++){f=Q[z]W(f=I(f,o-3))-Q[z])j=N[f],F[j/p]&&m[b+=a,j/p]<r?Q[++w]=j,m[j/p]=r:0,a+=a;}c=c<b?b:c;}o-=2;}C[P]=c W C[o++]-c);if(o>P){++P;if(n<M)goto G;++A;}F[N[e]/p]=0;--n;}W++e>s)F[N[e=O[n--]]/p]=0;}printf("%d",M<3?:A);}


Try it online!

Parameters p, q and n are supplied as command-line arguments.

The array C contains a characteristic value for each polyomino found, of the specified size and all smaller (down to n=2). For p less than 7 or 8 or so, this is probably the limit on how far you can go up a sequence. For larger p, you might well run out of edges, for example in the array N. There is also a constraint that the characteristic values need (p-1)*n bits, and here we are using 63ish bits. If we stick to unambitious values of p, q and n, I could save a handful of bytes with int and 9999. If you exceed any of the limits, there is no error trapping.

Here's the code unwrapped, but it's still fairly dense:

#define X
#define S o=p;F[z=N[O[++n]=e]/p]=1;while(o--)N[a=z*p+o]||v(c,s+=p,v(N[N[a]=c=s]=a))
#define T atoi(*++V)
long C X,a,b,c,p,q,o,N X,F X,s,O X,A,n,P,e,m X,Q X,d,f,j,w,z,r;
I(e,i){
return(e%p+p+i)%p+e/p*p;
}
v(e,l){
l=I(e,b=1);
while(N[e])++b,e=I(N[e],-1);
while(N[l])l=I(N[l],1),++b;
b-q||v(N[N[l]=e]=l);
}
main(M,V)int**V;
{
p=T;
q=T;
M=T;
N[e=s=2*p]=S;
e=S;
G:e=p;
while(~-n){
if(F[e/p]&!F[N[e]/p]){
S;
c=0;
o=4;
while(o){
d=s;
while(d--)if(F[d/p]&F[(Q[a=w=2]=j=N[d])/p]){
Q[b=1]=d;
m[d/p]=m[j/p]=++r;
z=0;
while(w-z++){
f=Q[z];
while((f=I(f,o-3))-Q[z])
j=N[f],F[j/p]&&m[b+=a,j/p]<r?Q[++w]=j,m[j/p]=r:0,a+=a;
}
c=c<b?b:c;
}
o-=2;
}
C[P]=c;
while(C[o++]-c);
if(o>P){
++P;
if(n<M)goto G;
++A;
}
F[N[e]/p]=0;
--n;
}
;
while(++e>s)F[N[e=O[n--]]/p]=0;
}
printf("%d",M<3?:A);
}


This post seems rather long already, so I'll link to TIO with a copy of ungolfed code with comments. The golfed code uses basically the same approach: Try it online!

In the ungolfed code, if you supply p, q, and n, it will try the size-n problem, then size-(n+1) and so on until something breaks. If limits are exceeded, the program is meant to terminate with a message. You can supply a fourth parameter to cap n.

The ungolfed code should handle finite grids (Platonic solids, or p=2 and/or q=2).

If you want to look at the different question of counts up to translation and rotation only, compile the ungolfed code with -DNOTFREE.

The golfed code can be checked against existing planar values as well as the values in the question. For example:

p=4,q=4,    n=11: a000105(11)=17073
p=3,q=99999,n=13: a000207(13)=24834
p=6,q=3,    n=10: a000228(10)=30490
p=3,q=6,    n=14: a000577(14)=26166
p=4,q=99999,n=10: a005036(10)=32721
p=4,q=5,    n=10: a119611(10)=21430


Most of these run within TIO's timeout of 60s. In other sequences, such as a330659 {3,7}, the golfed code can replicate and extend the current known values in OEIS:

p=3,q=7,    n=10: a330659(10)=  637
p=3,q=7,    n=11: a330659(11)= 1870
p=3,q=7,    n=12: a330659(12)= 5797
p=3,q=7,    n=13: a330659(13)=17866
p=3,q=7,    n=14: a330659(14)=56237


The ungolfed code has more checks, uses more memory and runs quicker (because, for example, it uses a hash table rather than a flat array for storing characteristics), so I'll add to the OEIS using that code. For example, extending a330659 to n=19 takes around 1000 seconds:

p=3,q=7,    n=14: a330659(14)=56237
p=3,q=7,    n=15: a330659(15)=177573
p=3,q=7,    n=16: a330659(16)=566904
p=3,q=7,    n=17: a330659(17)=1818527
p=3,q=7,    n=18: a330659(18)=5874180
p=3,q=7,    n=19: a330659(19)=19065038

• Suggest e+i instead of e%p+p+i Feb 24 at 6:30