\\ Author: Andrew Howroyd, Feb 2021.
\\ The following PARI code is largely based on the following reference.
\\ William G. Brown, Enumeration of Triangulations of the Disk, Proc. Lond. Math. Soc. s3-14 (1964) 746-768.

\\ As an introduction, it is first suggested to reads Tutte's paper.
\\ William T. Tutte, A census of planar triangulations, Canad. J. Math. 14 (1962), 21-38. See Eq. 5.12.
\\ Tutte's paper does not deal with symmetry and also deals with 3-connected triangulations.
\\ The triangulations enumerated here are 2-connected.

\\ Sequences covered by this script are the following arrays including columns, rows and other related sequences.
\\ A146305: rooted triangulations.
\\ A262586: up to orientation preserving isomorphism.
\\ A169809: unrooted with reflection symmetry.
\\ A169808: unrooted. A169808 = (A169809 + A262586)/2.

\\ Helper to make bivariate g.f. from a function.
MakeSquareBgf(fun, N, M, x, y)={sum(n=0, N, x^n*(O(y*y^M) + sum(m=0, M, y^m*fun(n,m)))) + O(x*x^N)}
\\ Helper to make matrix from bivariate g.f.
BgfToArray(gf, N, M)={matrix(N+1, M+1, n, m, polcoef(polcoef(gf, n-1, x), m-1, y))}


\\ Rooted disk triangulations.
\\ Triangulations are biconnected with n internal nodes and m+3 external nodes.
D(n,m)={2*(2*m+3)!*(4*n+2*m+1)!/(m!*(m+2)!*n!*(3*n+2*m+3)!)}

\\ Triangulations with Rotational Symmetry (Section II in Brown)

\\ Triangulations invariant under rotation of 1/d with d >= 3 about a central vertex.
\\ There are d*s+1 internal nodes and d*(p+1) external nodes.
\\ Independent of d.
\\ See 8.12 in Brown.
Er(s,p)={(2*p+2)!*(4*s+2*p+1)!/(p!*(p+1)!*s!*(3*s+2*p+2)!)}

\\ Triangulations invariant under rotation of 1/3 about a central triangle.
\\ There are 3*s internal nodes and 3*(p+1) external nodes.
\\ See 8.11 in Brown.
E3(s,p)={(2*p+1)!*(4*s+2*p)!/(p!*p!*s!*(3*s+2*p+1)!)}

\\ Triangulations invariant under rotation of 1/2.
\\ See 8.10 in Brown.
E2(s,j,p)=2*(2*p)!*(4*s+2*p+2*j-1)!/(p!*(p-1)!*s!*(3*s+2*p+2*j)!)

\\ Oriented triangulations.
\\ See 6.3 in Brown.
OrientedTriangs(n,m)={(D(n,m) + if(m%2==1, E2(n\2, n%2, (m+1)/2)) + if(gcd(m,n)%3==0, 2*E3(n/3, m/3)) + sumdiv(gcd(m+3,n-1), d, if(d>2, eulerphi(d)*Er((n-1)/d, (m+3)/d-1))))/(m+3)}


\\ Triangulations with Reflection Symmetry (Section III in Brown)

\\ J_0 function (sequence A002712 as g.f.)
\\ See 13.10: Satisfies J = 1 + x*J + x^2*J*(1 + x*J/2)*(J^2 - D(x^2,0)).
\\ Compute by iteratively growing precision.
Jgf(n)={my(q=Ser(vector(n+1, i, if(i%2, D(i\2,0)))), p=1+O(x)); for(n=1, n, p = 1 + x*p + x^2*p*(1 + x*p/2)*(p^2 - q)); p}

\\ Formula for K(x,y) + L(x,y) = (alpha(x,y) + beta(x,y))/gamma(x,y).
\\ See 13.1, 13.2, 13.3, 13.4.
\\ This has been rearranged to reduce number of ops on D and then J.
\\ Note that D here is D bar in the reference, where D bar = Dbar(x^2,y^2) and Dbar(x,y) = 1 + y*D(x,y) [see 12.4]
KPlusLRatFunc(D,J,x,y)={
  ((y^2 - x^2)*x*J + (-2*y^2*(x+y) + (x^2 + y^2 + 2*x*y*(1-y))*x*J + y*(2*x + 2*y + x*y*J)*x^2*J^2 -y^4*(2 + x*J)*D)*D)
    / (y*(x^2 - y^2 + x*(x^2-y^2)*J) + x*y^3*(2 + J*x)*D)
}


\\ Triangulations with reflective symmetry.
\\ Computes (K(x,y) + L(x,y))/2 to required precision.
\\ Returns bivariate generating function.
\\ The precision of x is N, precision of y varies and reduces by coefficient of x from N + M down to M at end.
AchiralTriangsXY(N, M)={
   KPlusLRatFunc(1 + y^2*MakeSquareBgf(D, N\2, (N+M)\2-1, x^2, y^2), Jgf(N), x, y)/2 + O(x*x^N)
}

\\ Debugging helper - shows precision of a bivariate g.f. by primary coefficient.
SubSerPrec(p)=[serprec(q,y) | q<-Vec(p)]

\\ Sequences for achiral triangulations
A169809Array(N,M=N)={BgfToArray(AchiralTriangsXY(N, M), N, M)}
A169809Triang(N)={my(v=Vec(AchiralTriangsXY(N, 0))); vector(#v, n, vector(n, k, polcoef(v[k], n-k)))}
A169809ColSeq(N,k)={Vec(polcoef(AchiralTriangsXY(N,k),k,y))}
A169809RowSeq(N,k)={Vec(polcoef(AchiralTriangsXY(k,N),k,x) + O(y*y^N))}

\\ Sequences for unrooted triangulations
A169808Array(N, M=N)={(A169809Array(N,M) + matrix(N+1, M+1, i, j, OrientedTriangs(i-1,j-1)))/2}
A169808Triang(N)={(A169809Triang(N) + vector(1+N, n, vector(n, k, OrientedTriangs(k-1,n-k))))/2}
A169808ColSeq(N,k)={(A169809ColSeq(N,k) + vector(N+1, i, OrientedTriangs(i-1, k)))/2}
A169808RowSeq(N,k)={(A169809RowSeq(N,k) + vector(N+1, i, OrientedTriangs(k, i-1)))/2}

\\ Faster versions for antidiagonal sum that avoids bivariate g.f.
A169809AntidiagonalSums(N)={
   my(p=sum(n=0, N\2, x^(2*n)*sum(k=0, n, D(n-k,k))) + O(x*x^(2*N\2)));
   Vec(KPlusLRatFunc(1 + x^2*p, Jgf(N+1), x, x)/2 + O(x*x^N))
}

A169808AntidiagonalSums(N)={(A169809AntidiagonalSums(N) + vector(N+1, n, sum(k=0, n-1, OrientedTriangs(k, n-1-k))) )/2}


\\ End