\\ Kevin Ryde, September 2022
\\
\\ Usage: gp <a356918.gp
\\
\\ This is some PARI/GP code to calculate terms of triangle A356918,
\\ Colijn and Plazzotta's distance metric d_1 between two rooted binary
\\ trees (every vertex 0 or 2 children), in their tree enumeration.
\\
\\ Requires PARI/GP 2.13 or higher for sqrtint() remainder and
\\ call-by-reference.

default(recover,0);
default(strictargs,1);

{
  \\ children(n), for n>=2, returns a vector [x,y] of the child subtrees of n,
  \\ x = A002024(n-1) and y = A002260(n-1)
  my(children(n) =
    my(r, s=sqrtint((n-1)<<1,&r));
    if(r>s,  r-=s;s++, r+=s);
    [s,r>>1]);

  \\ map_count_subtrees(~m,n,add) takes map m of subtree number => count
  \\ and applies +add at n and each subtree number within in n
  \\ (with multiplicity).
  my(map_count_subtrees(~m,n,add) =
     my(count=0); mapisdefined(m,n,&count);
     mapput(m,n,count+add);
     if(n>1, foreach(children(n),c, self()(~m,c,add))));

\\ A356918
T(n,k) =
  my(m=Map());
  map_count_subtrees(~m,n,1);
  map_count_subtrees(~m,k,-1);
  vecsum(abs(Mat(m)[,2]));
}

{
  print("A356918 rows:");
  for(n=1,9,
     printf("  n=%2d: ",n);
     for(k=1,n,
        printf(" %2d,", T(n,k)));
     print());
  print("  ...");
}

\\ The strategy in T(n,k) is to traverse the trees to find their subtree
\\ numbers and count occurrences in a Map().
\\ Subtree numbers are sparse, hence a map instead of a counts vector.
\\
\\ Subtrees in n count +1 each, and subtrees in k count -1 each, so that
\\ abs() of each net count value is the difference in number of occurrences
\\ in n and k.
\\
\\ Trees are traversed recursively.  Peak recursion depth is modest since
\\ the smallest n of height h is A108225(h) and that becomes very large very
\\ quickly.  Eg. height 30 is a roughly 10 megabyte number.
\\
\\ An alternative to recursion would be an explicit stack of pending work in
\\ a List().  But that measures barely different in speed.

print("end");