\\ Kevin Ryde, October 2022
\\
\\ This is some PARI/GP code to calculate rows of A357701 which is the
\\ pre-order depths vector of the rooted binary tree with Colijn-Plazzotta
\\ tree number n.
\\ 
\\ Requires PARI/GP 2.13 for sqrtint() remainder.


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]);

  \\ traverse() returns the depths vector of subtree n
  \\ with the root taken as "depth" (an integer)
  my(traverse(n,depth) =
     if(n==1, [depth],  \\ singleton
        my(c=children(n));
        concat([depth,  \\ pre-order
                self()(c[1],depth+1),
                self()(c[2],depth+1)])));

row(n) = traverse(n,0);
}

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

\\ traverse() is a recursive pre-order traversal of the tree.
\\ Peak recursion depth is peak vertex depth (tree height) and that's modest
\\ since the smallest n of height h is n = A108225(h) which is very large
\\ very quickly.
\\
\\ Subtree vectors returned by traverse() are concatenated which is quite a
\\ lot of vector copying when the tree is large.  This takes a measurable
\\ time (by comparing returning just a vertex count say).
\\ 
\\ But speed measures roughly the same for listput() to a call-by-reference
\\ target, or even just storing into a pre-sized vector in a global
\\ variable.


\\-----------------------------------------------------------------------------
\\ Vector of Rows
\\
\\ Initial rows of the sequence can be calculated by row copying and
\\ incrementing previous rows looping over x and y.  Each incr() is the
\\ "row+1" in the sequence formula
\\ 
\\     row(n) = 0, row(x)+1, row(y)+1

{
  \\ increment each term in vector v
  my(incr(v) = apply(d->d+1, v));

\\ return a vector of vectors
rows_vector(len) =
  my(ret=vector(len));
  if(len,
     ret[1] = [0];  \\ singleton
     my(n=1);       \\ last stored in ret[]
     for(x=1,oo,
        my(xrow=incr(ret[x]));
        for(y=1,x,
           if(n++ > len, break(2));
           ret[n] = concat([0, xrow, incr(ret[y])]))));
  ret;
}

{
  my(v=rows_vector(100));
  for(n=1,100,
     v[n] == row(n) || error());
}

\\ As noted in the sequence comments, rows are in lexicographically
\\ increasing order, which is lex().
{
  my(rows=rows_vector(100));
  for(n=2,#rows,
     lex(rows[n], rows[n-1]) > 0  || error());
}


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