\\ Kevin Ryde, September 2021
\\
\\ Usage; gp <a061773.gp
\\
\\ This is some PARI/GP code for calculating terms of A061773, being
\\ Matula-Goebel tree numbers sorted by number of vertices then numerically.
\\
\\ Rooted trees of n vertices are a new subtree under the root of a tree of
\\ < n vertices, or a new root above any tree of n-1 vertices.  This is a
\\ generic method and suits Matula-Goebel numbers since adding a subtree is
\\ a multiply and adding a root above is prime().
\\
\\ prime(t) of large t limits how much of the sequence can be generated.
\\ Initialize PARI "primelimit" to at or above the largest Matula-Goebel
\\ number to be generated (which is A005518(n)) to ensure fast prime().
\\ 
\\ On a 32-bit system circa PARI 2.13, n=13 vertices is the largest feasible
\\ size and is fast, inclusive of primelimit=997525853 startup, with enough
\\ memory.  n=14 would want prime(1 billion) = 22 billion which in 32 bits
\\ is beyond the maximum primelimit = 4 billion and so suspect more or less
\\ tests primes one-by-one from there to 22 billion.


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


\\ Return a vector of length "len" where ret[n] = A000081(n) = number of
\\ rooted trees of n vertices.
num_trees_vector(len) =
{
  my(ret=vector(len,i,1),
     c=vector(max(0,len-1)));   \\ c[n] = A209397(n)
  for(n=1,len-1,
     c[n] = sumdiv(n,d, ret[d]*d);
     ret[n+1] = sum(i=1,n, ret[i]*c[n-i+1]) / n);
  ret;
}


\\ A061773_rows(len) returns a "rows" vector of length "len" where rows[n]
\\ is a vector of the terms from row n of A061773, being the trees of n
\\ vertices.  For example,
\\ 
\\     A061773_rows(4) == [ [1],           \\ 1 vertex (singleton)
\\                          [2],           \\ 2 vertices (path-2)
\\                          [3,4],         \\ 3 vertices (path and star)
\\                          [5,6,7,8] ]    \\ 4 vertices
\\
\\ Trees of n vertices are created by two cases.  Firstly to make trees
\\ with 2 or more subtrees under the root, a tree P of p>=2 many vertices
\\ and just 1 subtree under the root is merged onto another tree S,
\\
\\         S   P                R
\\       / |   |     ==>      / | \
\\      X  Y   Z             X  Y   Z 
\\     /|  |\  |\           /|  |\  |\
\\
\\     R = S*P multiply to merge roots S,P of two trees 
\\     s + p = n+1 vertices, so total n after merge roots
\\
\\ In Matula-Goebel numbers, merging the roots S and P like this is as
\\ simple as multiply S*P.
\\
\\ To avoid creating duplicate trees, have new subtree Z >= X,Y, ie. >=
\\ all subtrees of S.  Comparison ">=" is by a "bigness" index (an integer)
\\ which is assigned to tree P when P is created and which represents the
\\ bigness of its subtree Z (sole subtree of P).
\\
\\ Planted trees P of n vertices (only 1 subtree under the root) are created
\\ by taking any tree Z of n-1 vertices and adding a root above.  In
\\ Matula-Goebel numbers, this means P = prime(Z).  (Repeated nested
\\ prime(prime(prime(...))) quickly become large numbers though.)
\\
\\ rows[] is a vector of vectors where rows[n] = vector of the Matula-Goebel
\\ numbers of trees of n vertices, sorted in ascending order of biggest
\\ child subtree.  (The final return sorts them to numerical order.)
\\
\\ bigs[] is a vector of vectors where bigs[n][i] is the bigness index for
\\ the biggest child subtree of tree rows[n][i].  bigs[n] is in ascending
\\ order.
\\
\\ Ascending biggest subtree means rows[n] starts with the non-planted trees
\\ (2 or more child subtrees at the root), followed by the planted trees.
\\ In the non-planted creation, loop "i" gets planted trees P by knowing
\\ they're the end-most num_trees[p-1] many entries of rows[p].


A061773_rows(maxn) =
{
  my(num_trees=num_trees_vector(maxn),
     rows=vector(maxn,n, if(n==1,[1], vector(num_trees[n]))),
     bigs=apply(Vecsmall,rows),
     bigidx=1);
  for(n=2,maxn,
     my(upto=0);  \\ where to store in rows[n] vector

     \\ non-planted trees
     for(p=2,n-1,       \\ number of vertices planted part
        my(s=n-p+1);    \\ number of vertices other part, p+s == n+1
        for(i=num_trees[p]-num_trees[p-1]+1, num_trees[p], \\ planted range
           my(P     = rows[p][i],   \\ planted tree P
              P_big = bigs[p][i]);  \\ bigness of its Z subtree
           for(j=1,#rows[s],
              P_big >= bigs[s][j] || break;
              rows[n][upto++] = P*rows[s][j];
              bigs[n][upto]   = P_big)));

     \\ planted trees
     for(i=1,num_trees[n-1],
        rows[n][upto++] = prime(rows[n-1][i]);
        bigs[n][upto]   = bigidx++));

  apply(vecsort,rows);
}

{
  \\ print some sample rows, abbreviated where necessary
  my(rows=A061773_rows(10),
     maxwidth=56); \\ characters
  print("A061773_rows()");
  for(n=1,#rows,
     printf("  n=%2d:  ",n);
     my(w=#Str(rows[n][#rows[n]]));  \\ width of last term
     for(i=1,#rows[n],
        print1(rows[n][i]);
        w += #Str(rows[n][i])+1;
        if(i<#rows[n],print1(","));
        if(i<#rows[n]-1 && w>maxwidth,  \\ too wide, abbreviate
           print1("..,",rows[n][#rows[n]]); break));
     print());

  print();
  print1("row widths ");
  for(n=1,#rows,
     print1(#rows[n],","));
  print("... (A000081)");
}

print("end");