\\ Kevin Ryde, June 2022
\\
\\ This is some PARI/GP code calculating rows of A354322, the
\\ distinct subtrees occurring in Matula-Goebel tree n.

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

\\-----------------------------------------------------------------------------
\\ Individual Row

\\ Return a vector of row n of A354322.
row(n) =
{
  my(ret=[n],
     pos=-1);  \\ offset from end of ret[]
  while(pos++ < #ret,
        ret=setunion(ret,
              [primepi(p) |p<-factor(ret[#ret-pos])[,1]]));
  ret;
}

{
  print("A354322 rows");
  for(n=1,12,
     printf("  n=%2d  %s\n",
            n, row(n)));
}

\\ "pos" workds down through ret[] from its end-most term to the start.
\\ ret[#ret-pos] is the largest subtree number not yet factored.
\\ Its child subtrees, which are primepi() indices of its distinct prime
\\ factors, merge into ret[].
\\ These subtrees are smaller numbers than their parent, so changes in ret[]
\\ are all before "#ret-pos".
\\ When #ret-pos goes below the start of ret[], all of ret[] has been
\\ factored and merged.
\\
\\ ret[] is a combination queue of pending work in the terms before pos, and
\\ result terms after pos.  In all cases ret[] terms are in ascending
\\ numerical order.
\\
\\ If a subtree occurs under multiple parents then this queue approach
\\ factor()s it just once.  In practice large common subtrees are rare
\\ though.
\\
\\ The main limitation in the code is primepi() on the prime factors of n.
\\ Circa PARI version 2.13, primepi() is good within the "primelimit" table
\\ size, but slows down rapidly beyond that.



\\-----------------------------------------------------------------------------
\\ Vector of Rows

\\ Return rows 1..maxn of A354322 as a vector of vectors.
\\ The return rows[] has rows[n] = row(n).
\\ Uses forfactored() which is new in PARI 2.11.x.
\\
rows_vector(maxn) =
{
  my(rows=Vec([[1]],maxn));
  forfactored(f=2,#rows,
     rows[f[1]] = concat(fold(setunion,
                              [rows[primepi(p)] |p<-f[2][,1]~]),
                         f[1]));  \\ n itself
  rows;
}

rows_vector(200) == vector(200,n,row(n))  || error();

\\ rows_vector() implements the union formula from the sequence,
\\
\\     row(n) = /      Union     row(primepi(p)) \
\\              \ p prime facor of n             /
\\              and followed by n itself
\\
\\ If a block of rows are wanted then this way is faster, since it re-uses
\\ child subtree rows already calculated.


\\-----------------------------------------------------------------------------
print("end");