\\ Kevin Ryde, July 2022
\\
\\ Usage: gp <a064002.gp
\\
\\ This is some PARI/GP code to calculate terms of A064002.
\\ a(n) is an individual term by descents, and suits large or small n.
\\ a_vector(len) is initial terms using just additions.

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

\\-----------------------------------------------------------------------------
\\ Individual Term

{
  \\ return a 2-element vector [i,j], for n>=2
  my(pair(n) =
     my(r, s=sqrtint((n-1)<<1,&r));
     if(r>s,  r-=s;s++, r+=s);
     [r>>1, s]);

  a(n) = if(n<=1, n,
            my(i,j); [i,j]=pair(n);
            self()(i) + self()(j) + 1);
}
{
  print1("A064002 = ");
  for(n=1,12, print1(a(n),", "));
  print("...");
}

\\ a(n) is a recursive implementation of the defining formula
\\
\\     a(n) = a(i) + a(j) + 1
\\
\\ pair(n) returns the relevant pair i,j for use at n
\\
\\     i = A002260(n-1)
\\     j = A002024(n-1)
\\
\\ calculated together by a sqrtint() and its "remainder".
\\
\\ Recursion is easy, and it measures faster than an explicit List() of
\\ pending descents.
\\
\\ Peak recursion depth is modest since the smallest n requiring depth k is
\\ A108225(k) and that becomes very large very quickly.  Eg. depth 30 is a
\\ roughly 10 megabyte number.


\\-----------------------------------------------------------------------------
\\ Vector of Terms

\\ A vector of initial terms from the sequence can follow the defining
\\ formula above, but iterating through i,j pairs
\\
\\     all 1 <= i <= j
\\     in order ascending j and within j ascending i
\\
\\ The first pair i=1,j=1 is for n=2.

a_vector(len) =
{
  my(v=vector(len,i,1),
     i=1,
     j=1);  \\ 1 <= j <= i
  for(n=2,len,
     v[n] = v[i] + v[j] + 1;
     if(j++ > i, i++;j=1));
  v;
}

{
  \\ check a_vector() against a()
  my(v=a_vector(500));
  #v == 500                || error();
  vector(#v,n, a(n)) == v  || error();
}

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