\\ 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");